1 Introduction

It is no longer news that the current global threat to human existence is the novel coronavirus disease (COVID-19), which is caused by a new strain of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The disease was first reported to the World Health Organization (WHO) in late December, 2019 in Wuhan, Hubei Province, China [1, 2]. COVID-19 has been declared as a pandemic (global epidemic) affecting over 200 countries and territories around the world. The disease is responsible for the current global socioeconomic disruptions affecting lives and livelihoods. As of July 15, 2020, there were 13,690,115 confirmed cases of COVID-19, including 5,071,482 active cases and 587,235 deaths worldwide [3].

Dry cough, fever and fatigue (tiredness) are the most common symptoms of the novel coronavirus disease. However, other symptoms include shivering (chills), pains, sore throat, difficulty in breathing, headache, skin rashes, runny nose, taste loss, diarrhea and fingers or toes dislocation [2, 4]. COVID-19 primarily spreads when an infected person expels respiratory droplets via coughing or sneezing. Some reports have shown that infected person with no symptoms can also transmit the virus [5,6,7]. There is no specific cure yet for COVID-19 due to its novelty. However, efforts are in place to develop vaccines and antivirals to cure the disease. Existing management of cases is targeted at relieving symptoms of the disease while the immune system of infected person fights the virus. It should be mentioned that older people and persons with underlying health conditions such as lung disease, asthma, cancer, diabetes, liver disease and other immunocompromised conditions are at higher risk of having complications from COVID-19 [4]. The novel coronavirus disease can be prevented by adhering to measures, including regular hand washing or use of alcohol-based sanitizer, maintaining social distancing, avoiding crowded places, use of face masks, and use of hand gloves with personal protective clothing by healthcare workers while giving care to COVID-19 patients [2, 4].

On the 27th of February 2020, one person with immediate travel history was diagnosed of COVID-19 in Lagos State. This marked the first confirmed case of the novel coronavirus to be reported in Nigeria [8]. As of July 15, 2020, the cumulative number of confirmed cases has risen to 34,259 as reported by the Nigeria Centre for Disease Control (NCDC). Of the confirmed cases, 13,999 individuals have recovered while 760 deaths were recorded, bringing the total number of currently infected cases being managed at isolation centres or health facilities (active cases) to 19,500. It is important to note that all 36 states including the Federal Capital Territory of the country have reported cases of COVID-19, with Kogi and Cross River States having the least, 5 confirmed cases each, and Lagos State recording the highest number with 12,941 confirmed cases as at July 15, 2020 [8].

Mathematical modelling continues to play significant role in demystifying complex systems whose behaviors or properties can be observed. Most recently, a number of models fitted to real data of COVID-19 transmission have been developed to identify factors driving the disease at different population levels [9,10,11,12,13,14,15,16,17,18]. For examples, in [9], the authors analyzed COVID-19 spread dynamics with environmental influence using the reported cases in Ghana. Chen et al. [10] developed a Reservior-People transmission network model to assess the transmissibility of SARS-CoV-2 by estimating the number of secondary infections that result from introducing a single infected reservoir or person into an otherwise susceptible population. A mathematical model for assessing the impacts of non-pharmaceutical interventions on the transmission dynamics of the novel coronavirus pandemic in Nigeria was presented in [11] by using the available cumulative COVID-19 mortality data. Jia et al. [12] presented a dynamical model of COVID-19 by taking into account the impact of policy intervention and meteorological factors based on officially published data from China. Authors in [13] proposed a conceptual mathematical framework to explore the role of government and individual behavioral reaction on the outbreak of COVID-19 using South Africa reported cases. In addition, Nabi [14] proposed and calibrated a deterministic compartmental model for describing the transmission dynamics of the novel coronavirus disease based on the publicly available epidemiological data for Brazil, Russia, India and Bangladesh. Ullah and Khan [17] formulated a mathematical model for COVID-19 with non-pharmaceutical interventions by considering real infected cases of Pakistan to estimate parameters, and also analyzed the corresponding optimal control model using optimal control theory. In another development, Srivastav et al. [16] applied optimal control theory to analyze a mathematical model for COVID-19 model using real data from Italy and Spain.

This study is motivated by the ongoing community spread of the novel coronavirus disease in the most populous African nation – Nigeria with a view to flattening the infection curve. To achieve this aim, a new nonlinear deterministic epidemic model for COVID-19 incorporating transmission routes from three infectious classes, including symptomatic infectious humans, asymptomatic infectious humans and hospitalized individuals, is parameterized and analyzed based on the cumulative number of reported active cases. The behavior of the model over a long period of time is qualitatively and quantitatively investigated. Sensitivities to changes in the parameters of the model with respect to the epidemiological threshold are evaluated to provide insights into formulation of policies required to curb the spread of the disease. Moreover, optimal control theory is applied to the epidemiology of COVID-19 transmission in order to assess the optimal levels of time-dependent prevention and management measures that will significantly minimize the number of infectious humans in the population. The cost-effective intervention strategy capable of flattening the virus infection curve over a definite time horizon is suggested.

2 Model formulation

In this section, an epidemic model which focuses only on the dynamics of COVID-19 outbreak in the absence of demographic effects (natural births and deaths) is considered (see, e.g. [19, 20]). Thus, the host population is stratified into seven mutually exclusive epidemiological classes: susceptible class S(t) (those who are prone to contract COVID-19 infection), exposed class E(t) (those who have been infected, but are not infectious), symptomatic infectious class I(t) (those who exhibit COVID-19 symptoms and are capable of spreading the disease), asymptomatic infectious class A(t) (those who are not symptomatic, but are infectious), hospitalized class H(t) (infectious individuals who are admitted to healthcare facility due to virus infection (active cases)), recovery class R(t) (those who have been recovered from the disease), and death class D(t) (those who have died as a result of COVID-19 infection).

It is pertinent to mention that over 800 healthcare workers have tested positive for COVID-19 in Nigeria due to inadequate protective equipment [21]. This necessitates the transmission route from hospitalized class as appeared in [20]. It is also important to state that the remains of persons (i.e., symptomatic infectious and hospitalized individuals) suspected or confirmed to have died due to COVID-19 are carefully handled in accordance with the recommended standard procedure [8]. Further, incidence of death due to the presence of the virus in asymptomatic individuals is unknown, thus the disease-induced death rate for asymptomatic class is neglected. Figure 1 shows a flowchart depicting the transmission dynamics of COVID-19 in the total host population size N, which is given by \(N=S(t)+E(t)+I(t)+A(t)+H(t)+R(t)+D(t)\). Consequently, the following system of nonlinear differential equations is obtained.

Fig. 1
figure 1

A flowchart of COVID-19 model with \(\lambda =\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)/N-D\)

Full size image
$$\begin{aligned} \dfrac{dS}{dt}= & {} -\dfrac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)}{N-D}\nonumber \\ \dfrac{dE}{dt}= & {} \dfrac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)}{N-D}-\alpha E\nonumber \\ \dfrac{dI}{dt}= & {} l_{1}\alpha E-(h_{1}+r_{1}+\delta _{1})I\nonumber \\ \dfrac{dA}{dt}= & {} (1-l_{1})\alpha E-(h_{2}+r_{2})A\nonumber \\ \dfrac{dH}{dt}= & {} h_{1}I+h_{2}A-(\gamma +\delta _{2})H\nonumber \\ \dfrac{dR}{dt}= & {} r_{1}I+r_{2}A+\gamma H\nonumber \\ \dfrac{dD}{dt}= & {} \delta _{1}I+\delta _{2}H. \end{aligned}$$
(1)

Table 1 provides the detailed descriptions of the parameters of the COVID-19 model (1).

Table 1 The parameters of the COVID-19 model (1)
Full size table

3 Theoretical analysis of the model

It should be mentioned that the parameters of the formulated model (1) are non-negative since the model describes disease outbreak in human population. Therefore, it is suffices to state that solutions of the model (1) are non-negative.

3.1 Basic reproduction number

The potential spread of disease outbreak in a population can be determined by finding an epidemiological threshold known as the basic reproduction number. In the context of the novel coronavirus, the basic reproduction number, denoted by \({\mathfrak {R}}_{0}\), can be defined as the average number of secondary cases produced by one COVID-19 infected person during the course of infectiousness in a wholly susceptible population. To obtain \({\mathfrak {R}}_{0}\) for the COVID-19 model (1), the next generation operator technique is employed following the notations in [22]. Matrices F and V evaluated at the disease-free equilibrium (N, 0, 0, 0, 0, 0, 0) are given, respectively, by

$$\begin{aligned} F=\left( \begin{array}{c@{\quad }c@{\quad }c@{\quad }c} 0&{}\beta &{}\beta \varepsilon _{1}&{}\beta \varepsilon _{2}\\ 0&{}0&{}0&{}0\\ 0&{}0&{}0&{}0\\ 0&{}0&{}0&{}0 \end{array}\right) \end{aligned}$$

and

$$\begin{aligned} V=\left( \begin{array}{cccc} \alpha &{}0&{}0&{}0\\ -l_{1}\alpha &{}h_{1}+r_{1}+\delta _{1}&{}0&{}0\\ -(1-l_1)\alpha &{}0&{}h_{2}+r_{2}&{}0\\ 0&{}-h_{1}&{}-h_{2}&{}\gamma +\delta _{2} \end{array}\right) . \end{aligned}$$

The spectral radius of \(FV^{-1}\), given by \({\mathfrak {R}}_{0}=\rho (FV^{-1})\), is obtained as

$$\begin{aligned} {\mathfrak {R}}_{0}=\frac{\beta l_{1}}{k_{1}}+\frac{\beta \varepsilon _{1} (1-l_{1})}{k_{2}}+\frac{\beta \varepsilon _{2}( l_{1}h_{1}k_{2}+(1-l_{1})h_{2}k_{1})}{k_{1}k_{2}k_{3}}, \end{aligned}$$
(2)

where, \(k_{1}=h_{1}+r_{1}+\delta _{1}\), \(k_{2}=h_{2}+r_{2}\) and \(k_{3}=\gamma +\delta _2\).

Thus, the following local stability result is claimed in the sense of [22].

Theorem 1

The disease-free equilibrium (DFE) of the COVID-19 model (1) is locally asymptotically stable if the basic reproduction number, \({\mathfrak {R}}_{0}<1\), and unstable if \({\mathfrak {R}}_{0}>1\).

The implication of Theorem 1 from epidemiological viewpoint is that a small number of infected individuals present in a population will not cause an outbreak as long as \({\mathfrak {R}}_{0}<1\). It must be emphasized that this result depends on the initial sizes of infected individuals in the population. Hence, to rule out this dependency, a global stability result becomes necessary.

3.2 Global asymptotic stability of DFE

The global asymptotic stability of the COVID-19 free equilibrium, given by (N, 0, 0, 0, 0, 0, 0), is explored as follows.

Theorem 2

The disease-free equilibrium (DFE) of the COVID-19 model (1) is globally asymptotically stable if \({\mathfrak {R}}_{0}\le 1\).

Proof

The following continuously differentiable and positive definite linear Lyapunov function is carefully constructed for the COVID-19 model (1) (see, e.g., [23, 24]).

$$\begin{aligned} {\mathfrak {F}}=f_{1}E+f_{2}I+f_{3}A+f_{4}H, \end{aligned}$$
(3)

where, \(f_{1}=[l_{1}k_{2}k_{3}+\varepsilon _{1}(1-l_{1})k_{1}k_{3} +\varepsilon _{2}(l_{1}h_{1}k_{2}+(1-l_{1})h_{2}k_{1})]/k_{1}k_{2}k_{3}\), \(f_{2}=[k_{3}+\varepsilon _{2}h_{1}]/k_{1}k_{3}\), \(f_{3}=[\varepsilon _{1}k_{3}+\varepsilon _{2}h_{2}]/k_{2}k_{3}\) and \(f_{4}=\varepsilon _{2}/k_{3}\).

The time derivative of (3) along the solution path of the COVID-19 model (1) is given by

$$\begin{aligned} \dot{{\mathfrak {F}}}= & {} f_{1}{\dot{E}}+f_{2}{\dot{I}}+f_{3}{\dot{A}}+f_{4}{\dot{H}},\nonumber \\= & {} f_{1}\left( \frac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)}{N-D}-\alpha E\right) +f_{2}(l_{1}\alpha E-k_{1}I)+f_{3}((1-l_{1})\alpha E-k_{2}A)\nonumber \\&+f_{4}(h_{1}I+h_{2}A-k_{3}H),\nonumber \\= & {} f_{1}\left( \frac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)}{N-D}\right) -I-\varepsilon _{1}A -\varepsilon _{2}H,\nonumber \\\le & {} f_{1}\beta (I+\varepsilon _{1}A+\varepsilon _{2}H)-(I+\varepsilon _{1}A +\varepsilon _{2}H);\;\;\;(\text {since}\;S\le N-D),\nonumber \\= & {} (I+\varepsilon _{1}A+\varepsilon _{2}H)(f_{1}\beta -1),\nonumber \\= & {} (I+\varepsilon _{1}A+\varepsilon _{2}H)\left( \frac{\beta l_{1}}{k_{1}}+\frac{\beta \varepsilon _{1} (1-l_{1})}{k_{2}}+\frac{\beta \varepsilon _{2}( l_{1}h_{1}k_{2}+(1-l_{1})h_{2}k_{1})}{k_{1}k_{2}k_{3}}-1\right) ,\nonumber \\= & {} (I+\varepsilon _{1}A+\varepsilon _{2}H)({\mathfrak {R}}_{0}-1). \end{aligned}$$
(4)

It follows that \(\dot{{\mathfrak {F}}}\le 0\) if \({\mathfrak {R}}_{0}<1\) with equality, \(\dot{{\mathfrak {F}}}=0\), if and only if \(E=0\), \(I=0\), \(A=0\) and \(H=0\). Therefore, it can be concluded from the LaSalle’s Invariance Principle [25] that the DFE is globally asymptotically stable whenever \({\mathfrak {R}}_{0}\le 1\). \(\square \)

The implication of Theorem 2 from the epidemiological viewpoint is that COVID-19 can be effectively controlled in the population if \({\mathfrak {R}}_{0}<1\) regardless of the initial sizes of the infected individuals. This theoretical result will be graphically illustrated later.

4 Parameter estimation and model fitting

This section is devoted to the estimation of parameters of the COVID-19 model (1) based on the reported data of cases in Nigeria [8]. In particular, the formulated model is parameterized to assess the burden of COVID-19 in terms of the number of individuals currently infected (active cases) in Nigeria. In this study, the values for parameters such as effective transmission coefficient (\(\beta \)), modification parameter (\(\varepsilon _{2}\)), and hospitalization rates (\(h_{1}\) and \(h_{2}\)) are estimated, while the values of other parameters are chosen from the well-established literature as appeared in Table 2. The estimation of parameters is carried out using the least squares method implemented in Excel Solver [26] with a view to minimizing summation of squared errors given by \(\sum (Y(t,p)-X_{real})^2\) subject to the COVID-19 model (1), where \(X_{real}\) is the real reported data, and Y(t, p) represents the solution of the model corresponding to the number of active cases over time t with the set of estimated parameters, denoted by p.

Moreover, initial conditions for the exposed humans, symptomatic and asymptomatic infectious individuals are estimated to account for the possibility that the actual number of cases is most likely much higher than reported in Nigeria as in other parts of the world, due to limited testing and the possible presence of COVID-19 infected individuals with no symptoms [27]. In order to accommodate the foregoing possibility, the estimations of parameters and initial conditions for the model fitting are allowed to commence from March 10, 2020, exactly a day after the second COVID-19 confirmed case was reported by NCDC [8], so that \(H(0)=2\). Keeping in mind that the total population of Nigeria is estimated at \(N=200\) million people [28], and that no report of recovery and deaths due to the disease was made as at March 10, 2020 (i.e., \(R(0)=D(0)=0\)), then it is easy to determine the initial susceptible population as \(S(0)=N-(E(0)+A(0)+I(0)+H(0)+R(0)+D(0))\).

Fig. 2
figure 2

Fitting of the COVID-19 model (1) with the available cumulative number of active cases from March 10 to July 15, 2020

Full size image
Table 2 Values of parameters and initial conditions
Full size table

The model is fitted to the real data along with the following initial conditions: \(E(0)=880\), \(A(0)=330\), \(I(0)=190\), so that \(S(0)=199,998,600\). The model fitting with reported active cases is depicted in Fig. 2. Consequently, the value of the basic reproduction number using the parameter values in Table 2 is \({\mathfrak {R}}_{0}=1.412\), which still exceeds the novel coronavirus pandemic threshold of one.

4.1 Sensitivity of the model to parameters

It is important to explore how sensitive the COVID-19 model (1) is to changes in each of its parameters in order to suggest intervention strategies that will help in bringing down the infection trajectory. In other words, carrying out sensitivity analysis will help in providing insights into what should be done or avoided to prevent the outbreak of the novel coronavirus [32]. Thus, as defined in [33], a relative change in the basic reproduction number, \({\mathfrak {R}}_{0}\), when each parameter changes is measured using the following

$$\begin{aligned} \Upsilon ^{{\mathfrak {R}}_{0}}_p=\frac{\partial {\mathfrak {R}}_{0}}{\partial p}\times \frac{p}{{\mathfrak {R}}_{0}}, \end{aligned}$$

where, \(\Upsilon ^{{\mathfrak {R}}_{0}}_p\) is the sensitivity index of a differentiable \({\mathfrak {R}}_{0}\) with respect to any parameter, p. Consequently, the sensitivity indexes for the model (1) are graphically shown in Fig. 3. It can be observed that of all the positive indices, the effective transmission coefficient, \(\beta \), is the highest, and therefore the most sensitive parameter. This means that an increase (or a decrease) of the value of \(\beta \) will increase (or decrease) \({\mathfrak {R}}_{0}\) by 100%. However, of all the negative indices shown in Fig. 3, the recovery rate for the hospitalized class (active cases), denoted by \(\gamma \), is the most sensitive parameter. This means that an increase (or a decrease) of the value of \(\gamma \) will decrease (or increase) \({\mathfrak {R}}_{0}\) by 42.1%.

In addition, contour plots of \({\mathfrak {R}}_{0}\) as a function of other parameters are shown in Fig. 4 to demonstrate how changes in these parameters affect the basic reproduction number, \({\mathfrak {R}}_{0}\). Figure 4a shows a decrease in \({\mathfrak {R}}_{0}\) with increasing recovery rates for both symptomatic and asymptomatic individuals. This spontaneous recovery from infection is due to the body’s innate immune response which is capable of inhibiting virus replication, promoting virus clearance, and triggering a prolonged adaptive immune response against reinfection [34]. Moreover, it can be confirmed, as shown in Fig. 4b, that the basic reproduction number can be brought lower than the threshold of one if efforts are geared towards lowering the transmission rate while simultaneously improving on the management of active cases.

Fig. 3
figure 3

COVID-19 model sensitivities to its associated parameters

Full size image
Fig. 4
figure 4

2-D contour plots of the basic reproduction number \({\mathfrak {R}}_{0}\)

Full size image
Fig. 5
figure 5

Projections with varying effects of parameters: a Effective transmission coefficient at values of \(\beta =[0.42\; ({\mathfrak {R}}_{0}=1.521>1),\; 0.38974 \;({\mathfrak {R}}_{0}=1.412>1),\; 0.36 \;({\mathfrak {R}}_{0}=1.304>1) \;\text {and}\; 0.19487\; ({\mathfrak {R}}_{0}=0.706<1)]\). b Recovery rate for hospitalized individuals at values of \(\gamma =[0.06667\; ({\mathfrak {R}}_{0}=1.412>1),\; 0.1 \;({\mathfrak {R}}_{0}=1.2>1),\; \text {and}\; 0.2\; ({\mathfrak {R}}_{0}=0.959<1)]\)

Full size image
Fig. 6
figure 6

Relationship between COVID-19 epidemiological threshold and parameters

Full size image

In Fig. 5a, the effects of effective transmission coefficient at different values are illustrated. It is projected that the cumulative number of active cases in Nigeria, given the current infection trajectory with \(\beta =0.38974\), may get up to \(5\times 10^{6}\) (equivalently, 2.5% of the total estimated population) unless drastic measures are taken to bring the basic reproduction number below one. As shown in Fig. 5b, the cumulative number of active cases reduces with increasing value of recovery rate for hospitalized individuals from \(\gamma =1/15\;\text {day}^{-1}\) to \(1/5\;\text {day}^{-1}\). In particular, if it takes 5 days or less to effectively manage infected individuals, the risk of infection transmission from the hospitalized class will reduce significantly. This indicates that the longer it takes to effectively manage COVID-19 patient, the heavier the burden of the disease on the healthcare facility and the population at large.

More explicitly, Fig. 6 reveals how the potential spread of the novel coronavirus in the population is affected by the model parameters in order to influence effective policies that can guide against the burden of the disease transmission. The blue rectangular box in Fig. 6a indicates a safe region where the effective transmission coefficient, \(\beta \), is below a value of 0.27613, so that the key epidemiological threshold, \({\mathfrak {R}}_{0}\), is less than one. Above this region, the basic reproduction number which determines the spread potential of the disease increases as the rate of transmission increases through effective contact between susceptible and infectious individuals in the population. Hence, to stay in the safe region, measures such as strict social distancing and maintenance of good personal hygiene are strongly advised. This result is consistent with previous findings in the literature [9, 11, 13, 15, 17]. Furthermore, the potential spread of the diseases decreases with increase in the recovery rate as shown in Fig. 6b. The safe region of \({\mathfrak {R}}_{0}<1\) in this case is attained at recovery rate above the value of 0.1727. Hence, given the high infectiousness of the novel coronavirus disease, treatment regime which ensures recovery within a short period of 5 days is strongly advised. In Fig. 6c, the effects of the modification parameters, \(\varepsilon _{1}\) and \(\varepsilon _{2}\), respectively for virus transmission from asymptomatic (black solid line) and hospitalized (blue solid line) individuals on the basic reproduction number are shown. The red rectangular box is the safe region where \(\varepsilon _{2}\) is below a value of 0.1. Above this value, \({\mathfrak {R}}_{0}\) will increase up to 3.6828 which may lead to a more serious pandemic situation. Hence, strict usage of personal protective equipment when attending to the hospitalized COVID-19 individuals is strongly advised to avoid further spread of the disease.

To validate the global stability result (Theorem 2), Fig. 7 illustrates the convergence to the COVID-19-free equilibrium irrespective of the initial sizes of the infectious individuals in the population.

Fig. 7
figure 7

Convergence of solution trajectories for a symptomatic humans, b asymptomatic humans and c hospitalized humans at different initial sizes in line with Theorem 2. Parameter values used are as given in Table 2 except for \(\beta =0.19487\) and \(\varepsilon =0.12\), so that \({\mathfrak {R}}_{0}=0.522<1\)

Full size image

5 Optimal control model and analysis

In the light of the sensitivity analysis result, optimal control theory is applied to the COVID-19 model (1), as applied to other infectious diseases models (e.g., [33, 35,36,37,38,39]) in order to mitigate the spread of COVID-19 in the population optimally. This is done by introducing two time-dependent control variables \(u_{1}(t)\) and \(u_{2}(t)\), which are described in details as follows.

  1. (i)

    \(u_{1}(t)\) denotes the preventive strategy targeted at inhibiting the virus transmission from symptomatic, asymptomatic and hospitalized humans. This can be achieved through public heath advocacy for social distancing, good personal hygiene, wearing face masks in public places, and provision of protective gear for healthcare workers. Noting that \(u_{1}(t)=1\) implies the strategy is effective for protection against infection, while \(u_{1}(t)=0\) implies strategy failure.

  2. (ii)

    \(u_{2}(t)\) represents control variable to step up the management of hospitalized individuals with a view to ensuring prompt recovery and preventing deaths due to complications. This can be achieved through prompt provision of supplemental oxygen or mechanical ventilation for hospitalized individuals with severe COVID-19. If \(u_{2}(t)=1\), then the control is effective in managing the disease, while \(u_{2}(t)=0\) means the absence of control.

Consequently, the optimal control model with the two aforementioned time-dependent variables is given by the following differential equations

$$\begin{aligned} \begin{aligned} \dfrac{dS}{dt}&=-(1-u_{1}(t))\dfrac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)}{N-D}\\ \dfrac{dE}{dt}&=(1-u_{1}(t))\dfrac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)}{N-D}-\alpha E \\ \dfrac{dI}{dt}&=l_{1}\alpha E-(h_{1}+r_{1}+\delta _{1})I\\ \dfrac{dA}{dt}&=(1-l_{1})\alpha E-(h_{2}+r_{2})A\\ \dfrac{dH}{dt}&=h_{1}I+h_{2}A-(\gamma +u_{2}(t)+\delta _{2})H\\ \dfrac{dR}{dt}&=r_{1}I+r_{2}A+(\gamma +u_{2}(t))H\\ \dfrac{dD}{dt}&=\delta _{1}I+\delta _{2}H. \end{aligned} \end{aligned}$$
(5)

The aim of introducing the two control variables is to seek the optimal solution required to minimize the numbers of symptomatic, asymptomatic and hospitalized individuals responsible for spreading the novel coronavirus in the population at minimum cost. Hence, the objective functional for this control problem is given by

$$\begin{aligned} \begin{array}{ll} {\mathfrak {J}}(u_1,u_2)=\displaystyle \min _{0\le u_{1},u_{2}\le 1}\displaystyle \int ^{T_f}_{0} \left( w_{1}I+w_{2}A+w_{3}H+\frac{1}{2}[w_4 u_{1}^2(t)+w_5 u_{2}^2(t)]\right) dt, \end{array} \end{aligned}$$
(6)

where, constants \(w_i,\;i=1,2,...,5\) are positive weights required to balance the corresponding terms in the objective functional. Following other literature on COVID-19 control problems [9, 16, 17, 40, 41], quadratic costs on the controls are chosen, where \(1/2w_4 u_{1}^2(t)\) is the total cost of implementing the preventive measure, and \(1/2w_5 u_{2}^2(t)\) is the total cost of managing active cases over the time interval \([0,T_f]\).

Precisely, the optimal control double \(u^*=(u^*_1, u^*_2)\) is sought such that

$$\begin{aligned} {\mathfrak {J}}(u_{1}^*,u_{2}^*)=\min \left\{ {\mathfrak {J}}(u_1,u_2): u_1,u_2\in {\mathcal {U}}\right\} , \end{aligned}$$
(7)

where, \({\mathcal {U}}\) is the non-empty control set defined by

\({\mathcal {U}}=\left\{ (u_1,u_2): (u_1(t),u_2(t)) \;\text {are measurable with}\; 0\le u_{1},u_2\le 1\;\text {for}\;t\in [0,T_f]\right\} \). Thus, to determine the necessary conditions that the optimal control double (\(u_{1}^*,u_{2}^*\)) must satisfy, Pontryagin’s maximum principle [42], which transforms the control problem (7) subject to the model (5) into a problem of minimizing pointwise a Hamiltonian \({\mathcal {H}}\), is used. This Hamiltonian is given by

$$\begin{aligned} {\mathcal {H}}= & {} w_{1}I+w_{2}A+w_{3}H+\frac{1}{2}[w_4 u_{1}^2(t)+w_5 u_{2}^2(t)]\nonumber \\&+\lambda _{1}\left[ -(1-u_{1}(t))\dfrac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)}{N-D}\right] \nonumber \\&+\lambda _{2}\left[ (1-u_{1}(t))\dfrac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)}{N-D}-\alpha E\right] \nonumber \\&+\lambda _{3}[l_{1}\alpha E-(h_{1}+r_{1}+\delta _{1})I]+\lambda _{4}[(1-l_{1})\alpha E-(h_{2}+r_{2})A]\nonumber \\&+\lambda _{5}[h_{1}I+h_{2}A-(\gamma +u_{2}(t)+\delta _{2})H]\nonumber \\&+\lambda _{6}[r_{1}I+r_{2}A+(\gamma +u_{2}(t))H]+\lambda _{7}[\delta _{1}I+\delta _{2}H], \end{aligned}$$
(8)

where, \(\lambda _{i}\), \(i=1,2,3,...,7\), represent the adjoint variables associated with the state variables of the model (5). The standard existence result for minimizing control problem as appeared in [43, 44] is adapted as follows.

Theorem 3

Given that \((u_{1}^{*},u_{2}^{*})\) minimizes the objective functional (6) subject to the corresponding state system (5), then the adjoint variables \(\lambda _{i}\), \(i=1,2,3,...,7\), satisfy the following system

$$\begin{aligned} \begin{aligned} \frac{d\lambda _{1}}{dt}&=\dfrac{\beta (1-u_{1}(t))(I+\varepsilon _{1}A +\varepsilon _{2}H)[\lambda _{2}-\lambda _{1}]}{N-D}\\ \frac{d\lambda _{2}}{dt}&=\alpha (\lambda _2-l_{1}\lambda _3-(1-l_1)\lambda _4)\\ \frac{d\lambda _{3}}{dt}&=\dfrac{\beta S(1-u_{1}(t))[\lambda _{2}-\lambda _{1}]}{N-D}+\delta _{1}(\lambda _3 -\lambda _7)+h_{1}(\lambda _3-\lambda _5)+r_{1}(\lambda _3-\lambda _6)-w_{1}\\ \frac{d\lambda _{4}}{dt}&=\dfrac{\varepsilon _{1}\beta S(1-u_{1}(t))[\lambda _{2}-\lambda _{1}]}{N-D}+h_{2}(\lambda _4-\lambda _5) +r_{2}(\lambda _4-\lambda _6)-w_2\\ \frac{d\lambda _{5}}{dt}&=\dfrac{\varepsilon _{2}\beta S(1-u_{1}(t))[\lambda _{2}-\lambda _{1}]}{N-D}+\delta _{2}(\lambda _5-\lambda _7) +(\gamma +u_2(t))(\lambda _5-\lambda _6)-w_3\\ \frac{d\lambda _{6}}{dt}&= 0\\ \frac{d\lambda _{7}}{dt}&=\dfrac{\beta (1-u_{1}(t))S(I+\varepsilon _{1}A+\varepsilon _{2}H)[\lambda _{1} -\lambda _{2}]}{(N-D)^{2}} \end{aligned} \end{aligned}$$
(9)

with the terminal (transversality) conditions

$$\begin{aligned} \begin{array}{l} \lambda _{i} (T_f)=0, \quad i=1,2,3,...,7. \end{array} \end{aligned}$$
(10)

Further, the optimal control double \((u^{*}_{1},u^{*}_{2})\) is given as follows

$$\begin{aligned} u^{*}_1= & {} \max \left\{ 0,\min \left\{ 1,\dfrac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)[\lambda _{2}-\lambda _{1}]}{w_4(N-D)}\right\} \right\} ,\nonumber \\ u^{*}_2= & {} \max \left\{ 0,\min \left\{ 1,\dfrac{H[\lambda _{5}-\lambda _{6}]}{w_5}\right\} \right\} . \end{aligned}$$
(11)

Proof

The existence of the optimal control double (\(u_{1},u_{2}\)) is based on the a priori boundedness of the state variables of model (5), convexity and boundedness of the Langragian of the control problem. The system of ordinary differential equations (9) governing the adjoint variables is derived by differentiating the Hamiltonian \({\mathcal {H}}\). Further, the control characterizations in (11) are derived by solving, on the interior of the control set \({\mathcal {U}}\), the partial differentials of the Hamiltonian \({\mathcal {H}}\) with respect to each of the controls \(u_{1}\) and \(u_{2}\). Hence, by standard arguments involving control bounds, it follows that

$$\begin{aligned} u_1^{*}=\left\{ \begin{array}{l@{\quad }l} 0&{}\;\text {if}\;\;\tau _{1}^*\le 0,\\ \tau _{1}^*&{}\;\text {if}\;\;0<\tau _{1}^*<1,\\ 1&{}\;\text {if}\;\;\tau _{1}^*\ge 1, \end{array}\right. \end{aligned}$$

and

$$\begin{aligned} u_2^{*}=\left\{ \begin{array}{l@{\quad }l} 0&{}\;\text {if}\;\;\tau _{2}^*\le 0,\\ \tau _{2}^*&{}\;\text {if}\;\;0<\tau _{2}^*<1,\\ 1&{}\;\text {if}\;\;\tau _{2}^*\ge 1, \end{array}\right. \end{aligned}$$

where,

$$\begin{aligned} \tau _{1}^*=\dfrac{\beta S(I+\varepsilon _{1}A+\varepsilon _{2}H)[\lambda _{2}-\lambda _{1}]}{w_4(N-D)} \end{aligned}$$

and

$$\begin{aligned} \tau _{2}^*=\dfrac{H[\lambda _{5}-\lambda _{6}]}{w_5}. \end{aligned}$$

This ends the proof. \(\square \)

5.1 Simulations of the control probem

Forward and backward Runge-Kutta method of order four implemented in MATLAB is used to solve the resulting optimality system which consists of systems (5) and (9) with the characterizations (11) within the time interval of [0, 100] days. The weight constants used for balancing the terms in the objective functional (6) are chosen to ensure that no term dominates the other. Thus, equal weight constant for minimizing the infectious classes is chosen, so that \(w_1=w_2=w_3=1\). On the other hand, the weight constants for measuring efforts or costs required to implement the controls are comparatively different. This results in values for \(w_4=50\) and \(w_5=100\). Details of the numerical procedure for simulating the obtained optimality system are contained in [45].

Fig. 8
figure 8

Control profile (\(u_{1}(t)\)) and its effects on the COVID-19 dynamics

Full size image
Fig. 9
figure 9

Control profile (\(u_{2}(t)\)) and its effects on the COVID-19 dynamics

Full size image
Fig. 10
figure 10

Combined effects of optimal controls \(u_1(t)\), \(u_2(t)\) on the COVID-19 dynamics

Full size image

Figure 8 demonstrates how single preventive measure, \(u_{1}(t)\), affects the spread of the novel coronavirus in the population. As shown in Fig. 8a, to minimize the objective functional (6), the optimal control \(u_{1}(t)\) is maintained at maximum level of 100% for about 93 days before relaxing to the minimum in final time. As expected, the sizes of COVID-19 infectious individuals are reduced when control is in place. Further, in Fig. 9a, b, the effects of single management control, \(u_{2}(t)\), are shown in minimizing the number of virus infection in the population, but not as much as when only preventive strategy is applied. More importantly, Fig. 10 shows the significance of combining the two optimal controls in bringing down the total number of infectious humans to zero. It is observed that optimal solution is attained when preventive measure is strictly adhered to at maximum level of 100% for 45 days, while the management control of the hospitalized individuals is at maximum level of 100% for 32 days. It can be seen that combination of the two controls is significantly more effective in curbing the spread of the virus than when each control is singly applied.

5.2 Cost-effectiveness analysis

It is important to determine the most cost-effective optimal control measure among the single and combined implementation of the two given control measures in order to optimally mitigate the spread of COVID-19 at the lowest cost possible. Hence, cost-effectiveness analysis is explored using the incremental cost effectiveness ratio (ICER) [9, 24, 46,47,48]. To avoid dissipation of available limited resources, ICER is used to compare any two competing measures for controlling the spread of disease or related problems. The ICER formula is given by

$$\begin{aligned} \text {ICER}=\frac{\text {Change in total costs between control measures}}{\text {Change in total number of infection averted by control measures}}. \end{aligned}$$
(12)

The total cost for each of the single implementation and combined effort of the control measures is derivable from the objective functional 6, while the infection averted is obtained by calculating the difference between infectious individuals without and with control measures. Let \(C_1\), \(C_2\) and \(C_{12}\) represent single implementation of the preventive measure \(u_{1}(t)\), single implementation of the management control \(u_{2}(t)\) and combined effort of the two measures, respectively. In what follows, Table 3 summarizes the ICER for each and combination of the control variables \(u_{1}(t)\) and \(u_{2}(t)\) in an increasing order of the total infection averted.

Table 3 ICER in the order of COVID-19 cases averted by control measures
Full size table

Using (12), the ICER for \(C_1\), \(C_2\) and \(C_{12}\) shown in Table 3 are calculated, respectively, as follows.

$$\begin{aligned} \text {ICER} (C_{1})= & {} \frac{4.7351\times 10^{4}}{3.5062\times 10^{6}}=0.0135,\\ \text {ICER} (C_{2})= & {} \frac{(8.2886-4.7351)\times 10^{4}}{(3.7311-3.5062)\times 10^{6}}=0.1580,\\ \text {ICER} (C_{12})= & {} \frac{(5.0064-8.2886)\times 10^{4}}{(3.8848-3.7311)\times 10^{6}}=-0.2136. \end{aligned}$$

Comparing \(C_{1}\) and \(C_{2}\), it is seen that ICER(\(C_1\)) is less than ICER(\(C_2\)). This means that \(C_{2}\) is more costly and less effective than \(C_1\). In other words, \(C_{1}\) dominates \(C_{2}\). Thus, single implementation of management control is removed from the list. As a consequence, \(C_1\) and \(C_{12}\) are assessed in Table 4 using similar procedure.

Table 4 Comparison between \(C_1\) and \(C_{12}\)
Full size table

It is revealed in Table 4 that \(C_{1}\) is dominated by \(C_{12}\) since ICER(\(C_{1}\)) is greater than ICER(\(C_{12}\)). This implies that \(C_{12}\) is less costly and more effective than \(C_{1}\). Hence, single implementation of preventive measure is excluded from the list. This shows that combined effort of the two control measures is the most cost-effective intervention capable of diminishing the burden of the novel coronavirus optimally in the host population.

6 Conclusion

This study presented a mathematical analysis of transmission dynamics of the novel coronavirus (COVID-19) pandemic with a view to providing further insights into the disease transmission, and to explore possible prevention and control measures capable of curbing the rising tide of the disease spread in the population. A deterministic mathematical model was formulated by subdividing the host population into susceptible, exposed, symptomatic infectious, asymptomatic infectious, hospitalized, recovered and dead classes. The nonlinear model has been developed by considering transmission routes from symptomatic infectious humans, asymptomatic infectious humans and hospitalized individuals. The theoretical analysis carried out, based on Lyapunov stability, showed that the model has a globally asymptotically stable disease-free equilibrium if the basic reproduction number of the novel coronavirus transmission is less than one. This is an indication that COVID-19 outbreak can be effectively controlled in the population irrespective of the number of the infectious individuals initially introduced into the completely susceptible population.

In particular, the formulated model was fitted to the reported data of cumulative number of hospitalized individuals (active cases) due to COVID-19 in Nigeria, and estimates for parameters such as, effective transmission coefficients (\(\beta \)), modification parameter for the transmission of virus infection from active cases (\(\varepsilon _2\)), hospitalization rates for both symptomatic and asymptomatic individuals (\(h_1\) and \(h_2\)) were determined using the least squares method. Sensitivity analysis of the fitted model was performed to find the parameters that drive the spread of the virus infection mostly in the population. It was revealed that \(\beta \) is the most sensitive parameter, followed by \(\varepsilon _{2}\) and the recovery rate for active cases denoted by \(\gamma \). In addition, given the high infectiousness of the novel coronavirus disease, safe regions where virus infection dies out at certain threshold values of the model parameters were derived in order to facilitate and strengthen formulation of effective policies that can help to avoid a more serious pandemic situation. Further insights into the projection of the disease transmission indicated that the cumulative number of active cases in Nigeria, given the current infection trajectory, may get up to 2.5% of the total estimated population of 200 million, unless drastic measures are taken to bring the basic reproduction number below one. It was shown by simulations that the basic reproduction number can be brought to a value less than one in Nigeria, if the current effective transmission rate of the disease can be reduced by 50%.

Moreover, the model was extended to include time-dependent optimal control variables (preventive and management measures) to examine their impacts in minimizing the burden of COVID-19 in the population using Pontryagin’s maximum principle. The optimal preventive measure was shown to be better than management control in reducing the burden of the disease, but the combined effort of the two controls has significant effect in reducing the number of infectious individuals in the population. This result was further supported by carrying out the cost-effectiveness analysis, and it was therefore established that the combined implementation of the two control measures is the most cost-effective when compared with the single implementation of each control measure.

In view of the above, it can be concluded that flattening of the novel coronavirus infection curve largely depends on two key parameters, namely effective transmission coefficient and recovery rate. Efforts geared toward lowering the transmission rate and improving on the treatment regime for quick recovery of cases are required to ensure virus-free region and effectively overcome the spread of the novel coronavirus in the host population. To achieve this optimally, it is recommended that adherence to the combination of preventive and management control measures should be strict. These measures include continued public advocacy for social distancing, good personal hygiene, wearing face masks in public places, provision of protective gear for healthcare workers caring for people with COVID-19 infection. If strict adherence is practiced, transmission routes from symptomatic, asymptomatic and hospitalized individuals will be effectively hindered, and availability of supplemental oxygen or mechanical ventilation for hospitalized individuals with severe COVID-19. In addition, given the ongoing community spread of the virus in Nigeria and almost all the countries of world, the need to scale up the testing capacity for detection and prompt treatment of asymptomatic individuals cannot be overemphasized in order to nip the transmission of the novel coronavirus in the bud. It is also pertinent to advise that healthcare facilities should not only be increased, but well equipped to accommodate the increasing spate of cases and effectively manage patients who may develop complications due the novel coronavirus infection.