key: cord-0891661-9rj2a7uk authors: Bassey, Bassey Echeng; Atsu, Jeremiah U. title: Global stability analysis of the role of multi-therapies and non-pharmaceutical treatment protocols for COVID-19 pandemic date: 2020-12-19 journal: Chaos Solitons Fractals DOI: 10.1016/j.chaos.2020.110574 sha: 922b518e69d2ebdcd913d06d2ec8a2b12007182d doc_id: 891661 cord_uid: 9rj2a7uk In this paper, we sought and presented an 8-Dimensional deterministic mathematical COVID-19 dynamic model that accounted for the global stability analysis of the role of dual-bilinear treatment protocols of COVID-19 infection. The model, which is characterized by human-to-human transmission mode was investigated using dual non-pharmaceutical (face-masking and social distancing) and dual pharmaceutical (hydroxylchloroquine and azithromycin) as control functions following the interplay of susceptible population and varying infectious population. First, we investigated the model state-space and then established and computed the system reproduction number for both off-treatment [Formula: see text] and for onset-treatment [Formula: see text]. We considered the model for off-treatment and thereafter by incorporating the theory of LaSalle's invariant principle into the classical method of Lyapunov functions, we presented an approach for global stability analysis of COVID-19 dynamics. Numerical verification of system theoretical predictions was computed using in-built Runge-Kutta of order of precision 4 in a Mathcad surface. The set approach produces highly significant results in the main text. For example, while rapid population extinction was observed by the susceptible under off-treatment scenario in the first [Formula: see text] days, the application of non-pharmaceuticals at early stage of infection proved very effective strategy in curtailing the spread of the virus. Moreso, the implementation of dual pharmacotherapies in conjunction with non-pharmaceuticals yields tremendous rejuvenation of susceptible population ([Formula: see text]) with maximal reduction in the rates of isolation, super spreaders and hospitalization of the infectives. Thus, experimental results of investigation affirm the suitability of proposed model for the control and treatment of the deadly disease provided individuals adheres to treatment protocols. Following the upsurge of what could be considered the most dreaded transmittable infectious disease in the history of virological infections in human race, the infectious disease known as coronavirus disease 2019 (COVID-19) has taken an unprecedented dehumanization of mankind with over 213 nations of the world falling prey to the deadly disease. COVID-19 in its nature, is a negative-gram of ribonucleic acid (RNA) virus that have both human and animals as its prey. Zoonotic scientists have shown that coronavirus infection have been identified as causing the disease of the type -severe acute respiratory syndrome (SARS-CoV-2), middle eastern respiratory syndrome (MERSS-CoV) and the most recent COIVID-19. Among these, ASRS-CoV-2 have been found to be the causative agent of the later -COVID-19 [1] . The human coro-to accommodate large number of asymptotic patients who could be infectious but with less clinical manifestations [ 3 , 6 ] . Generally, COVID-19 exhibits similar clinical symptoms as SARS-CoV-2 and MERS-CoV, which appear as typical pneumonia marked by cough, fever, headache, dry throat and subsequent onset acute respiratory syndrome -coronavirus 2 (SARS-CoV 2) with life-threating respiratory failure [7] . Victims of COVID-19 appears to cut across all human race with most vulnerable and severe cases occurring among the elderly (adults of age ≥ 65 years). Like most other infectious disease with peculiarities of nonavailability of outright medical cure and vaccines, varying nonpharmaceutical preventive and control intervention measures (in the range of hand-sanitizers, regular washing of hands, face-masks, social distancing, quarantization, isolation, contact tracing, hospitalization and lockdown) have been explored as useful control measures. For severe cases, the aspect of isolation and hospitalization have resulted to some level of clinical trials following the recommendation of pharmaceutical therapies (in the range of Dexamethasone, chloroquine, lopinavir/ritonavir, hydroxylchloroquine, azithromycin and erythromycin). Perturbingly, efforts to demystify the disease transmission dynamics and to propagate clinical methodological treatment protocols have attracted the attention of mathematical modeling. For instance, in the wake of COVID-19 pandemic, innovative findings have been formulated starting with the extensive evaluation of the impact of mathematics and mathematical models in understanding and controlling the 2019 novel coronavirus pandemic [1] . Noting the lack of vaccine for the control of the virus, [6] proposed the use of optimal quarantine strategies for the control of COVID-19, accounting for the long-term cost effect of the strategy. The outcome was massive even without vaccination provided survival level of viral load is kept less than 1. The model [8] proposed and studied the use of X-ray images characterized by hybrid 2D curvelet transform chaotic salp swarm algorithm and deep learning technique (CASSA) as a major alternative testing kits to the earlier proclaimed real-time reverse transcription-polymerase chain reaction (RT-PCR) by the World Health Organization (WHO). The model was analyzed using EfficientNet-BO method in conjunction with 2D curvelet transformation. Results showed that the model proved to be faster and low computation cost when compared to TR-PCR. The study [9] had considered the dynamical behavior of COVID-19 with carrier effect to outbreak epidemic. The model was formulated as a 5-Dimensional mathematical differential equations and studied using suitable Lyapunov function for the system global stability conditions. Numerical results obtained showed that the awareness as a single factor was not sufficient in reducing Covid-19 epidemic. The study proposed in addition to awareness, the inclusion of incidence rate, prevention rate and carrier as indices for the reduction of the spread of the virus. Incorporating Least-squares method into Lyapunov function, [10] investigated the global stability and cost-effectiveness of COVID-19, taking into account environmental factors and adopting six non-pharmaceutical control measures. Results indicated that the strategy involving practicing proper coughing etiquette, maintain distancing, covering cough and sneezing with disposable tissues and washing of hands is the most cost-effective strategy. Other mathematical models on COVID-19 with varying mathematical concepts can be found in [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] . Amidst these novel literatures on COVID-19, it is noticed that a standard mathematical model that accounted for the clinical and explicit combination of both pharmaceutical and nonpharmaceutical control strategies have not been given the desired attention. Thus, the present study considering human-to-human transmission mode, seek to formulate and analyze a set of standard mathematical model using dual-bilinear control treatment functions arising from dual non-pharmaceutical controls -use of face masks and social distancing and dual pharmaceutical therapieshydroxylchloroquine (HCQ) and Azithromycin (AZT). Resourcefully, with the introductory aspect in Section 1 , the mathematical and epidemiological presentation of the study is partition into seven sections. In Section 2 , we present the material and methods constituted by the system problem statement and derivation of mathematical model equations. In Section 3 , we discuss the mathematical analysis involving the investigation of the system state-space. Under system stability analysis in Section 4 , we investigated the existence of equilibrium states, system reproduction number, local stability in terms of reproduction number, system endemic equilibrium and the system global stability conditions. Numerical investigation of system theoretical predictions are presented in Section 5 . Section 6 is devoted to the discussion and analysis of obtained results. Finally, we domicile our succinct conclusion and incisive remarks in Section 7 . Notably, the present study is anticipated to unveil novel findings towards the annihilation of the deadly disease. The material and methods of this study is characterized by the problem statement and derivation of model equations as well as analysis of model basic mathematical properties. Following the advent of the unparalleled coronavirus pandemic, a number of notable mathematical models have been formulated in a move to mathematically present insight to the virus historical background, infection transmission dynamics and possible intervention and control measures with most models skewed to geographical locations. Taking advantage of the limitations of two seeming compactible models [ 12 , 16 ] in relation to the present investigation, we seek to formulate a holistic COVID-19 differential model, considered adaptable to any geographical location. For instance, the model by [16] had studied a mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan. Moreso, the model [12] had investigated the contribution of mathematical modeling of COVID-19 in the Niger Republic using 8-Dimensional differential equations. On critical review of these two models, we observed the following limitations: i Formulated models [ 12 , 16 ] are much peculiar to their identified localities. ii Since dead patients do not transmit the virus, it implies that they do not contribute to infection dynamics. Then, needless for the inclusion of dead compartment among the system model [16] . iii Infectious model that focuses only on varying exposed and infected variables may not give a true representation of disease transmission dynamics [12] . iv Lack of natural source rate in infectious model could lead to abrupt experimental result and untimely population extinction [16] . v Treatment interventions/control measures were not parametrically identified or valued for models [ 12 , 16 ] . vi The model [12] is optimal in nature (being treated as functions of time variation) but the study was not optimally analyzed. It is on this premise that we proposed what could be considered a standard mathematical COVID-19 dynamic model as anticipated in the next sub-section. Following the above aforementioned limitations, the present study in an attempt to formulate a broader COVID-19 dynamic model determined by specific state variables and parameters, is further guided by the following assumptions. i Only the infectives die due to infection, such that α 1 α 2 α 3 α 4 α 5 for all α i (i =1 ,..., 5) > 0 . ii Determination of aware infective is by screening method or any other clinical technique i.e., θ ≥ 0 . iii Contact rate of susceptible with super-spreaders is much less than isolated aware infectives, which in turn is much less than hospitalized infectives and subsequently much less than aware infected and much lesser than unaware asymptotic population (i.e., β 5 < β 4 < β 3 < β 2 < β 1 ). iv Only the hospitalized and isolated aware infectives use pharmaceutical control functions i.e., ( a 1 > a 2 0) . v Latent infection period and non-cytotoxic carrying process are ignored. vi Age-structure in transmission is ignored. vii Recovery are recruited to susceptible population i.e., σ 1 , σ 2 > 0 . Furthermore, suppose the population under study is denoted by N(t ) with population volume measure in cel l s/m l 3 and such that by subdividing the population, we let the S p (t) represent susceptible population who are not COVID-19 positive but may be infected if completely exposed, X p (t) depicts the exposed class, A u (t) defining the unaware asymptomatic infectious population, I a (t) -number of COVID-19 aware infectives, I s (t) -isolated infectious population, S s (t) -population of super spreaders, H i (t) -population hospitalized infectives and R p (t) representing recovered population, then by the existing limitations and assumption 1, the structured epidemiological differential equations of the present study is derived as: where If β i > 0 such that our control functions u i = 0 a i = 0 , where i = 1 , 2 , then system (1) completely represent a standard endemic COVID-19 infection dynamic model. Furthermore, the epidemiological descriptions of model (1) are as follows: from the first equation, the first three terms b p , σ 1 X p , σ 2 R p represent the birth rate, recovery rate of exposed and the recovery rate from varying severe forms of the infection, which reunite with the susceptible population. The fourth term β i ( ˆ N ) S p depicts the rate at which the susceptible interacts with the varying infectious compartments, which then serves as supply rate for the exposed class. The susceptible die at a rate μS p . The second term (1 − u 1 ) λX p from the second equation define the population of exposed that successively use face mask yet becomes asymptomatically infectious, while the last term (μ + σ 1 ) X p describe the sum of exposed that recovers due to coherent use of non-pharmaceutical control (face mask) measure couple with strong natural immune defense mechanism and natural clearance rate. From the third equation, the first term (1 − u 1 ) λX p denotes the transmutation of the exposed under face mask for unaware asymptomatic class, while the second term (1 − u 2 ) kθ A u is the proportional rate of unaware asymptomatic that are screen to become aware infectives with promising social distancing. The third term (1 − a 1 )(1 − a 2 ) ϕ 1 A u represent the proportion of asymptomatic infectious that progress to be hospitalized under pharmaceutical control functions, while the last term (μ + α 1 + ϕ 2 ) A u depicts the clearance differential sum from the asymptomatic infectious patients arising from natural and infection death rate as well as the proportion that becomes super-spreaders. In equation four, the first term (1 − u 2 ) kθ A u supply the source rate from asymptotic that becomes aware following the application of screening method. The second term ] I a defined the sum differential of aware Table 1 is extracted and modified from models [ 12 , 16 ] . infectives that move to H i and I s where they are subjected to treatment functions while the proportion (1 − ρ 1 − ρ 2 ) I a becomes superspreaders. The last term here depicts the clearance rates due to nature and infection respectively. The proportions of aware infectives and super-spreaders a 1 τ 1 ρ 1 I a , a 1 τ 2 γ s S s under initial treatment function are seen in the first and second terms of the fifth equation of the isolated infective compartment. The third term [(1 − a 1 )(1 − a 2 ) δ h ] I s represent severe isolated patients that progresses to hospitalization under multi-therapies. The last term (μ + α 4 + η 2 ) I s defines the removal rates arising from both natural death, death due to infection and a proportion that recovered from the infection. From equation six, the proportion of aware infectives that transmutes to super-spreaders is given by the first term (1 − ρ 1 − ρ 2 ) I a , which is busted by proportion of asymptomatic infective, ϕ 2 A u under no treatment function. The rate at which super-spreaders becomes isolated and then place under initial treatment is given by a 1 τ 2 γ s S s . The clearance rates due to natural death, infection rate and recovery rate are given by (μ Taking on equation seven, the first term (1 − a 1 )(1 − a 2 )[ ϕ 1 A u + ρ 2 I a + δ h I s ] describe the varying rates of severe infected population that are hospitalized under multi-therapies. The last term in this compartment (μ + α 3 + η 1 ) H i depicts the sum clearance rates due to deaths (natural and infection) and recovery rate. The final and eight compartment define the sum recovery population coming from hospitalized compartment, recovery from isolated compartment and from super-spreaders denoted by η 1 H i + η 2 I s + η 3 S s . Here, clearance rate is due to natural death rate and the proliferation of recovered to the susceptible denoted by (μ + σ 2 ) R p . Thus, following the above model description, model (1) can be schematically represented as in Fig. 1 , below: A critical review of system (1) shows that the present model is characterized by system force of infection rate depicted by Eq. (2) . This force of infection gives a biological meaning to the system reproduction number, (see Section 4 ). Moreso, the design of present model explicitly outline the methodological application of designated control functions following the introduction of dual pharmaceuticals at severe stages of COVID-19 infection. Furthermore, we perform the numerical simulations using generated initial values for the state-space and parameters of the system in relation to model (1) and Fig. 1 . That is, with modified data from related compactible models, Tables 1 and 2 below depicts the simulating data for the present investigation. Obviously, since model (1) represent a set of living organism, then we must verify the model well-posedness, which involves the mathematical properties of derived system. Rate at which I a progresses to I s and S s 0.3, 0.5 c i (i =1 ,.., 5) Rates of contact of susceptible with various infectious stages 0.5;0.4;0.3;0.2;0.1 Rates of recovery from H i , I s and S s 0.5; 0.27;0.13 Probability of interactions of susceptible with varying infectious classes 0.32;0.27;0.175; 0.125;0.05 Proportions of A u that progresses to H i and S s 0.3;0.18 Proliferations of recovered population to susceptible 0.14;0 Proportion of I a that progresses to S s 0.08 Rates of use of face mask and social distancing Treatment control functions (HCQ and AZT) Note: Tables 2 is clinically generated from certified data of [ 9 , 16 ] . Here, we quantitatively verify the characteristic properties of the model constituted by the system invariant region and non-negativity of system solutions. Notably, if the sum population under study is defined by Eq. (3) , then the differential sum of model (1) is given by where ˆ N = ( X p + A u + I a + S s + H i ) and N(t ) defined by Eq. (3) . Clearly, Eq. (4) shows that the differential sum of system (1) is a function of the system birth rate, natural death rate and death due to infection. We also use Eq. (4) as the index for the verification of the system invariant region. For the fact that the system represents a set of living organism, then it is sufficient to assume that for t > 0 , the state variables and parameters are non-negative at any given region. This region we shall investigate using the following theorem. Theorem 1. Let the system state variables be bounded in a closed set D such that Then, the closed set D is positively invariant and attracting with respect to system (1). We invoke results from [ 21 , 22 ] , then by Eq. (6) , we have In the absence of mortality due to COVID-19 infection, the population is said to be free from the virus i.e., ˆ which gives a first order homogeneous differential inequality. The integration in the presence of its initial conditions yields Then, in line with Birkhof and Rota's theorem on differential inequality as applied by [23] , when t → ∞ , we arrive at 0 ≤ N(t ) ≤ b p μ for all t ≥ 0 . But from model (1), The resulting integral gives where C, is the constant of integration. We know from Eq. (3) that It follows that C = 1 , implying that the population is constant, positive and equal 1. Hence, all the feasible solutions of the system (1) enter the invariant region Therefore, the region is positive and attracting implying that the model is both mathematically well-posed in the region D and epidemiologically meaningful. The following theorem is use to show that the solutions of the system state variables remain positive for all t > 0 Proof. Invoking existing results from [ 22 , 24 ] , then system (1) can be confined in compact subset as ) be any solution with positive initial conditions such that Then, the time derivative of N(t ) along solution of system (1) with zero mortality rate is dN dt Applying the theorem of differential equation, we obtain Clearly, it has been proved that all the solutions of system (1), which is initiated in 8 + is confined in the region , implying the solutions are bounded and positive in the interval [0 , ∞ ) . In this section, we quantitatively investigate the model equilibria and explicitly study the model local and global stability conditions. For the existence of an equilibrium state, it is assumed that model (1) is at its steady state i.e., This gives the system Then, from Eq. (5) , the COVID-19 free equilibrium (C-19FE) for system (1) exists if u i = 0 , a i = 0 for all i = 1 , 2 with other controls held constant. Therefore, in computing the C-19FE, we let E 0 denotes C-19FE such that at C-19FE, there is no infection and thus, no recovery i.e., Then, we define We solve Eq. (7) using Eqs. (5) and (6) to yield . From the second equation, substituting S * p we have . In a similar approach, we solve for A * u , I * a , I * s , S * s , H * i and R * p , which gives the comprehensive result for E 0 as: where, Now, since β * i = 0 define the point at which C-19FE exists, then from Eqs. (5) and (8) , the equilibrium point given by Eq. (7) is obtained as: Clearly, Eq. (10) depicts C-19FE, where no infection exists. But with the introduction of COVID-19 into the system, we are the required to investigate the pattern of transmission dynamics. This process invariably demands that we establish the system reproduction number denoted by 0 . Focusing on infectious diseases, the basic reproduction number is a critical mathematical quantity considered paramount to the public health sector. For a typical infectious disease like the COVID-19, we defined the reproduction number as the average number of new cases of the infection generated by an infected individual (i.e., super-spreaders) following the interaction with the susceptible population. Here, we adopt the approach by [25] , which explore the Next generation matrix defined by where the notations F i and V i , denotes the matrices of new infections in compartment i and the transfer terms at COVID-19 free equilibrium into compartment i while E 0 is the C-19FE for finding 0 . Now, since reproduction number is about the rate of infection, then we required to rewrite system (1) starting with the exposed class and substituting Eq. (9) to have, Substituting Eq. (10) , we have Computing for V i , we have Applying the linearization method, we have Thus, the COVID-19 reproduction number as applied by [ 26 , 27 ] corresponds to the spectral radius of where j = 1 , ..., 5 represent the reproduction numbers for the infectious state variables defined in Here, we explore the eigenvalues of the linearized Jacobian matrix, which is evaluated in terms of the negative trace and a positive determinant or as having negative eigenvalues [28] . Of note, Eq. (14) is significant in the sense that it affirm the existence of Eq. (10) for a COVID-19 free population provided 0 < 1 . Moreso, the endemicity of COVID-19 can be effectively control, if ( β i c i ) < 1 and under significant φ. A situation that will lead to reduction in 0 and subsequently, the elimination of COVID-19 from the population. The following theorem justify the role of 0 in COVID-19 transmission. Theorem 3. Whenever the C-19FE for model (1) exists, it is locally asymptomatically stable provided 0 < 1 and unstable if 0 > 1 . Proof. Results of existing theorem by [29] is invoke for this prove. Let J be the Jacobian matrix for the system (1). Then, J at C-19FE point ( E 0 ) is derive as: 8) , the differential sum of the system (1) at C-19EE, is derived as: [24] b which corresponds to the fact that at equilibrium, β * i = 0 . If β * i > 0 , then there exists disease endemicity, which in terms of 0 , we have provided 0 > 1 with ˆ as the disease constant derived from Eqs (9) to (13) . Hence, proof completed. We shall consider in this section the global stability of C-19FE and that of C-19EE. The following notations are necessary. with G ( X 1 , X 2 ) = 0 , where X 1 ∈ m denotes the uninfected population and X 2 ∈ n denotes the infection population; X 0 = (X E 1 , 0) denotes the C-19EE of the system. Also, assume that the following conditions holds: 0) , X E 1 is globally asymptomatically stable. has all non-negative off-diagonal elements and X is the region where the model makes biological meaning. Then, the C-19EE, X 0 = (X E 1 , 0) is globally asymptomatically stable provided that 0 < 1 and unstable otherwise. We prove our global stability by invoking the theorem from [30] . 0 , 0 , 0 , 0 , 0 , 0 , 0) is globally asymptotically stable equilibrium for the system (1) provided 0 < 1 with Lemma 2 satisfied. (1), then which has the solution Clearly, S p (t) → b p μ as t → ∞ regardless of the initial value S P (0) . Therefore, it shows that condition C 1 (of Lemma 2 ) holds for our model. Taking on the infectious sub-systems and using Eq. (13) with condition C 2 (of Lemma 2 ), we have and ω 5 = φβ 5 c 5 b p μN with 1 , 2 , 4 , 8 and m 6 defined by Eqs (13) and (9) respectively. This implies that We also notice that the matrix A is an M-matrix since its off-diagonal elements are non-negative. Hence, this proves the global stability of C-19FE ( E 0 ). Here, with the incorporation of the LaSalle's invariant principle, we shall adopt the Lyapunov function for the investigation of the global stability of system endemic equilibrium (C-19EE). This approach, which had been explored by [ 24 , 31 ] is found useful in compartmental epidemic models. Then, the following theorem establishes the system global stability of C-19EE. Theorem 6. Let V be the Lyapunov function defined for system (1), then the global stability of the system endemic equilibrium (C-19EE) holds if its time derivative dL dt ≤ 0 . which is clearly positively invariant in 8 + . We also noted that at ˆ Then, to prove for the global stability result, we construct the following Lyapunov function where w i > 0 is a Lyapunov constant, N i is the population of i th compartment, N * i is the equilibrium value of N i and V , a continuous and differentiable Lyapunov function. Computing the time derivative of V , along the trajectories of the system (1), we obtain Substituting the derivative of system (1) into Eq. (21) and accounting for Eq. (9) , we have Or Then, Eq. (23) can be written in compact form as: where Now, we discuss the global asymptotic stability of Q by showing that the matrix Y of Eq. (25) is Lyapunov stable or −Y is diagonally stable. The following lemmas and theorem yields the required proof. Therefore, using the LaSalle's invariant principle, the proof for the global stability of the system endemic equilibrium is completed from the following theorem. (1) is globally asymptotically stable in E 0 , provided 0 > 1 and unstable otherwise. Proof. Invoking Lemmas 5 , 6 and Theorem 7 , we obtain dV dt < 0 when E 0 = E * and E 0 is not on S − axis (a set of measure zero). Therefore, the largest invariant set in E 0 such that dV dt < 0 is a singleton E * , which is our endemic equilibrium point. Then, by LaSalle's invariant principle [32] , it follows that the endemic equilibrium of model (1), E * is globally asymptotically stable (GAS) in D , if E 0 = E * . This complete the proof. Here, attempt is made to verify the viability of derived theoretical predictions. That is, we numerically illustrate the system force of infection β i ( ˆ N ) under 0(1) = 10 . 94 . Clearly, the inclusion of both system force of infection and the reproduction number in system (1) is an integral component of the model novelty. These two key components determine the dynamic behavior of the system (1) at u i = 0 , a i = 0 for all i = 1 , 2 , which we shall simulate. Fi-nally, we simulate the system stability endemic equilibrium following the application of control functions u i > 0 , a i > 0 for all i = 1 , 2 under reproduction number at onset of treatment computed to be 0(2) = 3 . 224 . The entire simulations is feasible under in-built Runge-Kutta software in a Mathcad environment. Invoking Tables 1 , 2 and Eq. (2) , we illustrate the system force of infection β i ( ˆ N ) against the susceptible S p (t) and the exposed X p (t) components. That is, Fig. 2 (a)-(b) depicts the dynamics of the system force of infection for untreated COVID-19 scenario. Notably, Fig. 2 (a)-(b) portrait the consequential rate of force of infection on the susceptible and the exposed by the infectious classes denoted by β i ( ˆ N ) . Clearly, for an untreated COVID-19 dynamics, which is vindicated by high system reproduction number 0(1) = 10 . 94 , the virulence ingress required to aggressively contaminate both the susceptible ( 0 . 5 ≤ S p (t) ≤ 5 . 335 × 10 12 cel l /m l 3 ) and exposed sub- for all t f ≤ 90 days . This explains why the practical exponent spread of the virus is observed under off-treatment scenario world-wide. Furthermore, the consequential of system force of infection on varying subpopulations observed under offtreatment is seen in the next Section 5.2 . The impact of both the force of infection and system reproduction number under off-treatment protocol is explicitly illustrated as depicted by Fig. 3 (a) days. This rapid decline is evident by the exponential increase in the rate at which population become exposed with value at 0 . 3 ≤ X p (t) ≤ 49 . 785 cel l s/m l 3 for all t f ≤ 18 days and then decline slightly to stability at X p (t) ≤ 41 . 024 cel l s/m l 3 after 18 ≤ t f ≤ 90 days (see Fig. 3 (b) ). Fig. 3 (c) depicts rapid increase in the population that becomes unaware asymptomatic infectious class i.e., 0 . 1 ≤ A u (t) ≤ 11 . 675 cel l s/m l 3 for the first 10days and then decline slightly to stability at A u (t) ≤ 9 . 8 cel l s/m l 3 for all 20 ≤ t f ≤ 90 days. Fig. 2. (a-b) Dynamics of model force of infection against S p (t) , X p (t) with 0(1) = 10 . 94 . In Fig. 3 (d) , the screen aware infectives is simulated. We observe unsteady slow inclination curve of 0 . 15 ≤ I a (t) ≤ 1 . 312 cel l s/m l 3 but far lower compared to Fig. 3 (c) . From Fig. 3 (e) , following the non-availability of treatment measures, it is clear that isolation compartment, which constitute control/treatment measure remain at zero i.e., I s (t) = 0 . 0 for all t f ≤ 90 days. Of interest, the system is bound to experience commensurate increase in the proportion of super spreaders under off-treatment scenario. This can be attributed to varying natural adaptive immune response. Having theoretically investigated the endemicity for the spread of COVID-19 viral load and the methodological application of dualbilinear control functions, we devote this sub-section to the illustration of the viability of our findings. Obviously, with the introduction of pair non-pharmaceutical control functions (face masks and social distancing) and pair pharmaceutical intervention functions (hydroxyl-chloroquine -HCQ and Azithromycin -AZT), we compute as represented by Fig. 4 (a) days. Notably, the introduction of non-pharmaceutical preventive measures as in Fig. 4 (b) , shows decline in the rate of exposed class with value at 0 . Fig. 3 (F) . Finally, following the significant decrease in the rate of isolation, super spreaders and hospitalization, the proportion of recovery equally decline to stability R p ≤ 0 . 047 cel l s/m l 3 for all 40 ≤ t f ≤ 90 days -see Fig. 4 (h) . Of note, the choice of Lyapunov function in conjunction with LaSalle's invariant principle in analyzing the global stability conditions of design model have clearly given insight to COVID-19 infection dynamics and as well, portrait a novel methodological control approach in the quest to mitigate COVID-19 pandemic. These results provides unique treatment approach when compare with results of system motivating models [ 12 , 16 ] , where the aforementioned technique were not considered. Moreso, the method implemented by this study have affirm the high effectiveness of the combination of designated non-pharmaceutical and pharmacotherapies in controlling the deadly disease. The study by [16] had been formulated as an adhoc compartmental COVID-19 model to account for the peculiarities of COVID-19 infection dynamics in Wuhan, China, while the study [12] considered contact parameter as its cardinal point in its mathematical formulation of COVID-19 transmission in Niger Republic. On account of earlier highlighted limitations of these two models, which borders on non-availability of mathematical model for the combination of medical pharmaceuticals and non-pharmaceuticals in the treatment protocols of spiking COVID-19 pandemic, the present study sought to determine and analyze the global stability conditions that could lead to predictions of COVID-19 infection dynamics and treatment control protocols. The study using nonlinear differential equations, formulated its model as an 8-Dimensional deterministic compartmental mathematical COVID-19 dynamic model involving the interactions between healthy and COVID-19 shrouded infectious population. Designated control functions are categorized into dual non-pharmaceuticals (face-masking and social distancing) and dual pharmacotherapies (hydroxylchloroquine and azithromycin). The study adopted classical method of Lyapunov function in conjunction with LaSalle's invariant principle for the analysis of the system global stability conditions. The investigation was conducted for both untreated and onset treatment scenarios and numerical validations performed using in-built Runge-Kutta of order of precision 4 in a Mathcad surface. Results obtained indicated that under off-treatment scenario, the logistic dose response curve of heavy droplet of finer aerosol viral load of COVID-19 required to extensively contaminate the susceptible population of S p (t) ≤ 5 . 335 × 10 12 cel l s/m l 3 is in the range of β i ( ˆ N ) ≤ 2 . 12 × 10 11 under exponential reproduction number of 0(1) = 10 . 94 . This result clearly buttress why within short space of time, the entire world was highly ravaged by the deadly virus. Moreso, under offtreatment scenario as depicted by Fig. 3 (a) -(h), we observed rapid extinction of the susceptible population leading to excruciating endemic spread of the virus for all 18 ≤ t f ≤ 90 days. Following the introduction of methodological control approach under dual-bilinear control functions, the system reproduction number declined drastically to 0(2) = 3 . 224 . This result is in agreement with COVID-19 estimated reproduction number (2-3) by [1] . Moreso, we found from Fig. 4 (a)-(h) that the application of only non-pharmaceuticals at onset of infection (mild) indicated immense reduction in the spread of the virus, while for severe cases, the application of non-pharmaceuticals enhanced by induced multi-therapies lead to significant reduction in the rate of isolations, super spreads and hospitalized population with attained stability of ( I s (t) ≤ 5 . 8 × 10 −3 , S s (t) ≤ 0 . 061 , H i (t) ≤ 0 . 031 ) cel l s/m l 3 . Of note, the positive impact from dual-bilinear treatment protocols is vindicated by the enhanced reduction in the rate of recovery and the accelerated rejuvenation of the susceptible population as against population extinction under off-treatment scenario. That is, the application of the method of Lyapunov function incorporating LaSalle's invariant principle and methodological dual bilinear control protocols has a very desirable effect upon the susceptible population. Comparatively, the present results when compared with those of system motivating models [ 12 , 16 ] , shows that our controls does behaves somewhat different from control functions used in motivating models not explicitly studied under dualbilinear protocols. Following the non-availability of mathematical model for the combination of medical pharmaceuticals and non-pharmaceuticals in the treatment protocols of spiking COVID-19 pandemic, the present study sought to determine and analyze the global stability conditions for the role of dual-bilinear control functions in the control and treatment dynamics of novel COVID-19 pandemic. The model adopted human-to-human transmission mode with population under-consideration partitioned into 8-Dimensional deterministic compartments: Susceptible population S p (t) , Exposed class X p (t) , Unaware asymptomatic infectious population A u (t) Aware infectives I a (t) , isolated infectious population I s (t) , Super spreaders S s (t) , Hospitalized infectives H i (t) and the recovered population R p (t) . First, we investigated the model state-space and established the system reproduction number for both off/on treatment protocols using modified generated data from certified models. Using classical method of Lyapunov functions in combination with the theory of LaSalle's invariant principle, we have discussed the global stability conditions of COVID-19 model. Obviously, the method of Lyapunov function has been widely applied to varying dynamical systems, but the essential part of this analysis is based on the incorporation of the theory of LaSalle's invariant principle under dual-bilinear control protocols. The analytical predictions of the global stability analysis are provided and numerically simulated. We found that under designated control functions, COVID-19 transmission was drastically reduced to insignificant threshold of ( 0 . 1 ≤ A u (t) ≤ 0 . 247 , 0 . 15 ≥ I a (t) ≥ 0 . 01 , 0 . 05 ≤ S s (t) ≤ 0 . 061 and 0 . 01 ≤ H i (t) ≤ 0 . 031 ) cel l s/m l 3 for all 18 ≤ t f ≤ 90 days, leading to tremendous rejuvenation of the susceptible population, 0 . 5 ≤ S p (t) ≤ 3 . 143 cel l s/m l 3 . It is presumed that the insignificant persistence of the virus can be attributed to possible reinfection of recovered population upon integrated into the susceptible population. Furthermore, the present results in comparison with those of motivating models, have projected classical methodological and epidemiological concept under designated clinical conditions. Thus, for enhance optimal result, the study highly suggest possible application of optimal control theory to the existing model. Bassey Echeng Bassey: Conceptualization, formulations, writing, investigation, editing, review, programing and analysis. Jerimiah U. Atsu: Supervision, methodology, funding, editing and validation. This undertaking certify that onbehalf of the authors, the corresponding author wish to state clearly that there exist no conflict of interest in the submission of this manuscript. I also would like to declare that the work described was an original research that has not been published previously and not under any consideration for publication elsewhere. Using mathematics to understand and control the coronavirus pandemic Mathematical assessment of the impact of non-pharmaceutical interventions on curtailing the 2019 novel Coronavirus T cells found in coronavirus patients 'bode well' for long-term immunity World Health Organization. WHO Coronavirus disease (COVID-19) Dashboard MW2rjYNlHlaLtHevv9rld75E0 _ 61mSQ-pEL-26Netr _ MaAoexEALw _ wcB Nigeria Centre for Disease Control. Coronavirus disease (COVID-19) pandemic Optimal quarantine strategies for COVID-19 control models COVID-19 research in Africa Recognition of COVID-19 disease from X-ray images by hybrid model consisting of 2D curvelet transform, chaotic salp sqram algorithm and deep learning technique A dynamical behavior of COVID-19 virus model with carrier effect to outbreak epidemic Global stability and cost-effectiveness analysis of COVID-19 considering the impact of the environment: using data from Ghana Mathematical modeling and analysis of COVID-19 pandemic in Nigeria Contribution to the mathematical modeling of COVID-19 in Niger Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts Early dynamics of transmission and control of COVID-19: a mathematical modelling study The reproductive number of COVID-19 is higher compared to SARS Coronavirus Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan Will an imperfect vaccine curtail the COVID-19 pandemic in the U.S A novel coronavirus from patients with pneumonia in China Modeling of a super-spreading event of the MERS-corona virus during the Hajj season using simulation of the existing data A pneumonia outbreak associated with a new coronavirus of probable bat origin Stability analysis and modeling of listeriosis dynamics in human and animal population Dynamic optimal control for multi-chemotherapy treatment of dual listeriosis infection in human and animal population Invariant region, endemic equilibria and stability analysis The Volterra-Lyapunov matrix theory for global stability analysis of a model of the HIV/AIDS Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Age-structured malaria transmission model A mathematical model for the dynamics of cholera with control measures Stability analysis of the dynamics of tungiasis transmission in endemic areas Modeling and analysis of cholera dynamics with vaccination On the computation of R 0 and its role on global stability in mathematical approaches for emerging and re-emerging infectious disease: an introduction. IMA Volumes in Math Approaches Global properties of basic virus dynamics models The stability of dynamical systems. CBMS-NSF regional conference series in applied mathematics