Application of optimal control to the dynamics of COVID-19 disease in South Africa.

SARS-CoV-2 (COVID-19) belongs to the beta-coronavirus family, which include: the severe acute respiratory syndrome coronavirus (SARS-CoV) and the Middle East respiratory syndrome coronavirus (MERS-CoV). Since its outbreak in South Africa in March 2020, it has lead to high mortality and thousands of people contracting the virus. Mathematical analysis of a model without controls was done and the basic reproduction number ( R 0 ) of the COVID-19 for the South African pandemic determined. Permissible controls were introduced and an optimal control problem using the Pontraygain Maximum Principle is formulated. Numerical findings suggest that joint implementation of effective mask usage, physical distancing and active screening and testing, are effective measures to curtail the spread of the disease in the human population. The results obtained in this paper are of public health importance in the control and management of the spread for the novel coronavirus, SARS-CoV-2, in South Africa.

Introduction

The recent pandemic, which is commonly known as novel coronavirus (COVID-19), is a lower respiratory infection which originated from Wuhan China in 2019 and it has spread all over the world [1]. The first case for COVID-19 occurred in early December, 2019 and by January 29, 2020 more than 7000 infected cases had already been reported in China [2]. The virus spreads through coughing, sneezing and most people who fall sick with COVID-19 usually experience mild to moderate symptoms and recover without special treatment [3]. The burden of COVID-19 has been massive worldwide, and governments are finding ways to control the pandemic and protect the health systems of the countries while safeguarding their economies. There are no antivirus drugs yet to fight the virus and scientists around the world are working on potential treatments and vaccines for the pandemic [4]. The only available mode of minimising the disease burden has been to control its transmission in the population using non-pharmaceutical interventions. Mathematical models have played a major role in increasing the understanding of the underlying biological and psycho-social processes that control the spread of infectious diseases and how the disease can be managed [5]. Various control strategies have to be deployed to manage the spread of COVID-19. It is thus important to consider standard mathematical approaches to determine the optimality of the control strategies. Optimal control is a standard method for solving dynamic optimisation problems, where these problems are defined in continuous time [6]. To formulate an optimal control problem, a mathematical model of the system is required to control specifications of all the boundary conditions on states and constraints to be satisfied by forms and controls [7]. The first case of the COVID-19 disease in South Africa was reported on March 5, 2020 as an imported case from Italy [8]. Due to the increase in new number of cases, the South African government introduced strategies to control the spread of the pandemic such as social distancing to limit close contact among individuals, contact tracing, testing and screening [9]. On March 26, 2020, the South African government enforced a 21-day national lockdown to restrict contact of individuals and to delay the spread of the virus. As of August 3, 2020, the number of infected cases in South Africa for all the provinces was 516 862 out of which 358 037 recovered [10]. When compared to July 3, 2020, the number of cases had increased by 192% and recoveries by 315% in a period of four weeks.

Since the outbreak of the pandemic, researchers have used mathematical and statistical models to predict and model the extent to which the disease can be contained, or controlled. Nyabadza et al. [11] considered the impact of social distancing on the spread of COVID-19 in South Africa. Their results showed that an increase in social distancing significantly decreases the disease spread. Masuku et al. [12] used a system of differential equations to quantify the early COVID-19 outbreak transmission in South Africa and exploring vaccine efficacy scenarios. Their findings suggested that the COVID-19 outbreak in South Africa had a basic reproductive number of 2.95 and a highly efficacious vaccine would have been required to contain COVID-19 in South Africa. Chandra et al. [13] looked at mathematical models with a social distancing parameter for the early estimation of COVID-19 spread in India. Their results showed that mathematical modelling and simulations on social distancing play an important role in spread estimation. The effect of social distancing was discussed with different social distancing rates and they found that social distancing can reduce spread of COVID-19. A normalised Susceptible, Exposed, Infected and Recoveries mathematical model with optimal control techniques was considered in [14] with the aim of establishing vaccination and treatment plans that reduce the spread of an infectious disease. Their findings revealed that after stopping vaccination and/or treatment, the percentage of the infected people increased. The study in [15] proposed a generalised model of COVID-19 to study the behaviour of the disease transmission under five different control strategies. The results from the numerical simulations reflected that quarantining and better medical treatment of infected individuals reduced the critically infected cases, and the transmission risk. Emerging evidence and predictions on COVID-19 allude to the fact that the infection may become endemic and may never go away [1]. Any intervention strategies would need to adjust to this new reality. Unlike most of the current models which ignore the vital dynamics, models need to accordingly be adjusted to capture the nature of the pandemic. The current model incorporates the vital dynamics to capture the dynamics of COVID-19 infection using the South African setting in addition to optimizing the control strategies. The objective of this study is to develop dynamics COVID-19 model for the South African setting incorporating current infection prevention efforts that include screening and testing, the use of masks and physical distancing. Optimal control theory is used to quantify the effectiveness of these strategies. The paper is arranged as follows: in the next section, the relevant model is formulated. Section 3 presents the model properties and analysis without controls. In Section 4, the model with optimal control is presented, analysis is done using Pontryagin’s Maximum Principle. Numerical results are presented in section 5 and the last section concludes the paper.

Methods

Exploratory data analysis of the COVID-19 in South Africa

The Exploratory Data Analysis (EDA) helps to reveal some hidden features in data sets, a good practice to understand the data first and gather as many insights from it [16]. In this research, EDA is used to indicate hidden features in COVID-19 data of South Africa, and data to be modelled is collected by the Data Repository and Dashboard for South Africa [10] and official government press release [17], from March 03, 2020, when the first case for South Africa was identified. This data set has daily and cumulative information on the number of COVID-19 infected patients, deaths, and recoveries. On March 26, 2020, the government implemented a 21-day national lockdown due to increased daily confirmed cases. The restrictions like social distancing, wearing of masks were implemented to control the spread of the disease.

Table 1 shows the basic features of the COVID-19 data, where represents the first quartile from box plots and represents the third quartile from the box plots. This study presents COVID-19 data between March 26, 2020, and July 19, 2020. On March 26, the total number of cases reported infected 927, and out of these, 218 were new cases. Between March 26 and July 19, the minimum number of cumulative infected patients is 927, and the maximum is 364 328; the mean is 68 278.77 with the standard deviation of 94 256.84. The coefficient of variation (CV) is equal to 1.38, which indicates a relatively high variation, thus causing high variability in the number of COVID-19 confirmed cases in South Africa. For the recoveries population, the minimum number is 0, and the maximum is 191 059; the mean is 34040.64 with the standard deviation of 47 407.55. The CV is 1.39 and also indicates a relatively high variation. For the reported deaths, the minimum number is 1, and the maximum is 5033; the mean is 1152.82 with a standard deviation of 1422.69. The CV is 1.23 and also indicates a relatively high variation. The skewness for all three population is greater than 1, this means that the data sets are highly skewed.

Table 1

Basic features of South African COVID-19 Data.

Statistics	Infected	Recoveries	Deaths
Minimum	927.00	0.00	1.00
Maximum	364328.00	191059.00	5033.00
Q_1	4290.50	1473.00	82.50
Median	21343.00	10104.00	407.00
Q_3	94991.50	50967.00	1903.50
Mean	68278.77	34040.64	1152.82
Range	363401.00	191059.00	5032.00
Standard Deviation	94256.84	47407.55	1422.69
Skewness	1.63	1.67	1.25
Kurtosis	1.69	2.05	0.40

Fig. 1 show the time plots of daily infected individuals, daily recoveries and daily deaths cases. EDA is primarily used to see what data can reveal beyond the formal modeling or hypothesis testing task and better understand dataset variables and the relationships between them. It can also help determine if the statistical techniques you are considering for data analysis are appropriate [18].

Fig. 1

Plots representing EDA of the COVID-19 in South Africa for: (a) Daily infected cases. (b) Daily recovered class. (c) Daily deaths cases.

Mathematical model without control

In this section, a brief model description of the transmission dynamics of COVID-19 without control measures is given. Considering a non-constant population, the model is described using a set of five compartments of the susceptibles, who are recruited at a constant rate These individuals can get infected [19], and move into the exposed compartment Individuals in this compartment do not show symptoms and are assumed not to transmit the infection. After a period of 5–6 between 5 to 14 days, an exposed individual progresses to become symptomatic and infectious [20]. The exposed individuals are assumed to become symptomatic at a rate Of those that become symptomatic, a proportion is assumed to be undiagnosed due to system capacity failure and the non-development of severe disease and move to the compartment while the remainder is diagnosed and move to the compartment The individuals in are assumed to undergo active screening and testing and when confirmed positive, they move to the class at a rate The diagnosed and undiagnosed infectious individuals recover at rates and respectively and move into the recovered compartment Individuals in all compartments are assumed to have a natural mortality rate The diagnosed and undiagnosed individuals have additional disease induced death rates and respectively. The force of infection is defined bywhere is the infection rate, and is the modification parameter used to compare the infectiousness of individuals in relative to those Here we assume that the class infectivity is less that that for the since individuals who know their status are more likely to limit contact than those who do not know. The total human population for the COVID-19 transmission dynamics at any time is given byThe flow diagram of individuals between compartments is shown in Fig. 2 .

Fig. 2

Model flow chart depicting the dynamics of COVID-19 transmissions.

Considering the model description, assumptions, model parameters and Fig. 2 we have the following system of ordinary differential equations:subject to the initial conditionsSince the recovered humans do not contribute to COVID-19 disease transmission, class can be eliminated from the system using the relationshipand still adequately represent the dynamics of system (1) by

Mathematical analysis

Boundedness of solutions

The region in was proposed for model Eq.  (3) to be well-posed. These solutions are shown to be uniformly bounded. Thus, the following theorem is established:

The solutions of model system(3)are contained in the regionwhich is given byfor the initial conditions(2)and the regionis attracting.

The change in total human population of Eq.  (3) is given byby the positivity of The solution of (4) yieldswhere is a constant representing the initial condition. Thus, This result implies that the model describing the human population is epidemiologically well-posed, and the solutions remain in the region  □

Positivity of solutions

Here, we show that the COVID-19 disease model (3) solutions remain positive in the positive orthant given any non negative initial conditions.

Given the initial conditions, the solutionsof the system(3)are all positive for all

From the first equation of the system (3), we have thatwhose solutions givesThus, will always be positive for all Also, from the second equation of the system (3), we have thatwhose solution yieldsThus, will always be positive for all Similarly, the solutions of the remainder of the equations in (3) is straight forward and are all positive; that is, and over the modelling time. □

Equilibrium points

The disease free equilibrium

Assuming that at the disease free equilibrium (DFE) point there are no infectious individuals. Setting the right hand side to zero, we obtain the DFE denoted by

The reproduction number

For simplicity, we set in (3). The basic reproduction number denoted by () is the average number of secondary cases generated by an infectious individual who introduced in a wholly susceptible population. We calculate the COVID-19 for model system (3) using the next generation matrix method in [21]. We consider three compartments and which contribute to new infections or the secondary cases and construct the column matrices;where, represents the rate of generation of new infections in compartments and and represents the transition in and out of the compartments. We define matrices and to the Jacobian matrices for and respectively evaluated at the DFE. The basic reproduction number is the spectral radius of the product of and i.e, where in the inverse matrix of the matrix V. Computing the matrices and we obtain;The spectral radius is thus the maximum eigenvalue of the second generation matrix and is given bywhereThe expression represents the number of new COVID-19 infections generated by infectious undiagnosed individuals, It contains the product of the transmission rate the fraction of exposed individuals moving to undiagnosed class and is the average period spend in the class. represents the number of new COVID-19 infection cases generated by diagnosed infectious individuals, is the product of infectious rate due to diagnosed COVID-19 individuals, fraction of exposed individuals moving to diagnosed infectious stage and average period spent in the compartment.

It is easy to show that our model satisfy the hypotheses of the Theorem 2 in [21] and hence, the stability of the DFE is ensured and summarised as follows:

The DFE is locally asymptotically stable wheneverand unstable otherwise.

The endemic equilibrium

After some algebraic manipulations, the endemic equilibrium denoted as (EE) is given bywhere withWe have the following result on the existence of

Note that if we set reduces to

The EE exist if and only if

Global stability analysis of the equilibria

In this section, we show that the disease-free equilibrium is globally asymptotically stable if We also prove that the endemic equilibrium state is globally asymptotically stable if

The disease-free equilibrium,of the model system(3)is globally asymptotically stable if

Consider the following Lyapunov function defined bywhere and are the positive constants we seek to determine. Finding the time derivative of the function and upon substitution of Eq.  (3) at the DFE we have that

Eliminating and the linear terms, we solve for and to getSubstituting for and we obtain

We see that if and if Therefore, the largest compact invariant set in such that   when is the singleton Thus, by LaSalle’s Invariance Principle [22], [23], the steady state is globally asymptotically stable whenever  □

Ifthe endemic equilibrium of the model Eq.  (3) , is globally asymptotically stable in .

Given that the endemic equilibrium exists if and only if Consider Lyapunov functionwhere and are constants to be determined. Substituting the expressions for the derivatives of the state variables yields,

At the steady states we haveGiven that and setting the time derivative of the Lyapunov function isThe derivative of thus reduces towhereAfter some algebraic manipulation, we obtain and so thatGiven that the arithmetical mean is always greater or equal to the geometric mean, the expressionsSo for all and if and only if and The largest invariant set where is the singleton of By Lasalle’s invariance principal [23], is globally asymptotically stable with respect to the invariant set  □

Optimal control problem

Mathematical model with optimal control

We introduce optimal control measures to the system (3). We formulate an optimal control problem with the following control variables: the prevention effort of COVID-19 disease such as mask usage, physical distancing, and , active screening and testing. The new force of infection becomesand the control variables to model Eq.  (3) gives rise tosubject to the initial conditionsThe objective is to reduce the number of undiagnosed humans with the COVID-19 disease. Hence, we formulate a minimization problems with the following objective function where is the initial time, is the terminal time, is the weight associated with reducing the undiagnosed infected human class , are the associated cost weight with the control variables and respectively. We define the Hamiltonian function by applying the Pontryagin’s Maximum Principle [24] as follows:where are the adjoint variables. We now give the theorem for the existence and uniqueness of the optimal control problem as in [6] as follows;

There exist an optimal controlsuch thatsubject to the optimal control model system(5)with initial conditions(1).

Optimality of the system

Given that and are the solutions of the optimal control model system (5) and (6) associated with the optimal control variables Then, there exist an adjoint system which satisfies where

with transversality boundary conditions defined asgiven bywhere the permissible control functionsand, are obtained by settingThusFrom Eq.  (11) , we solve for and to obtain their respective permissible control, solutions, which are With the upper and lower constraints on the admissible controls and we now have that the optimal controls can be characterised as

Model fit with controls and parameter estimation

The method of least square curve fitting is used to fit system (3) to the data. The method is useful in obtaining the parameters values that are then used in the numerical simulations. Fig. 3 shows that model (3) fits well to the COVID-19 new daily cases of South Africa as from the 26 of March up to the 19 of July 2020, with data obtained from [10]. The curve fitting process generates the following parameters with some of the parameters having been obtained from the cited research.

Fig. 3

Model fit to data for COVID-19 early growth of new cases in South Africa for the first wave from 26 March to 19 July 2020. Optimal parameter values:

Sensitivity analysis

Sensitivity analysis of model (3) is performed using the methodology described in [27]. In this paper, the objective of sensitivity analysis is to identify critical parameters driving the COVID-19 pandemic in South Africa.

Fig. 4(a) shows the effects of parameters and on undiagnosed infectious population. is largely negatively correlated to . The increase in the rate of active testing is associated with the decrease of undiagnosed infectious class. is mostly negative correlated to . The increase in the rate at which exposed become infectious is associated with the increase in the undiagnosed class. is mostly positive correlated to . The increase in proportion of exposed individuals which becomes undiagnosed is associated with the increase in the undiagnosed class. is mostly positive correlated to . The increase in the infection rate is associated with the increase in the undiagnosed class. is negatively correlated to . The increase in the mask usage and physical distancing is associated with the decrease in the undiagnosed class. is negatively correlated to . The increase in the active testing and screening is associated with the decrease in the undiagnosed class. Fig. 4(b) shows the effects of parameters and on diagnosed infectious population. is largely negatively correlated to . The increase in the rate of active testing is associated with the increase of diagnosed infectious class. is mostly positive correlated to . The increase in the rate at which exposed become infectious is associated with the increase in the diagnosed class. is mostly positive correlated to . The increase in proportion of exposed individuals which becomes undiagnosed is associated with the decrease in the diagnosed class. is mostly positive correlated to . The increase in the infection rate is associated with the increase in the diagnosed class. is negatively correlated to . The increase in the mask usage and physical distancing is associated with the decrease in the diagnosed class. is negatively correlated to . The increase in the active testing and screening is associated with the decrease in the diagnosed class.

Fig. 4

PRCC plots showing the effects of the dynamics of parameter values on COVID-19 daily new cases in South Africa for (a) parameters and on (b) parameters and on

Effects of parameters on the

The effect of parameter on were analysed using a contour plot. We chose two significant parameters and and give the contour plot as a function of Fig. 5 shows that an increase in the parameter results in a decrease in the numerical value of This implies that increased screening has the propensity to decrease the spread of COVID-19. It is thus important to increase the rate of screening and testing in order to control the pandemic. On the other hand, the parameter does not significantly impact

Fig. 5

Contour plot of as a function of versus

Numerical results

The optimal control strategy is found by the iterative method of the fourth order Runge-Kutta iterative scheme. The state equations are initially solved by the forward Runge-Kutta method of the fourth order. Then, by using the backward Runge-Kutta method of the fourth order (backward sweep), we solved the co-state equations with the transversality conditions. The permissible controls are updated by using a convex combination of the previous controls and the value from the characterisations of and This process is re-iterated and the iteration is ended if the current state, the adjoint and the control values converge sufficiently following [6]. Also, the numerical simulations were performed based on the control variables, to investigate the effects of introducing time dependent controls. The combination of the following three different scenarios is considered, enumerated as follows: (1) the use of only one permissible control (2) the use of only one permissible control and (3) the use of a combination of permissible controls and Using the model parameter values in Table 2 , the initial conditions as defined in Section 4.3 and the associated weights from Tables 3 to simulate the model with controls and obtain the following results.

Table 2

Parameters of (COVID-19) model for South Africa and their values.

Symbol	Parameter description	Value [Range] day^−1	Source
Λ	Recruitment rate	11244 [10000 - 25000]	Estimated in[11]
μ	Mortality rate	0.0091 [0.0050 - 0.010]	[25]
σ	Rate of active testing	0.2070 [0.4000 - 1.000]	Data fit
δ_1	Diagnosed induced death rate	0.1500 [0.0000 - 0.100]	[26]
δ_2	Un-diagnosed induced death rate	0.0089 [0.0000 - 0.100]	Data fit
η	Modification parameter	0.5663 [0.0000 - 1.000]	Data fit
γ	Rate at which exposed become infectious	0.5732 [0.2000 - 0.600]	Data fit
ρ	Proportion of E which becomes undiagnosed	0.9799 [0.0000 - 1.000]	By definition
ω_1	Recovery rate of diagnosed infected humans	0.1800	Data fit
ω_2	Recovery rate of undiagnosed infected humans	0.0206 [0.0000 - 0.100]	[26]
β	Infection rate	1.3295 [1.2100 - 1.500]	Data fit
u_1	Mask usage and physical distancing	0.7557 [0.0000 - 1.000]	Data fit
u_2	Active screening and testing	0.3728 [0.0000 - 1.000]	Data fit

Table 3

Costs associated with optimal control variables used for simulations.

Symbol	Parameter description	Base line Values (day^−1)	Source
α	Weight associated with reducing the undiagnosed infected human	1	Estimated
α_1	Weight on the cost of the prevention	0.9	Estimated
α_2	Weight on the cost of the active screening and testing	0.5	Estimated

Scenario (1)

This scenario shows the dynamics of COVID-19 using the control variable which consist of mask usage and physical distancing to optimise the objective function while setting the control variable to zero. We observed that in Fig. 6 (a), that maintaining a strict masks usage and physical distancing have significant influence on the peaks of undiagnosed cases of COVID-19. In the absence of active screening and testing, masks and social distancing has the propensity to minimise the spread of COVID-19 disease in South Africa. Fig. 6(b), shows the control profile for this particular first scenario, where is given the maximum of approximately 20 days to drop to its lower bound, translating to a maximum effectiveness of masks usage and social distancing.

Fig. 6

The dynamics on COVID-19 using control with for (a) undiagnosed humans (b) the control profile for

Scenario (2)

This scenario presents the dynamics of COVID-19 using the control variable for screening and testing to optimise the objective function while is set to zero. We observed from Fig. 7 (a), screening and testing reduces the peak of undiagnoseed infectives. Fig. 7(b) shows the control profile for , where should be allowed to run at its highest level during the entire duration of the pandemic. This means that the scaling of screening and testing should always be done at capacity to control COVID-19.

Fig. 7

The dynamics on COVID-19 using control , with for (a) undiagnosed humans (b) the control profile for

Scenario (3)

This scenario presents both control variables and being run concurrently. For these scenarios, we observed in Fig. 8 (a) that the control strategies of wearing masks, physical distancing and as well screening and testing resulted in a significant decline in the number of undiagnosed COVID-19 cases. Thus, this shows that to reduce the spread of COVID-19, these controls should be used concurrently. Fig. 8(b) shows the control profile for both and .

Fig. 8

The dynamics on COVID-19 using controls and for (a) undiagnosed humans (b) the control profile for and

Conclusions

In this paper, a continuous time model was used to study the effects of the dynamics of COVID-19 pandemic on the South Africa population. The model takes into account the susceptible, exposed, undetected, detected and recovered humans. Mathematical analysis for the model without controls was carried out and showed that the steady steady states are locally asymptotically stable if the basic COVID-19 reproduction number is greater than unity and otherwise unstable. Further, the sensitivity analysis of the model equations without control was performed and numerical simulations based on the highly sensitive model parameters on as shown in Fig. 4. On the other hand, a time dependent optimal control problem was formulated by introducing control variables. The model with control system (5) was solved analytical and numerically to predict the possible and necessary controls to contain the spread of the disease. The results on the numerical simulations suggest that the best strategy is to wear masks and physical distancing followed by control of testing and screening to reduce COVID-19 infection and these strategies will be more effective in controlling the rate at which the uninfected humans contract the disease.

While many research papers on COVID-19 have been published to date, the work presented here focuses on the dynamics of the infection in the presence to time dependent optimal controls. Given that new scientific evidence on the virus continues to emerge, the model presented here is certainly neither the ultimate dossier not without culpability. The results should then be treated with circumspection. Some of the major challenges, other than the continuous evolution of scientific evidence, are the none inclusion of some of the dynamics of COVID-19 such as the role played by contaminated surfaces, asymptomatic infectives and the age dependant disease dynamics. Aspects of stochasticity in the infection rates, the role of social-economic factors in COVID-19 infection spread and the evolution of government policies can be incorporated in the model to improve its usefulness. Also, in [20], it is documented that exposed individuals in compartment can be contagious. Therefore, the assumption that individuals in compartment can not transmit COVID-19 can be relaxed. This means that transmission from individuals in compartment can occur before the onset of symptoms. Such an adjustment requires the inclusion of compartment in the force of infection. Despite these fallibilities, the model still presents some interesting results on the introspections of COVID-19.

Data Availability

The data are sourced from public and official government press release [17].

Funding Statement

No funding or grant were utilised for this research work.

Declaration of Competing Interest

No conflict of interest.