A non-linear mathematical model for analysing the impact of COVID-19 disease on higher education in developing countries

This study proposes a non-linear mathematical model for analysing the effect of COVID-19 dynamics on the student population in higher education institutions. The theory of positivity and boundedness of solution is used to investigate the well-posedness of the model. The disease-free equilibrium solution is examined analytically. The next-generation operator method calculates the basic reproduction number (R0). Sensitivity analyses are carried out to determine the relative importance of the model parameters in spreading COVID-19. In light of the sensitivity analysis results, the model is further extended to an optimal control problem by introducing four time-dependent control variables: personal protective measures, quarantine (or self-isolation), treatment, and management measures to mitigate the community spread of COVID-19 in the population. Simulations evaluate the effects of different combinations of the control variables in minimizing COVID-19 infection. Moreover, a cost-effectiveness analysis is conducted to ascertain the most effective and least expensive strategy for preventing and controlling the spread of COVID-19 with limited resources in the student population.


Introduction
The first human cases of coronavirus disease 2019 (COVID- 19), the disease caused by the novel coronavirus, subsequently named SARS-CoV-2 were first reported by officials in Wuhan City, China, in December 2019 [1,2]. From all available investigation, the ecological origin of SARS-CoV-2 is in bat populations. All available evidence to date suggests that the virus has a natural animal origin and is not a manipulated or constructed virus [2]. From all indication, the source of COVID-19 suggests that SARS-CoV-2 has a zoonotic source [2,3]. It is not likely that bat is the source of transmitting it to human, since there exists limited close contact between humans and bats, more likely that transmission of the virus to humans happened through another animal species, one that is more likely to be handled by humans [1,2,4]. This intermediate animal host or zoonotic source could be a domestic animal, a wild animal, or a domesticated wild animal and, as of now, has not been identified [3,4]. COVID-19 is transmitted through direct contact with respiratory droplets of an infected person (generated through coughing and sneezing). Individuals can also be infected from and touching surfaces contaminated with the virus and touching their face (e.g., eyes, nose, mouth) [5]. COVID-19 virus may survive on used for the estimation of the final size of the coronavirus epidemic. In [21], it was reported that isolation of the infected human overall can reduce the risk of future COVID-19 spread. In addition, the model of [21] shows that the coronavirus can spread through contact, and they described how fast thing changes by counting the number of people who are infected and the likelihood of new infections. The authors in [21] recommended that state and territory governments have different restrictions in place for public gatherings, and citizens need to follow the directions from time to time to minimize the health risk. Tang and co-workers in [22] developed and analysed a deterministic model which incorporates quarantine and hospitalization to estimate the transmission risk of COVID-19 and its implication for public health interventions.
Mathematical modelling of infectious diseases in the presence of limited resources have been studied by different researchers [23][24][25][26], but none of them has applied it to . Similarly, to the best of our knowledge, the impact of online facilities to aid teaching in the era of COVID-19 pandemic were less or not studied in the literature using mathematical modelling approach. Therefore, this study presents a mathematical model of COVID-19 examining the impact of COVID-19 disease on tertiary education in developing countries in the presence of limited healthcare facilities and online facilities. To the best of our knowledge, this is the first mathematical modelling paper presenting COVID-19 with these features. The rest of this study is sectioned as follows: Section 2 discusses the model formulation. Analysis of the autonomous version of the model is carried out in Section 3. Formulation and analysis of the non-autonomous version of the model is carried out in Section 4. In addition, the numerical simulations of both the autonomous and non-autonomous models as well as the costeffectiveness analysis are performed in Section 5. Section 6 presents the concluding remarks of the study.

Model formulation
This section considers the formulation of an epidemic model which focuses on the impact of COVID-19 disease on tertiary education. Thus, the host (that is, student) population, denoted as ( ), is stratified into nine mutually exclusive compartments, namely, susceptible class, ( ) (students who are prone to contract COVID-19 infection), exposed class, ( ) (students who have been infected, but not infectious yet), quarantined class, ( ) (those who are in self-quarantine due to the imposed lockdown), infectious class, ( ) (students who exhibit the symptoms of COVID-19 and are capable of the disease spread), online facilities-accessed class, ( ) (those who have access to the online facilities), online facilities-unaccessed class, ( ) (those who cannot access the online facilities), learned class, ( ) (those who have learned through the available online facilities), hospitalized class, ( ) (those who are infectious and admitted to healthcare facility) and recovered class, ( ) (individual students who have recovered from COVID-19 infection). Thus,

( ) = ( ) + ( ) + ( ) + ( ) + ( ) + ( ) + ( ) + ( ) + ( ). (1)
In order to understand the effect of inadequate online facilities (such as Language learning apps, online tutoring platforms, virtual tutoring, video conferencing tools, or online learning software and others) in connection with transmission dynamics of COVID-19, we introduce a learning rate, denoted as , which incorporates the effect of capacity and inadequate online facilities in terms of internet facilities-population ratio, represented by > 0. Also, depends on the number of students who have access to internet ( ), such that is a function of both and . Obviously, the per capita learning rate, ( , ), is an increasing function of and a decreasing function of , which is bounded above and below by respective 0 and 1 for every > 0. In this case, 1 represents the maximum per capita learning rate due to the sufficient online facilities and large students with access to the online facilities, while 0 is the minimum per capita learning rate due to the function of basic online facilities.
Based on the above assumption, the impact of online facilities on the learning rate is modelled as It implied from (2) that as tends to zero, tends to 0 , and as tends to infinity, tends to 1 .
In a similar manner, the recovery rate for individuals in class incorporating the impact of capacity and inadequate health facilities is formulated as a nonlinear function as where > 0 represents the hospital bed-population ratio, 0 is the minimum per capita recovery rate owing to the function of basic healthcare system, and 1 accounts for the maximum per capita recovery rate corresponding to the sufficient clinical resources and few hospitalized humans. A similar function to the ones given in (2) and (3) was previously suggested in [23]. Consequently, the model governing the COVID-19 dynamics in the student population according to the flowchart illustrated in Fig. 1 is derived as with initial conditions given as where = ( , , , , , , , , ). The definitions of the parameters used in model (4) are given in Table 1, while the schematic diagram of the model is as shown in Fig. 1.

Positivity of solutions
We now consider the positivity of model (4). We are to prove that all the state variables remain non-negative and the solutions of model (4) with positive initial conditions remain positive for all time > 0. We thus state the following theorem.

Invariant region
From (6) we have ≤ 0 which implies that is a positively invariant set. We also note that by solving (6) we have where (0) is the initial condition of ( ). Thus, 0 ≤ ( ) ≤ 0 as → ∞ and hence is an attractive set. □

Disease-free equilibrium and the basic reproduction number
The COVID-19 model (4) has a disease-free equilibrium given by  0 = ( * , * , * , * , * , * , * , * , * ) = ( 0 , 0, 0, 0, 0, 0, 0, 0, 0 a scenario depicting an infection-free state in the society. The basic reproduction number,  0 , defined as the expected number of secondary cases produced by a single infectious individual in a completely susceptible population over the duration of its infectious period, is a threshold quantity that allows us to predict whether the disease will die out or persist in the population [29]. Generally,  0 < 1 means that the disease (COVID-19) cannot invade the population,  0 > 1 means that each infected individual produces more than one secondary infected individual and  0 = 1 requires further investigation. The determination of  0 is done using the next generation matrix approach [23,27,[29][30][31].
It is worthy of note that at the disease-free equilibrium,  0 , ( , ) = 0 and ( , ) = 1 . Thus, the  0 of model (4) is given by Sensitivity index of  0 to each of the parameters of model (4) evaluated at the baseline parameter values given in Table 3. COVID-19 disease can be eliminated from the population (when  0 < 1) if the initial sizes of the population of the model are in the basin of attraction of the disease-free equilibrium. This will be established by the following Theorem 3.3.

Theorem 3.3.
The disease-free equilibrium  0 is locally asymptotically stable when  0 < 1 and is unstable when  0 > 1.

Sensitivity analysis
In this section, sensitivity analysis of model (4) is explored. To do this, we follow the approach in previous studies [27,[32][33][34]. This helps to identify the parameters of the model that most influence the dynamical behaviour of the model. Thus, to determine the sensitivity of model (4) using the COVID-19 threshold quantity, the basic reproduction number  0 given in (10), as the response function with respect to the model parameters, the normalized forward-sensitivity index of  0 that depends on a parameter is defined as Given the explicit formula (10) for the basic reproduction number  0 , the analytical expressions for the sensitivity of  0 in respect of the parameters defining it are computed. In particular, the analytical expression for the sensitivity of  0 with respect to in view of (11) is a constant value, while that of is a complex expression, and are both obtained as ( ℎ 1 ( 1 + 2 + ℎ 2 + + 0 ) + ℎ 2 (ℎ 1 + 0 + 1 )) ( 1 + 2 + ℎ 2 + + 0 )( ℎ 1 + 0 + 1 + 1 ) + ℎ 2 (ℎ 1 + 0 + 1 ) .
In a similar manner, the analytical sensitivity indices of  0 for the other parameters are computed. But the results are omitted due to their complexity. However, the sensitivity index (SI) of  0 for all the parameters comprising it are evaluated at the baseline parameter values given in Table 3. The signs and values of SI are presented in Table 2.
From Table 2, it is observed that the sign of SI is positive for some parameters ( , , ℎ 2 , ), while it is negative for the others ( , ℎ 1 , 0 , 1 , 2 , , 1 , 1 ). From the set of parameters with positive SI sign, , and are most positive. Whereas, 1 , and ℎ 1 are most negative from the set of parameters with negative SI signs. The epidemiological insight from the positive sign of SI of the COVID-19 threshold quantity,  0 , is that increasing or decreasing the value of any of the parameters in this category will generate an increase or decrease in the threshold  0 of COVID-19. The negative sign of SI on the contrary suggests that increasing the value of each of the parameter set in this category will lead to a decrease in the  0 value, and vice-versa. For example,  0 = +1 indicates that increasing the effective transmission rate of   Table 2, considerations of any control strategies that reduce the effective transmission rate of COVID-19 ( ), modification parameter relative to the infectiousness of hospitalized students ( ), progression rate from exposed to become symptomatic ( ), and increase the maxima per capital recovery rate associated with the sufficient healthcare resources and few hospitalized humans ( 1 ), detection rate for exposed students for isolation or self-quarantine ( ) and hospitalization rate for symptomatic infectious students (ℎ 1 ) are needed to ensure an effective control of COVID-19 transmission and spread in the student population.
Furthermore, Latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCC) are used to carry out a global sensitivity analysis to further ascertain the parameters that mostly influence the dynamics of model (4) with respect to the basic reproduction number,  0 , given in (10). LHS matrices are generated by assuming all model parameters obey a uniform distribution. Then, using a similar idea of the authors in [37,38], multiple runs of 1000 sample size per LHS run are performed based on the model parameter baseline values provided in Table 3 with the ranges of 30% in either direction from the baseline values. Thus, the sensitivity of model parameters (4) is evaluated by finding the PRCC between each parameter and  0 .  (4). It is observed from Fig. 2 that the PRCC value is positive for some parameters, that is, those that are positively correlated ( , , ℎ 2 and ) and negative for others, that is, those that are negatively correlated ( , ℎ 1 , 0 , 1 , 2 , , 1 and 1 ). The top PRCC-ranked parameters that drive the response function,  0 , most are the effective transmission rate, , maximum per capita recovery rate corresponding to the sufficient clinical resources and few hospitalized individuals, 1 , detection rate for exposed individuals, , modification parameter for infectiousness of hospitalized individuals, , and rate of infectious of exposed individuals, (see Fig. 2 and Table 2). Therefore, the result from global sensitivity analysis aligns with the result arising from the local sensitivity analysis.

Optimal control problem formulation and analysis
In view of the results obtained from the sensitivity analysis carried out on the parameters of COVID-19 model (4), we apply optimal control theory to the model so as to hamper the community spread of the disease optimally. To do this, four time-dependent control variables (i) 1 ( ) represents the personal protection strategy (precautionary measures) introduced to inhibit the transmission of the virus by symptomatic and hospitalized students. Compliance to different precautionary measures such as social distancing, good personal hygiene, wearing face masks in public places, and provision of protective gear for healthcare workers may be helpful to achieve this. Thus, the effective transmission rate in the presence of this control intervention becomes ∶= (1− 1 ( )). Notably, 1 ( ) = 1 suggests the efficacy of the strategy to protect against COVID-19 infection. On the other hand, 1 ( ) = 0 implies the failure of the strategy for protection against infection. (ii) 2 ( ) denotes the effort to scale up the contact-tracing policy and testing capacity for detection to quarantine the exposed students. Consequently, the detection rate for exposed students to become self-quarantined is modified as ∶= + 2 ( ). (iii) 3 ( ) is considered as treatment (or hospitalization) control variable, which is used to enhance the hospitalization of diagnosed COVID-19 infected cases following testing. It permits to scaleup the total number of hospitalized students and reduce the mortality rate induced by the disease for the symptomatic infectious students [39], so that the hospitalization rate of symptomatic infectious students becomes ℎ 1 ∶= ℎ 1 + 3 ( ), and the disease-induced mortality rate becomes 1 ∶= (1 − 3 ( )) 1 .
(iv) 4 ( ) represents the control variable to enhance the management of hospitalized students 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 students with severe COVID-19. It allows to increase the total number of recovered students and decrease the mortality rate induced by the disease [39]. To this aim, the maximum per capita recovery rate associated with the sufficient clinical resources and few hospitalized students becomes 1 ∶= ( 1 + 4 ( )), and the mortality rate due to COVID-19 for hospitalized students becomes 1 ∶= (1 − 4 ( )) 1 . If 4 ( ) = 1, then the control is effective in managing the disease, while 4 ( ) = 0 means the absence of control.
Consequently, incorporating the aforementioned four time-dependent control variables into the COVID-19 model (4) leads to the optimal control model given by where = ( , , 4 ( )) = 0 + ( 1 + 4 ( ) − 0 ) + , The four control variables are introduced in line with the aim of seeking the optimal solution required to minimize the sizes of symptomatic and hospitalized students subpopulations responsible for COVID-19 spread in the population at a minimum cost. Thus, the objective functional for the optimal control problem is designed as constrained by the system (12). In precise, we define The use of personal precautionary measures 1 ( ) protect both infected and susceptible humans, the contact tracing and testing policy 2 ( ) acts on the exposed individuals for detection and prompt quarantine, the treatment effort 3 ( ) acts on diagnosed infected cases following testing, while the management control 4 ( ) acts on the hospitalized students with a view to ensuring prompt recovery and prevent deaths due to complications. Thus, L 1 ( , ), L 2 ( , , ) and L 3 ( , ) appearing in (14) are written in their explicit forms as where , > 0, for = 1, 2, 3, 4, are weight constants needed for balancing the associated terms in the objective functional and is the final time such that 0 ≤ ≤ . The function L 2 ( , , ) accounts for the linear state-dependent costs associated with the personal protection control, quarantine (or self-isolation) control for exposed students, treatment control for symptomatic students and management control for hospitalized infected students, while the function L 3 ( , ) defines the quadratic state-independent costs related to the absolute or background costs such as ordering, shipment and distribution, storing among others. It is worthy to note that, unlike some previous works on mathematical models of biological processes involving optimal control (see, for example, [40][41][42][43][44][45][46][47][48][49] and some of the references cited therein), we combine linear state-dependent costs with quadratic state-independent costs (similar ideas can be found in [39,50]). The quadratic terms in (13) permit to make the Lagrangian convex with respect to the control variables ( ), = 1, 2, 3, 4, to guarantee the existence of a unique solution. Di Liddo in [51] gives more explanations about the choice of costs type in the objective functional.
The main goal is to seek an optimal control quadruple * = ( * ), where = 1, 2, 3, 4, such that where U is a non-empty bounded Lebesgue measurable control set given by

Characterization of optimal control
Pontryagin's maximum principle [54] provides the necessary conditions that an optimal control quadruple of the COVID-19 model (12) must satisfy. The principle helps in converting the COVID-19 model (12) together with the objective functional (13) into a problem of minimizing pointwise a Hamiltonian , with respect to the time-dependent controls ( ), for = 1, … , 4. The Hamiltonian for the non-autonomous system (12) is determined as where are the respective adjoint variables for the states ∈ { , , , , , , , , } , and is the right-hand side for th state of the system (12). The control characterization result is claimed in Theorem 4.2 below.

Numerical simulations
The present section explores the dynamical behaviours of COVID-19 model (4) and the derived two-point boundary problem of sixteen-  Table 3.

Autonomous system
In this model, we present the estimated parameter values in Table 3, the value of  0 associated with model (4) is approximately  0 = 1.679, which is above the COVID-19 threshold quantity of one. It can be observed from Fig. 3 that as the value of the effective transmission rate, , increases the magnitudes of exposed, hospitalized and infectious classes increase. It means that this parameter is having a serious negative impact on the population and re-affirming the results obtained from the sensitivity analysis.
In Fig. 4, the effect of varying hospitalization rate for symptomatic infected students, ℎ 1 , on the dynamics of COVID-19 transmission and spread is demonstrated. It is observed from Figs. 4(a) and 4(c) that the number of hospitalized students and recovered students increases as ℎ 1 increases, and thereby leading to a fewer number of symptomatic infectious students in the population as Fig. 4(b) shows.
Similarly, illustration of the effect of hospitalization rate for selfquarantined students, ℎ 2 , on COVID-19 dynamics is graphically displayed in Fig. 5. Figs. 5(a) and 5(c), respectively, show that as ℎ 2 increases, the population size of symptomatically infected and hospitalized students increase. Moreover, in Figs. 5(b), 5(d)-5(f), the magnitude of quarantined, learned class, online facilities-unaccessed class and online facilities-accessed class decreases as ℎ 2 increases. Furthermore, Figs. 6(b)-6(c) reveal that the size of hospitalized class and infectious class, respectively, rise as the disease progression rate from the exposed class, , rises. Whereas from Fig. 6(a), it is observed that the size of the exposed class reduces as increases, with more students become learned as decreases (see Fig. 6(d)). Fig. 7 shows the effect of modification parameter for the relative infectiousness of hospitalized students, , on the population. It is shown that as increases, the numbers of hospitalized and symptomatically infected students increase with fewer number of students having access to online facilities (see Figs. 7(b)-7(d), respectively). In addition, Fig. 7(a) suggests that an inverse relationship exists between the exposed class and , as the number of exposed students reduces (increases) with an increase (decrease) in .
The time series plots of the effect of detection rate for exposed students, , on the COVID-19 dynamics is given in Fig. 8. It is observed that as increases, the magnitude of exposed, hospitalized and symptomatically infected student population decreases (as Figs. 8(a)-8(c), respectively, show), with more number of students get quarantined and recovered (see Figs. 8(d) and 8(e), respectively). Fig. 9 shows how maximal per capita recovery rate corresponding to the sufficient clinical resources and few hospitalized students ( 1 ) affects the dynamics of COVID-19 of governed by model (4). It is seen that the number of hospitalized students reduces (increases) as 1 increases (reduces) as Fig. 9(a) shows, leading to the recovery of more (less) infected students in the population as revealed in Fig. 9(b). This suggests that provision of a health care setting with increased well-equipped healthcare facilities will reduce the disease burden in the population.

Non-autonomous system
An iterative forward-backward sweep method (FBSM) based on the fourth-order Runge-Kutta method implemented in MATLAB is applied to the sixteen-dimensional optimality system, which is a two-point boundary problem, consisting of the state system (12) and the adjoint system (25) with the control characterizations (27) within the time interval of [0, 50] days in order to find the optimal control strategies required to minimize the spread of COVID-19 pandemic in a population at minimum cost. Because the optimality system has different time orientations, the solutions of equations of the non-autonomous system (12) with initial conditions and the control initial guess are sought forward in time. Whereas the equations of the adjoint system (25) with terminal conditions (26) are solved backward in time using the state system's current iteration solution. We refer the readers who are interested in the numerical procedure for the simulation of optimality system with different time orientations using FBSM to Lenhart and Workman [55].
The weight constants values and (where = 1, 2, 3, 4) appearing in the objective functional (13), in addition to the parameter values given in Table 3, are taken as 1 = 0.001, 2 = 0.015, 3 = 0.025, 4 = 0.015, 1 = 250, 2 = 150, 3 = 130 and 4 = 250. It is important to state that these weight constants values are theoretical as they are taken only to carry out the numerical illustrations of the control problem proposed and theoretically analysed in this work.
Strategy A: optimal use of preventive measure ( 1 ( )), quarantine (or selfisolation) control ( 2 ( )) and treatment control ( 3 ( )) In Fig. 10, the simulation of optimal control problem with the effort of personal protective measures 1 ( ) combined with quarantine or case tracing control 2 ( ) and treatment or hospitalization control 3 ( ) (i.e., 1 ≠ 0, 2 ≠ 0, 3 ≠ 0 and 4 = 0) is illustrated. Fig. 10(a) shows that the magnitude of symptomatic infected student subpopulation ( ) with control diminishes sharply at earlier time when compared with the case without controls. The size of hospitalized student subpopulation ( ) with controls slightly increased over the simulation time sub-interval [0, 7) days, and slightly reduced over the remaining time sub-interval (8,50] days as shown in Fig. 10(b). Further, it is seen in Fig. 10(c) that the magnitude of the disease prevalence in the absence of control strategy is slightly higher than the magnitude in the presence of control, which is more graphically perceptible within the time subinterval (7,50] days. The control profile presented in Fig. 10(d) shows that the optimal control 1 ( ) should be sustained at the maximum coverage (100%) in the first 37 days before dropping to the minimum level at the final time of control implementation, whereas optimal controls 2 ( ) and 3 ( ) are at the upper bound for 7 days (counting from the commencement of intervention strategy) before consistently maintained near the lower bound up till the final time.
Strategy B: optimal use of preventive control ( 1 ( )), quarantine (or selfisolation) control ( 2 ( )) and management control ( 4 ( )) The impact of synergy of the combination of personal protective measures 1 ( ) with quarantine (or self-isolation) control 2 ( ) and management control of hospitalized students 4 ( ) (i.e., 1 ≠ 0, 2 ≠ 0, 4 ≠ 0 and 3 = 0) on the dynamics of COVID-19 is investigated by simulating the optimal control problem. It is observed that the magnitude of symptomatic infected and hospitalized students subpopulations, and the disease prevalence with control are lower than the magnitude without control as shown in Figs. 11(a)-11(c). The control profile given by Fig. 11(d) shows that the optimal control 1 ( ) should be sustained at the upper bound for 18 days (counting from the commencement of control implementation day) before declining gradually to the lower bound of 0% at 25th day and constantly held at 0% throughout the remaining simulation period, control 2 ( ) is maintained at maximum level of 100% implementation in the first 6 days of control intervention before declining gradually to the minimum coverage of 0% at day 30, and consistently retained at this minimum coverage for the remaining 20 days of intervention period, meanwhile optimal control 4 ( ) is sustained at the maximum coverage of 100% throughout the entire 50 days of intervention.
Strategy C: optimal use of preventive measure ( 1 ( )), treatment control ( 3 ( )) and management control ( 4 ( )) To demonstrate the impact of combination of personal protective measure 1 ( ), treatment control ( 3 ( )) and management control ( 4 ( )) on the community spread of COVID-19, the optimal control problem is simulated. It is seen from Figs. 12(a)-12(c) that the total numbers of symptomatic infected students ( ), hospitalized students ( ) and the disease prevalence ( ) + ( ) in the case of control strategy are lower when compared with no control scenario. The control profile given by Fig. 12(d) shows that the optimal control 1 ( ) is at the maximum coverage of 100% in the time horizon [0, 18) days before declining sharply to the minimum coverage of 0% during the 23rd day of intervention and constantly retained at this minimum coverage throughout the remaining intervention period, control 3 ( ) is at the upper bound (100%) in the first 14 days before declining gradually to the lower bound (0%) during the 40th day and retained consistently at the lower bound in the last 10 days of the period of intervention, while control 4 ( ) is held at the maximum level of 100% throughout the entire 50 days of intervention.
Strategy D: optimal use of quarantine (or self-isolation) control ( 2 ( )), treatment control ( 3 ( )) and management control ( 4 ( )) Fig. 13 illustrates how the combination of quarantine (or selfisolation) control ( 2 ( )) with treatment control ( 3 ( )) and management control ( 4 ( )) halts the transmission dynamics of COVID-19 in the student population. In Figs. 13(a)-13(c), it is revealed that sizes of symptomatic infected student subpopulation ( ), hospitalized student subpopulation ( ) and the disease prevalence ( ) + ( ) with no control strategy scenario are higher when compared with the case of control intervention strategy implementation. The control profile given by Fig. 13(d) shows that the optimal controls 2 ( ) and 3 ( ) are at the maximum coverage of 100% in the first 15 and 13 days of intervention period before declining gradually to the minimum coverage of 0% during the 30 and 40 days of intervention, respectively, and constantly retained at this minimum coverage throughout the remaining intervention period, meanwhile control 4 ( ) is held at the maximum level of 100% implementation throughout the entire 50 days of intervention period.
Strategy E: optimal use of preventive control ( 1 ( )), quarantine (or selfisolation) control ( 2 ( )), treatment/hospitalization control ( 3 ( )) and management control ( 4 ( )) In Fig. 14, the graphical illustration of how the combination of personal protective measures 1 ( ) with quarantine (or self-isolation) control 2 ( ), treatment control 3 ( ) and management control 4 ( ) (i.e., 1 ≠ 0, 2 ≠ 0, 3 ≠ 0 and 4 ≠ 0) hampers the transmission dynamics of COVID-19 in the student population is considered. As shown in Figs. 13(a)-13(c), the sizes of symptomatic infected student subpopulation ( ), hospitalized student subpopulation ( ) and the disease prevalence ( ) + ( ) with no control strategy scenario are higher compared to when control strategy is in use. The control profile given by Fig. 14(d) shows that the optimal control 1 ( ) should be sustained maximally for 18 days (counting from the beginning of control implementation period) before declining sharply to the minimum coverage of 0% during the 23rd day of control intervention period and continuously retained at the minimum coverage throughout the remaining intervention period, 2 ( ) and 3 ( ) are both at the maximum coverage of 100% in the first 5 days of intervention period before decreasing gradually to the minimum coverage of 0% during the 30 and 36 days of intervention, respectively, and consistently retained at this minimum coverage throughout the remaining intervention period, whereas control 4 ( ) is held at the maximum level of 100% throughout the entire 50 days of intervention period.

Cost-effectiveness analysis
With the help of ideas from previous studies [33,40,41,43,44,[56][57][58][59], we explore the cost-effectiveness analysis (economic evaluation) of the five control combination strategies under investigation using three cost analysis methods, namely, infection averted ratio (IAR), average cost-effectiveness ratio (ACER) and incremental cost-effectiveness ratio (ICER), in this section. This is necessary as there is a need to determine the most effective and least costly intervention strategy among the various combinations considered based on the results arising from the simulations of the optimal control problem.

Infection averted ratio
The IAR is defined as IAR = Total number of infection averted Total number of recovered and protected .
The total number of infection averted, denoted as IA, is given in view of the objective functional (13) as where and are the total number of symptomatic and hospitalized individuals without control respectively, and and are the respective total number of symptomatic and hospitalized individuals with control strategy. According to this cost analysis approach, the strategy with the highest IAR value is the most cost-effective strategy [56,57]. Using (29), the IAR for each strategy is calculated. The numerical results for the implemented five strategies are demonstrated in Fig. 15 (see also Table 4). Strategy E involving the use of all the four optimal controls, namely, personal protective control ( 1 ( )), quarantine (or selfisolation) control ( 2 ( )), treatment control ( 3 ( )), and management control ( 4 ( )) gives the highest IAR. Thus, Strategy E is the most costeffective according to this cost analysis approach. This is followed by Strategy B which combines 1 ( ), 2 ( ) and 4 ( ) with 3 ( ) = 0. Strategy D consisting the combination of quarantine control, treatment control, and management control ( 2 ( ), 3 ( ) and 4 ( ) with 1 ( ) = 0) is the next Table 4 Intervention strategy, total infection averted, total cost, IAR and ACER.

Strategy
Total cost-effective strategy. This is followed by Strategy C, which involves the combination of personal protective control, treatment control, and management control ( 1 ( ), 3 ( ) and 4 ( ) with 2 ( ) = 0). The least cost-effective strategy is Strategy A involving the simultaneous use of personal protective control, quarantine control, and treatment control ( 1 ( ), 2 ( ) and 3 ( ) with 4 ( ) = 0). Strategy A being the least costeffective strategy is because the use of the strategy averts the lowest number of infections in the population as shown in Fig. 16 (see also Table 4). It is also important to note that the IAR values for Strategy B to Strategy E are relatively close to each other, with high discrepancy between each value and that of Strategy A as shown in Fig. 15 (see also Table 4). This indicates the importance of considering management control ( 4 ( )) in the control combination strategies to effectively reduce COVID-19 prevalence in the population. This further reaffirms the results demonstrated in Figs. 10 to 14.

Average cost-effectiveness ratio
The next cost analysis method considered in this paper is the ACER. This method deals with a single intervention strategy and weighing the intervention against its baseline option. The formula to calculate ACER is given as ACER = Total cost produced by the intervention Total number of infection averted .
In view of the objective functional (13), the total cost produced by an intervention is expressed mathematically as TC = ∫ 0 ( ) .
A strategy with the least ACER value is the most cost-effective according to this cost analysis approach [35,56,57]. Thus, the ACER is calculated for each of the five strategies using the formula (31). The numerical results obtained are presented by bar chart in Fig. 17 (see also Table 4). It is shown that Strategy B has the least ratio, and therefore is the most cost-effective strategy using ACER cost analysis approach. This is followed by Strategy D, Strategy E, then Strategy C. The least cost-effective strategy (a strategy with the highest ratio) is Strategy A.
To further affirm this result, we compute the ICER values for the various control combination strategies under consideration.

Incremental cost-effectiveness ratio
On the other hand, the ICER has to do with the comparison of the differences between the costs and health outcomes of the alternative intervention strategies (usually two or more) competing for the same resources. To compare competing intervention strategies (usually two or more) incrementally, one intervention strategy is compared with the next less-effective alternative strategy. Simply put, the ICER is stated as ICER = ▵ in costs expended on control strategies implementation ▵ in total number of infections averted by the competing strategies .
In Table 5, the total numbers of infections averted by strategies A, B, C, D and E are arranged from the least to the highest alongside the respective associated costs expended on the implementation of the five strategies. Thus, using (33), we obtain the numerical results in Table 5 following the computation of ICER as follows:   Comparing Strategy A and Strategy C in Table 5, it is seen that ICER(C) is less than ICER(A). This indicates that Strategy C strongly dominates Strategy A. Thus, Strategy C has greater effectiveness at cheaper cost compared to Strategy A. Therefore, Strategy A is left out from further consideration. Next, Strategy C is further compared with Strategy D. Hence, using (33), the summary of ICER is given in Table 6. A look at Table 6 reveals that Strategy D (from the negative ICER value) strongly dominates Strategy C as ICER(C) is greater than ICER(D). This means that strategy C is less effective and more expensive to implement than Strategy D. Therefore, it is best to remove Strategy C from the set of intervention strategies to implement in order to preserve the limited resources. Consequently, Strategy C is excluded and Strategy D is further compared with Strategy B. Thus, using (33), the summary of ICER is given in  Strategy E. Consequently, the most cost-effective strategy according to ICER cost analysis approach is Strategy D. Arising from the foregoing cost analysis, Strategy E (combination of optimal preventive control, optimal quarantine control, optimal treatment control and optimal management control) is the most cost-effective strategy according to IAR method, Strategy B (combination of optimal preventive control, optimal quarantine control and optimal management control only) is the most cost-effective strategy in view of ACER approach, while Strategy D which combines the efforts of optimal personal protective measures with optimal quarantine (or self-isolation) control of exposed individuals and optimal management control of hospitalized individuals (active cases) is confirmed to be the most cost-effective intervention strategy capable of diminishing the COVID-19 burden optimally in the sense of ICER cost analysis technique. It is worthy of note that Strategy D produces the least total cost required to achieve the optimal control of COVID-19 in the population as given in Table 5. This is clearly presented in Fig. 18, which also reaffirms the results obtained from the IAR, ACER and ICER methods that Strategy D is the most cost-effective strategy.

Conclusion
When the number of COVID-19 cases started to rise, closure of educational institutions, including higher institutions, was one of the measures adopted by the governments to contain the community spread of the global epidemic. In developing (low-and middle-income) countries, almost all educational institutions were physically closed for exceptionally long periods, ranging from several months to years. In particular, the long-term closure of higher institutions induced worse learning outcomes, physical and mental issues, and crises among students. Alternatively, different countries adopted online education via different media such as television and internet. However, this alternative education system could add minimal value to students' learning in many developing countries due to the limited available resources [60]. Thus, this study presents a mathematical analysis of transmission dynamics of COVID-19 in a population of student of higher institution of learning via online with limited available online and healthcare facilities to gain further insights into the impact of the pandemic on higher education and to explore the optimal strategy required for minimizing the number of infected students at minimum cost. The first part of this work formulates an autonomous model with a system of ordinary differential equations by subdividing the student population into nine mutually exclusive compartments of susceptible, exposed, quarantined, symptomatic infectious, online facilities-accessed, online facilities-unaccessed, learned, hospitalized and recovered. By analysing the fundamental properties of solutions of the model, it is shown that the model is mathematically well-posed. The disease-free equilibrium of the model is obtained, and proved to be locally asymptotically stable whenever the basic reproduction number associated with the model,  0 , is below unity. Sensitivity analysis is performed to explore the model parameters that drive the transmission dynamics of the disease in the student population. It is revealed that the effective transmission rate ( ), modification parameter for reduced infectiousness of hospitalized students ( ), rate of infectious of exposed students ( ), maximum per capita recovery rate corresponding to the sufficient healthcare facilities and few hospitalized students ( 1 ), detection rate for exposed students ( ) and hospitalization rate for symptomatic infectious students (ℎ 1 ) are the most sensitive parameters.
Moreover, the autonomous model is extended to a non-autonomous system of ordinary differential equations in view of the epidemiological insights arising from the sensitivity analysis in the second part of the work. The new optimal control framework captures four time-varying control variables, namely personal protective measure, quarantine (or self-isolation), treatment, and management measure to assess the best optimal strategy needed to minimize the spread of COVID-19 in a population of student of tertiary institutions at cheapest cost. The non-autonomous system is analysed by employing optimal control theory with the aid of Pontryagin's maximum principle. The effects of five strategies arising from different combinations of the four control variables are investigated on the dynamics of COVID-19 spread in the population. It is found that, if strict adherence is practised, each of the combination strategies has positive effect in reducing the number of COVID-19 infections, however any combination strategy that includes management control of hospitalized students is highly effective in stemming the transmission and spread of the disease in the student population. This result showcases the relative importance of enhancing the management measures (such as prompt provision of supplemental oxygen or mechanical ventilation) for hospitalized students in order to ensure their prompt recovery and prevention of deaths due to complications in the course of fighting against COVID-19. Furthermore, the results arising from cost-effectiveness analysis suggest that a control intervention strategy that combines optimal use of quarantine (or self-isolation) control for exposed students, treatment control for symptomatic infectious students and management control for hospitalized students is the most effective strategy with minimal cost that could be suggested for the optimal control of COVID-19 spread among the five different combinations of optimal control variables analysed in this paper.
Further work can be conducted by modifying the nonlinear mathematical models of COVID-19 dynamics developed and analysed in this paper by introducing a random mass vaccination programme to enhance the four control variables examined herein.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
All data generated or analysed during this study are included in this article.