Abstract

Novel coronavirus (COVID-19) has been spreading and wreaking havoc globally, despite massive efforts by the government and World Health Organization (WHO). Consideration of partially recovered carriers is hypothesized to play a leading role in the persistence of the disease and its introduction to new areas. A model for transmission of COVID-19 by symptomless partially recovered carriers is proposed and analysed. It is shown that key parameters can be identified such that below a threshold level, built on these parameters, the epidemic tends towards extinction, while above another threshold, it tends towards a nontrivial epidemic state. Moreover, optimal control analysis of the model, using Pontryagin’s maximum principle, is performed. The optimal controls are characterized in terms of the optimality system and solved numerically for several scenarios. Numerical simulations and sensitivity analysis of the basic reproduction number, , indicate that the disease is mainly driven by parameters involving the partially recovered carriers rather than symptomatic ones. Moreover, optimal control analysis of the model, using Pontryagin’s maximum principle, is performed. The optimal controls were characterized in terms of the optimality system and solved numerically for several scenarios. Numerical simulations were explored to illustrate our theoretical findings, scenarios were built, and the model predicted that social distancing and treatment of the symptomatic will slow down the epidemic curve and reduce mortality of COVID-19 given that there is an average adherence to social distancing and effective treatment are administered.

1. Introduction

Coronavirus is one of the major pathogens that primarily targets the human respiratory system. The previous outbreak of coronaviruses include the severe acute respiratory syndrome (SARS), which occurred in 2003 in Mainland China [1], the Middle East respiratory Syndrome (MERS) in 2012 in Saudi Arabia [2], and the MERS outbreak in 2015 in South Korea [3]. In late December 2019, a cluster of patients was admitted to hospitals with an initial diagnosis of pneumonia of an unknown etiology. These patients were etiologically linked to seafood and wet animal wholesale market in Wuhan, the capital city of the Hubei province, China [4]. Early reports predicted the onset of a potential coronavirus outbreak and estimated the reproduction number to be significantly larger than 1 (ranges from 2.24 to 3.58) [5].

A report published in Nature revealed that Chinese health authorities concluded that as of February 7, 2020, there have been 31,161 people who have contracted the infection in China and more than 630 people have died of the infection [6] [http://www.nature.com/articles/d41586-020-00154]. Also, the World Health Organization (WHO) reported 51,174 confirmed cases including 15,384 severe cases and 1,666 death cases in China. At the time of preparing this manuscript, the number of cases (March 7, 2020) has passed 100,000 worldwide, reported by Euronews.

The most common symptoms at the onset of COVID-19 illness are fever, cough, and fatigue, while other symptoms include sputum production, headache, hemoptysis, diarrhea, and lymphopenia [7].

Mathematical models for dynamics of the transmission and control of COVID-19 have been developed to gain insights on the disease dynamics. Also, as recognized by the WHO [8], mathematical models, especially those who are timely, play a key role in informing evidence-based decisions by health decision and policy makers. Only a few mathematical models have so far been publicly released, to the best of our knowledge. Notably among these studies are [918].

This study has an important difference from those reported above, in that it presents a mathematical model that considers a proportion of infectious patients who recover partially. These proportions may still be virus carriers. According to a new study developed by Lan et al. [19], they followed four medical professionals ages 30 to 36 years who developed COVID-19 and were treated at Wuhan University’s Zhongnan Hospital in China between Jan. 1 and Feb. 15. All of the individuals recovered, and only one was hospitalized during the illness. The patients were treated with oseltamivir, better known under the brand name of Tamiflu, an antiviral drug. The patients were considered recovered after their symptoms. After recovery, the patients were asked to quarantine themselves at home for five days. They continued to undergo throat swabs for the coronavirus after five days for up to 13 days postrecovery. The results showed that every test between day 5 and day 13 was positive for the virus. Also, a paper published by Pappas [20] gave an instance of a Japanese patient who recovered from COVID-19 and then became ill with the disease for a second time. Pappas [20] explained further that one possibility is that the patient’s immune system did not fight off the virus completely as it began to replicate inside her lung. These findings suggested that at least a proportion of infectious individuals did not recover fully from the disease because a small concentration of the virus might still be in their system which began to replicate. Hence, these groups are classified as partially recovered carriers, as they are virus carriers. This is a new finding that has not been modelled mathematically. Therefore, the present study bridges this knowledge gap.

The remaining part of this paper is organized as follows: in Section 2, model description and formulation are presented. Detailed stability of disease-free equilibrium points is analysed in Section 3. Sections 4 presents the sensitivity analysis of the model parameters while in Section 5, we provide the optimal control problem, its formulation, optimal strategies, proof of existence, and necessary conditions for optimal control. Section 6 provides the numerical simulations of the formulated model. Finally, Section 7 gives the overall concluding remarks of the study.

2. Mathematical Formulation of the COVID-19 Model

In this section, a model for the spread of COVID-19 in human and vector population is formulated. The total human population denoted by , is partitioned into five classes, namely, the susceptible individual , the exposed individual , the infectious individual , the partially recovered carriers , and the fully recovered individuals , so that .

2.1. Assumptions of the Model

The following assumptions were made in order to formulate the equations of the model:(a)Containment rate, , is introduced into the susceptible population(b)SARS-COV-2 can persist in the body for at least two weeks after symptoms of the disease clear up [20](c)Infectious patients whose immune system fights off the virus completely with no single virus left progress to the fully recovered class(d)Infectious patients whose immune system does not fight off the virus completely (i.e., low level of the virus are still present in their system) progress to the partially recovered carrier class(e)Reverse transcription polymerase chain reaction (RT-PCR) may not be able to detect the presence of a very low level of the virus in the system of partially recovered carriers at the time that they are discharged from the hospital(f)Virus begins to replicate inside the lungs of the partially recovered carriers again [20](g)Partially recovered carriers are contagious but show mild or no symptoms(h)Individuals in the partially recovered carrier class can also recover from COVID-19, but with a significantly slow rate

Our model includes a net inflow of susceptible individuals into a region at a rate, per unit time. This parameter includes new births, immigration, and emigration. The susceptible individuals acquire COVID-19 through contact with either the infectious individuals or partially recovered carriers at rates and . Also, the susceptible population is reduced by containment rate, . The proportion of the infectious individuals become partially recovered carriers at a rate . Individuals in the partially recovered carrier class recover fully from the disease at a slower rate .

Applying the assumptions, nomenclature of parameters, and definitions of variables, the following system of ordinary differential equations is formulated:with initial conditionswhere the model parameters are nonnegative.

Remark 1. We restructure and modify the model to have more insightful information. However, it has not in any way changed the research focus.
The first four equations, i.e., (1), (2), (3), and (4), are independent of the compartments , i.e., (5). Therefore, after decoupling the equations for from models (1), (2), (3), (4), and (5), we devote analyses on the remaining equations of (1), (2), (3), and (4) which becomeswhere the disease-free equilibrium point is .

3. Stability Analysis

3.1. Local Stability of Disease-Free Equilibrium Solution

Using the next-generation operator approach of Diekmann et al., [21], the effective basic reproduction number associated with disease-free equilibrium, , and denoted by , is obtained aswhere

Clearly, the disease-free equilibrium is locally asymptotically stable if . To see this, we obtain the Jacobian matrix of systems (7), (8), (9), and (10) evaluated at :where

One of the eigenvalues is . The other three are eigenvalues of the matrixwhose characteristic equation iswhere

The Routh-Hurwitz stability criteria for a matrix is that all values of the determinant of Hurwitz matrices are positive. The determinants of Hurwitz matrices are:where and

From (22), we have or

(23) gives

Clearly, is positive in (19). Also, or , in (17), (18), and (19). Finally, we shall show that . Since is positive, then, we need to show that is also positive so that .

Hence, if . Therefore, if . Hence, by Routh-Hurwitz theorem [22], all the eigenvalues of the Jacobian matrix have negative real parts when and the disease-free equilibrium solution is locally asymptotically stable if . The following theorem summarizes the above result in what follows:

Theorem 2. The disease-free equilibrium solution is locally asymptotically stable if .

Remark 3. The biological meaning of the theorem above is that COVID-19 eradication depends on the initial number of the infectious individuals in the population. This implies that a small invasion of infectious individuals into a completely susceptible population will not lead to an outbreak of the disease.
We present in Figures 1, 2, 3, and 4 the time series solution of the susceptible, exposed, infected, and partially recovered population of COVID-19, the prevalence against time for the unforced COVID-19 model with associated wavelet spectrum, and the estimated wavelet spectrum of the first week of the year 2 and year 4.5 for the unforced COVID-19 model.

3.2. Global Stability of Disease-Free Equilibrium Solution

The global asymptotic stability of models (19), (20), (22), and (23) is explored in what follows.

Theorem 4. The disease-free equilibrium is globally asymptotically stable if and unstable if .

Proof. Consider the Lyapunov function.Its time derivative isTherefore, if with if and only if or . It follows that the largest invariant subset in is the singleton , and hence by LaSalle’s invariance principle [23], the disease-free equilibrium will be approached by all solution trajectories, and hence, the disease-free equilibrium solution is globally asymptotically stable.

Remark 5. The epidemiological implication of the above theorem is that COVID-19 can be eradicated irrespective of the initial sizes of the subpopulation of the model.

3.3. Global Stability of the Endemic Equilibrium Point

The endemic equilibrium solution satisfies the following equations:

Finding , , and in (27), (28), (29), and (30) gives

Adding (27) and (28), we get

Substituting (31) and (32) in (33) into (34) gives the following

From (35), we can obtain to be

Substituting in (31), (32), and (33) gives

Thus, the endemic equilibrium solution exists whenever .

Furthermore, we establish the global asymptotic stability of the endemic equilibrium solutions converging to the endemic equilibrium point for . We shall carry out this, by constructing a suitable Lyapunov function of Goh-Volterra type; see [24]. The result below establishes the global stability of the endemic equilibrium solution .

Theorem 6. The unique endemic equilibrium is globally asymptotically stable whenever .

Proof. Given the following equations which are satisfied by the endemic equilibrium point ,Consider the following Goh-Volterra Lyapunov functionwhere
and
with the Lyapunov time derivative obtained asUsing (38), we haveFurther simplification givesReplacing and by their values and exploiting (38), (39), (40), and (41) giveUsing (38), (39), (40), (41), (46), and (47), we haveUsing arithmetic-geometric means inequality, i.e., , where and , it follows that with if and only if , , , and .
Hence, the largest compact invariant subset of the set where is and by classical stability theorem of Lyapunov and LaSalle’s Invariance Principle, it follows that every solution in approaches for as .

Remark 7. The epidemiological implication of the above result is that COVID-19 will establish itself whenever in the population.

4. Uncertainty and Sensitivity Analysis

4.1. Local Sensitivity Analysis

In this section, we carried out sensitivity analysis of parameters of model systems (7), (8), (9), and (10) in order to determine the relative importance of the model parameters on the disease infection. To determine how best to reduce the infection, it is necessary to know the relative importance of the different factors responsible for the infections.

Sensitivity indices could be computed numerically so as to figure out parameters that have high impact on basic reproduction number and which of the parameters should be given preferential treatment by intervention strategies.

Analytically, sensitivity analysis on all parameters which account for disease dynamics is done using the Chitnis et al. [25] approach; we compute sensitivity indices of the which measures initial disease infection and allows us to measure the relative change in a state variable when a variable changes.

The normalized forward sensitivity index of a variable to a parameter is the ratio of the relative change in the variable to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives.

Definition. The normalized forward sensitivity index of a variable, , that depends differentially on a parameter, , is defined asfor .

Consequently, we derive analytical expression for the sensitivity index of aswhere denotes each parameter involved in .

Usingwhere .

We compute the sensitivity index of each parameter with respect to the , for instance

We have Table 1 which summarizes the sensitivity indices of with respect to parameters

4.1.1. Interpretation of Sensitivity Indices Obtained in Table 2

The computed sensitivity indices on with respect to the involved parameters give insights to the model system proposed. Provided all parameters remain constant, most sensitive parameters are (the contact rate between susceptible and infectious individuals, natural birth rate). The implication is that the increase in contact rate between susceptible and infectious individuals and natural birth rate have high tendency to increase in COVID-19 epidemic due to positive variation on basic reproduction number . (COVID-19-induced rate) is the only parameter that has no effect on the spread of the disease. Interestingly, improvement in slower rate of recovery () and increase in transition rate from the infectious disease class to the partially recovered () one are major parameters that reduce the disease epidemic. In the same vein, the sensitivity indexes of reveal that the increase in containment rate of susceptible individuals, recovery rate of infectious individuals, and recovery rate of partially recovered carriers contribute to the reduction in the spread of COVID-19 infection.

In light of this, we recommend that the World health Organization (WHO) should put in place health interventions and strategies that will not only improve the slower recovery rate of the partially recovered individuals but also increase the transition rate from the infectious class to partially recovered one.

4.2. Global Sensitivity Analysis

The global sensitivity analysis helps to investigate the change in the output values resulting from changes in all parameter values over the ranges in parameters ([26]. To perform the global sensitivity analysis of the reproduction number and the range of values in Table 3. We aim here to establish the most influential parameters in and also give some insightful ideas from the PRCC plot by running a sample size of 1,000. It was observed that the scatter plot in Figures 58 shows a positive relation between and parameters ,, , , , , , , , , and have a negative relation between them and . We understand that the positive relation here implies that a high rate of either of these parameters ,, , , , , , and will generate a higher transmission rate during an outbreak, while the negative relation for , , and will aid in decreasing the severity of the disease-induced death rate. Figure 8(a) presents well defined simulation results of the scatter plots in Figures 58 and that of the numerical signs in Table 2. Figure 8(a) also indicates that the most influential parameters that can make the disease spread fast are ,, , and . Here, any effort to reduce the relative relevance or impact of , , , and will reduce the spreading rate of COVID-19. Figure 9(b) shows that while (reproduction number) is decreasing with an increase in , it increases with an increase in .

5. Analysis of Optimal Control Strategies against COVID-19

There are several measures of prevention and control of infectious diseases, but in this paper, we employ the nonpharmaceutical intervention (NPI), like social distancing, , and the nonspecific treatment effort, .

5.1. Social Distancing

This is a nonpharmaceutical intervention (NPI), also known as physical distancing, which refers to the measures taken to prevent the spread of COVID-19 disease by maintaining a specific physical distance between individuals and reducing the number of time persons come in close contact with one another.

5.2. Treatment

This involves the use of an agent, procedure, or regimen, such as a drug to cure or mitigate the COVID-19 disease. This treatment employed here is a nonspecific regimen.

Here, we introduce these time-dependent interventions where is the time-dependent preventive effort which we refer to as social distancing, while is the time-dependent treatment effort which is an unspecific regimen to slow the spread of COVID-19. These two control functions are expected to be bounded and Lebesgue integrable on the closed interval , where is the duration of time for which we apply our control measures. The COVID-19 models (1), (2), (3), (4), (5), and (7) becomewith initial conditions

The model equation can be written in the vectorial form

, for where the model parameters are nonnegative. The goal here is to determine the optimal control strategy which minimizes the number of symptomatic infectious humans and the cost of the control measures and to predict the impact of these control measures as a means to advice public health officials on the best control and/or elimination policies. The control functions and are defined on the closed interval , where and The objective functional is defined mathematically which is given by

where , , and denote the weight constants (it is a parameter which describes the comparative importance of the two terms in the functional) for the relative costs of and the control interventions. Hence, there is need to obtain an optimal control pair such that

where

5.3. The Optimal Control Problem Statement

Here, we present the optimal control problem for the COVID-19 disease which is stated as follows:subject to the dynamicswith initial conditions

The terminal or final time is fixed and subjected to the control constraints.

There is a need to solve the following optimal control problem such that we find a control , which minimizes the objective functional, that is

If there exist such , it is referred to as optimal control. The optimal control along side with the corresponding solution gives the optimal control pair

The first question that we must address is to confirm whether an optimal control pair exists. Thus, this question of existence is settled by the following Lemma.

Lemma 8. (Filippov-Cesari existence theorem). For all , define the setSuppose that(1) is convex for every (2) is compact(3)There exists a constant such that for all and all admissible pairs Then, is an optimal pair where

Lemma 9. Suppose that there exists a solution of the optimal control problem (68)–(75).

Proof. We rewrite systems (68)–(73) as followswhere . We also define a set according to the Filippov-Cesari theorem [27], [26].
Here, we first show that is a convex for
Suppose that , we show that is a convex for each by proving that the line connecting and is completely in Thus, we need show that . For and that the control function such that for Let , and .
Thus, we haveSuppose that , it is observed that and if we let , the first part of the convex combination belongs to The second component will be checked as follows:Suppose that . It is observed that belongs to . Hence, we say finally that the convex combination of is in and is clearly compact. Secondly, we show that the solution of (69)–(73) is bounded.
For all initial conditions in , the trajectory obtained from the initial condition remains in a bounded domain included in
The total population is given by with the derivative .
It is observed thatWe have that where is the solution of equationThus, If then If any solution exists, we can find it with the Pontryagin’s maximum principle (PMP). This is done first, by incorporating a time-varying Lagrange multiplier , whose elements are referred to as the adjoint or costate variables. We hereby analyse the model equations (68), (69), (70), (71), (72), (73), and (75) with control variables . We seek an optimal control by applying the Pontryagin’s maximum principle (PMP) ([28]) and reformulating (68), (69), (70), (71), (72), (73), and (75) into a problem of minimizing pointwise a Hamiltonian , with respect to the control variables . Hence, the Hamiltonian is defined for all bywhere , , , , represent the adjoint variables or costate variables. Next, we obtain the following using the Pontryagin maximum principle (PMP) and the existence result of the Filippov-Cesari theorem.

Lemma 10. Suppose an optimal control and and the corresponding trajectories of state systems (68), (69), (70), (71), (72), and (73) which minimize over . Then, there exist adjoint variables , , , , and which satisfyand with transversality conditions and controls and satisfy the optimality conditions

Proof. We obtained the nonlinear system of differential equations for the adjoint variables by differentiating the Hamiltonian function , which is evaluated at the optimal control such that the adjoint/costate system is given byand with transversality conditions and the controls.
From the optimality conditionwe obtainAnd we obtainThen we haveTo obtain the optimal control and the corresponding prevalence , we therefore solve the following systemwith initial conditionswhereand with transversality conditions where and are given by (91) and (92). We cannot solve systems (95), (96), (97), (98), and (99) manually but numerical methods must be employed.

We hereby present the numerical solution of the optimality system and its corresponding optimal control pairs, the parameter values, and the various scenarios and interpretation in the next section.

6. Numerical Results

In this section, we present our numerical findings based on the optimal transmission parameter control measures for the COVID-19 dynamic model. We obtained the optimal control by solving the optimality system which involves the nonlinear differential equations (state equation) and the adjoint or costate equations. We applied the iterative scheme to solve the optimality system, and using the fourth-order Runge-Kutta method, we start to solve the state equations while we guess for the control measures over the simulated time. We solved the adjoint equations by a backward fourth-order Runge-Kutta scheme using the present iterative solutions of the state equation because of the transversality conditions. We employed a convex combination of the previous controls to update the controls and the value from the characterizations of the optimality conditions. We repeat the process, and iterations are truncated if the values of the unknowns at the preceding iteration are near to the present iterations ([29]). A COVID-19 model with preventive and treatment as control measures was investigated to predict the effects of control practices and the transmission of COVID-19. Building various scenarios of the two control measures either one control measure at a time or two control measures at a time, we examine and compare numerical results from simulations using the following scenarios.(1)Strategy A: employing social distancing without treatment (2)Strategy B: treatment of the symptomatic individuals without using social distancing (3)Strategy C: employing all the two control measures

The parameter values are given in Tables 1 and 3, and we obtained the optimal control strategy by the forward-backward sweep Runge-Kutta method of order 4. The weight constants , , and . We ran the simulations for days, and the results from the simulations with parameters in Table 3 are presented in Figures 10, 11, and 12.

In strategy A, we use only the control measure to optimize the objective functional , while the control measure on treatment is set to zero. The results in Figure 10(a) reveal that the control strategy provide a decrease in the number of symptomatic human as against the increase in the uncontrolled case. The control profile is presented in Figure 10(b), where the optimal social distancing control remains constant till the final time.

Applying this strategy B, we employed the treatment control measure only while we set the control on social distancing to zero. We observed that Figure 11(a) showed that effective treatment only has significant impact in decreasing COVID-19 incidence in the population, while the control profile in Figure 11(b) showed also that the control remains constant till the final time.

In strategy C, we applied the two controls , to optimise the objective functional . For this strategy, in Figure 12(a), it was observed that the control strategy had significant impact in reducing the number of symptomatic humans as against the increased number of cases under the uncontrolled case. The control profile in Figure 12(b) showed also that the control remains constant till the final time.

7. Conclusion

Hence, we presented a COVID-19 model with partially recovered carriers. We investigated the stability of the model which revealed that the disease-free equilibrium is locally and globally asymptotically stable. Scatter plots and the tornado plots of the parameters in (reproduction number) revealed that the contact rate between susceptible and infectious individuals () has a major impact on the transmission of COVID-19 disease, followed by the transmission rate from the exposed to infectious class, . This supports the fact that since the transmission rate is very high, the disease can spread fast in the population. We also observed from our results that the influence of the partially recovered contact rate is lower compared to the contact rate between susceptible and infectious individuals. The scenarios built on the optimal control strategies showed that the use of social distancing and treatment are the best option to control the disease in that they reduce the impact of the epidemic in the community and slows down the epidemic curve. Other interesting modeling works can be found in [30, 31]. The developed models predicted the reduction and control of COVID-19 through incorporating multiple control interventions.

Data Availability

The data used to support the findings of this study are included within the article. The authors used a set of parameter values whose sources are from the literature and others are estimated depending on the epidemiology of COVID-19 in Nigeria as shown in Table 1.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The second author would like to thank the Faculty of Science in the Federal University Oye Ekiti, Ekiti State, Nigeria. All authors would like to thank the Faculty of Basic Medical and Applied Sciences in Lead City University, Ibadan, Oyo state, Nigeria.