Control and elimination in an SEIR model for the disease dynamics of COVID-19 with vaccination

1.   Introduction

The COVID-19 pandemic has imposed a massive global challenge in public health for many countries, following its discovery early in December 2019 in Wuhan, China. Since then, the virus has rapidly spread to other provinces in China and the rest of the world. The first case of the COVID-19 disease in South Africa was reported in March 2020, ^[29]. Due to the increase in COVID-19 incidence, 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. On 26 March 2020, in order to curb the spread of the corona virus, the South African government enforced a national lockdown to restrict contact of individuals ^[7,12,17].

Prior to COVID-19, the world has already seen some serious outbreaks of disease from corona viruses, see, for instance, ^[23]. Mathematical models of COVID-19 disease dynamics, e.g., ^[6,14,18,36], have been far more complex than, for instance, the 4-compartment model in ^[23]. Recent published work, by way of mathematical models have shown that short term interventions such as physical distancing, hand washing and mask wearing are useful in reducing the COVID-19 case incidence ^{[14,18,22,30]}. Such non-pharmaceutical intervention were implemented world-wide during the heavy impact of COVID-19, see, for instance, in ^[2,3,25].

However, long term interventions including vaccination are highly desirable, given the socio-economic costs of lockdowns and physical distancing ^[11]. Since, vaccine was introduced, there has been much research on its effect in the control of COVID-19 in the population. Mathematical modeling studies in this regard can be found in ^{[1,4,5,6,18,26,36]}. The current paper is another contribution. We here present a compartmental model in ordinary differential equations (ODEs) for COVID-19 with a generalized contact rate and vaccination. Our paper focuses on the intervention that took place in South Africa for the reduction of the COVID-19 cases, which includes education, lockdown, law enforcement, etc., and vaccination. A detailed analysis of the model is included, and we also include an optimal control problem.

Regarding the form of the model, the closest work to ours seems to be ^[5]. The model in ^[5] has a nonlinear recovery rate. In the current paper the contact rate is more general. Our more general contact rate is inspired by the work of ^[10,13,24] and others. Our results have a slightly different focus from that in ^[5], and we also study an optimal control problem. Into the model we introduce two control functions: a control on the rate of vaccination of susceptibles and another control on the contact rate. The latter control amounts to introduction and enforcement of measures such as lockdown, social distancing, mask-wearing, sanitizing of hands and of objects and education on these preventive measures. Optimal control theory determines the best scheduling of the two strategies for the most effective results in combating COVID-19, balanced against the cost of these interventions. Our paper is one of only a few papers on optimal control of COVID-19 in South Africa ^[8,18,25]. These papers consider control on only non-pharmaceutical interventions. The current paper is the only one that includes vaccination in the model.

The paper is set up as follows. In Section 2 we present the model, and we prove existence of the global positive solutions. Section 3 focuses on the disease-free equilibrium point and its stability, and Section 4 deals with the endemic equilibrium point. The optimal control problem is presented in Section 5. In Section 6 we illustrate our theoretical results by way of numerical simulations with reference to South Africa. The paper is concluded in Section 7.

2.   The model

2.1.   Model description

We consider a population of size $N(t)$ at time $t$, which is subdivided into four classes. These classes are the susceptibles $S(t)$, the infectious but asymptomatic cases $E(t)$, the symptomatic infectious cases $I(t)$, and recovered individuals $R(t)$. Therefore, we have the following equation:

Members of the symptomatic COVID-19 class $I$ are indeed infectious. However, they are assumed to be so effectively isolated that eventually they do not contribute significantly to the transmission of the disease. The model is displayed in the flow diagram below, as Figure 1:

We use the following notation for the (mostly constant) parameters:

The contact rate is a function $\theta(I)$, of the infection class $I$. The vaccination rate is a non-negative constant. The rest of the parameters listed above are assumed to be positive constants. The function $\theta(I)$ is assumed to be positive-valued, and to have a continuous derivative $\theta'(I)\leq 0$. We shall write:

$\theta_0: = \theta(0)$ and $\theta_1: = \theta'(0)$.

The motivation for such a contact rate is that it reflects the response of the susceptible population with respect to lower or higher values of $I$. The idea is that as the susceptibles become aware of the magnitude of the higher numbers of symptomatic infections, they become more cautious in terms of their contact with people, especially those suspected of being infected. This leads to a reduction in the contact rate. The model is described by the system of ODEs (2.1).

The initial conditions are:

We write $N(t)$ as the human population size at time $t$, i.e.,

$S(t)+E(t)+I(t)+R(t) = N(t)$.

We shall adopt the following notation:

$\mu_0 = \mu+z, \mu_1 = \mu+\delta_1+\alpha+\omega_1$, $\mu_2 = \mu+\delta_2+\omega_2$ and $(A_0+A_1)/\mu = P,$

and we introduce the set $\Gamma$,

The following theorem is proved by contradiction, following an argument that has been popularly used, for instance, in ^[19,20].

Theorem 2. Given any $y\in \Gamma$, then there exists a unique solution $X(t)$ of the system (2.1) with $X(0) = y$ and with $X(s) \in \Gamma$ for each $s\in [0, \infty)$.

Proof. Consider any point $y\in \Gamma$. Then, there exists a local solution $X(t)$ in $\Gamma$ which is a $C^1$ function and with initial value $X(0) = y$. Suppose that $t_1$ is the time of exit from $\Gamma$. Since $X(t)$ is continuous, and $\Gamma$ is an open set, it follows that $t_1 > 0$. We shall show by contradiction that $t_1 = \infty$. Thus, suppose that $t_1$ is finite.

Let us define the following functions $L_0(X(t))$ for $0\le t \le t_1$ (and note that for notational simplicity, we suppress the function's dependence on $t$).

Each term on the right-hand side of equation (2.3) is a non-negative valued function of time $t$ for $0\le t < t_1$. We note that for each of the first four terms on the right-hand side of equation (2.3) above, if any of the $X_i(t)$ tends to $0$ as $t\to t_1$, then $L_0(X(t)) \to \infty$. We also note that for a positive constant $c$ and a variable $z$ with $0 < z < c$,

Therefore, if the explosion time is finite, then

Now, we calculate the derivative (and we suppress the time-argument $(t)$):

Note that $\theta E < \theta_0P$. After cancellation of some terms which are negative, we obtain an inequality:

Thus, $\frac{dL_0}{dt}\leq B_0$, with $B_0$ being the constant

Therefore, over the bounded interval $[0, t_1)$, $L_0(t)\leq B_0t_1$, and so $L_0$ is bounded over the interval $[0, t_1)$. This is a contradiction to the statement of equation (2.4). Therefore, we can conclude that $X(t)$ never exits the set $\Gamma$. This concludes the proof.

Following the methodology as in ^[31], we obtain the basic reproduction number as follows:

We note that the variable $R$ does not appear in the first three equations of the system (2.1). Therefore, we can consider the reduced system (2.7) below, derived from the system (2.1) by ignoring the $R$-equation, i.e., the fourth equation.

3.   Stability of the disease-free equilibrium

In this section we focus on the case $A_1 = 0$. In this case the model (2.7) permits a disease-free equilibrium

Theorem 3. If $A_1 = 0$ and $\mathcal{R}_0 < 1$, then the disease-free equilibrium of the model (2.7) is globally asymptotically stable.

Proof. Suppose that $\mathcal{R}_0 < 1$. Then, $\frac{\theta A_0}{\mu_0}-\mu_1 < 0.$ There exists $c > 0$ such that

Fix such a number $c$ and define a function $L_1(X(t))$ as:

The time derivative of $L_1$ is computed as

where $Q = \frac{A_0\theta_0 }{\mu_0}-\mu_1+c\alpha.$ Thus, by Equation (3.1) we have $Q < 0$. So, $L_1$ is Lyapunov and the proof is complete.

4.   Endemic equilibrium

In this section we investigate endemic equilibrium points, $X^* = (S^*, E^*, I^*)$. The general technique of the proof of Theorem 4.1 is similar as in the papers ^{[10,13,24,27]} and others.

Theorem 4.1. The system (2.7) has an endemic equilibrium point $X^* = (S^*, E^*, I^*)$ if and only if one of the following conditions (a) or (b) holds, and then the endemic equilibrium point is unique.

(a) $A_1 > 0$,

(b) $A_1 = 0$ and ${\cal R}_0 > 1$.

Proof. For an equilibrium point $X^*$, the following relations hold.

These identities give rise to the equation that must be satisfied by $I^*$ if $I^*$ exists:

(a) Suppose that $A_1 > 0$. Let us introduce the constant

Let $F:(0, \zeta)\to {\mathbb R}$ be the function defined as

Now, since $A_1 > 0$, we note that

A routine calculation gives

with

and the discriminant of $Q(x)$ is $-4\alpha^2\mu_1\mu_2A_0 < 0$. Therefore, $Q(x) > 0$ for all $x\in (0, \zeta)$, and consequently $F'(x) > 0$ for all $x\in (0, \zeta)$.

Since $\theta$ and $F$ are continuous, and their domains are connected (and overlapping), boundedness of $\theta$ together with Eq (4.3) implies that there exists $x_0\in (0, \zeta)$ such that $\theta(x_0) = F(x_0)$. Such $x_0$ is unique since $\theta$ is a decreasing function, and $F$ is an increasing function. The number $x_0$ happens to be $I^*$. This completes the proof of part (a) of the theorem.

(b) Suppose that $A_1 = 0$. Similarly as above, let $\xi = \alpha A_0/(\mu_1\mu_2)$ and let $G:[0, \xi)\to {\mathbb R}$ be the function

Now, $G$ is an increasing continuous function on the interval $[0, \xi)$ with

Since $\theta$ is continuous and decreasing, there exists $x_1\in (0, \xi)$ such that $\theta(x_1) = G(x_1)$ if and only if $\theta(0) > G(0)$. The latter condition is equivalent to $\theta_0 > \mu_0\mu_1/A_0$, which is equivalent to the condition ${\cal R}_0 > 1$. Such $x_1$ is unique, and $I^* = x_1$. This proves (b).

Theorem 4.2. If an endemic equilibrium $X^* = (S^*, E^*, I^*)$ for the model (2.7) exists, then $X^*$ is locally asymptotically stable.

Proof. Let us write $\theta(I^*) = \theta^*$ and $\theta'(I^*) = \theta_3$. At the point $X^*$ we have the Jacobian matrix:

$J(X^*) = \left(\begin{array}{ccc} -\theta^* E^*-\mu_0 & -\theta^* S^* & -\theta_3E^*S^* \\ \theta^* E^* & \theta^* S^*-\mu_1 & \theta_3E^*S^*\\ 0 & \alpha & -\mu_2\\ \end{array}\right)$

and then

where

We note that $\mu_1-\theta^* S^* = A_1/E^* > 0$. Therefore, $g_1 > 0$. Also, $\theta_3\le 0$. Therefore, $g_2 > 0$ and $g_3 > 0$. Now, let us test $g_1g_2-g_3$ for positivity. We note that:

which means that $g_1g_2 - g_3 > 0$. Therefore, $\det(\lambda I_3-J(X^*)$ satisfies the Routh-Hurwitz criterion. Thus, the endemic equilibrium point $X^*$ is locally asymptotically stable.

5.   Optimal control problem

This section is devoted to formulating and solving the optimal control problem. The idea is to design a vaccination schedule and another intervention by way of education, law enforcement, etc., which aim to minimize the number of COVID-19 cases on the one hand and the overall cost of these two strategies on the other hand, over a certain time horizon $[0, T]$.

We introduce two control functions $u_1(t)$ and $u_2(t)$ in such a way that the vaccination rate will be replaced by the function $z_0 u_1(t)$, with $z_0$ being the maximum affordable rate of vaccination, and the contact rate $\theta$ is replaced by $\theta [1-u_2(t)]$. For a finite time horizon $[0, T]$, we have

Thus, $u_2$ is an intervention, by way of education, law enforcement, lockdown, etc., towards reducing the contact rate, while $u_1$ is proportional to the vaccination rate.

The optimal control problem has an objective functional, as in Eq (5.1) below, over a fixed time interval $[0, T].$

which is to be minimized subject to model system (5.2), below.

We use squares in the first two terms of the objective functional for convexity. The constants $C_1$ and $C_2$ represent unit cost. The last term of the integrand in equation (5.1) represents the compartment $I(t)$ that we wish to decrease, and $C_3$ is a weight constant. The control variables $u_1(t)$ and $u_2(t)$ are bounded, and are assumed to be measurable functions of time. We solve our optimization problem by using well-established control theory, such as in ^[15] or the very detailed presentation in ^[16]. We construct the Hamiltonian function, and to this end we introduce Lagrange multipliers $\lambda_1(t), \lambda_2(t)$ and $\lambda_3(t)$, also referred to as the costate variables. In the theorem below, the control variable, the state variables and the costate variables are functions of time, but this dependence is suppressed except where required explicitly.

where

Theorem 5. An optimal solution for the problem of minimizing the objective functional in Eq (5.1), subject to the system (5.2), does exist, and it satisfies the following system of differential equations:

with transversality conditions

Furthermore, the optimal controls are given by

and

Proof. Sufficient conditions for existence of an optimal control, such as, for instance, in ^[9], can readily be checked, and we omit the detail. We continue by applying the Pontryagin maximum principle. The partial derivatives of the Hamiltonian with respect to the different state variables are calculated, and then we obtain the time derivatives of the costate variables as in Eq (5.3).

The functions $u_1^*$ and $u_2^*$ must optimize $J$. Therefore, we start our search for $u_1^*$ and $u_2^*$ by making:

and out of this, $u_1^*$ and $u_2^*$ emerge in the stated forms.

6.   Numerical simulation

6.1.   Parameter values

For a specific application we must declare the form of the function $\theta$, and in the current application we choose a function of the form

The model is applied to South Africa, over time periods when there are (assumed to be) no internal or external disturbances on the model. Let us assume that for South Africa, such a time horizon is the period from 1 February to 1 June of 2021. Another such period would be from 1 August until 30 November of each the two years 2021 and 2022. Parameter values for the model are declared in Table 1. These are mostly obtained from the published literature. Some values are roughly estimated, when information or data is not available. The last parameter to be valuated is $\beta$, which is the most important one in the contact rate. We choose $\beta$ so as to fit the COVID-19 incidence data as recorded at Worldometers, ^[35].

6.2.   Two periods

For South Africa, the active infections show peaks at June/July and December/January for 2020 and 2021, see, for instance, ^[35]. One would attribute these peaks to people clustering densely together, more so than during other times of the year. Such unusual behavior, if it persists for a significant period of time, is not accounted for in a simple model with constant parameters such as the model that we have presented here. When the conditions and assumptions of the model are disturbed, the model with parameters constant and fixed, is no longer accurate or is not applicable. Such disturbances can take the form of external factors such as regulations by governments restricting the movement of people, it can be in the form of influx of carriers of the disease, or it can even be internally produced, for instance, when people tend to have much closer contact than usual, such as at a party, feast or another form of gathering. A disturbance may be such that the model will be able to detect it. When such disturbances occur, it may be possible to modify the model by allowing for seasonality or by some other intervention on the model ^[21,27]. An instance, where ideas from the model were utilized for studying a system during the period of the disturbance, is found in ^[34].

We utilize the model over the following periods,

Period 1: 24 July 2020 – 05 November 2020,

Period 2: 08 January 2021 – 29 April 2021.

Initial state: The coordinates of the initial point are determined as follows. The daily incidence (i.e., the number of daily new cases) and $I_0$ (active cases) are obtained from ^[35]. From these we can get an estimate for $E_0$, using the equation

Daily New Cases $= \ \alpha E_0$.

Then, we make an estimate for $S_0 < P_0 -(E_0 + I_0).$ These values were adjusted slightly to get a better fit on the rest of the daily new cases. We then take our initial state to be

$I_0 = 135 110$, $E_0 = 11500$, $S_0 = 0.998P_0$.

Remark 6.2.1. During Period 1, there were no vaccines in circulation in South Africa. Vaccination of the public only started during the course of Period 2. Thus, the effect of vaccination during the periods in point is regarded as insignificant.

Let us first consider Period 1. The observed incidence of COVID-19 is obtained from the source ^[35]. We note that the model with parameters as in Table 1 makes a good fit, see, Figure 2, on the incidence cases over the first 35 days of Period 1. Then, for the rest of that period, we note that the model under-estimates the observed data.

Our explanation for this under-estimation is as follows. During the extraordinary rise of daily cases, there was awareness and caution. After the peak of incidence, when incidence dropped below a certain level, human behaviour became more risky. If the same alertness level had been maintained at least for the given period, the daily new cases would have dropped below 9 000, versus the observed 19 000. Now, quite remarkably, exactly the same model also fitted the data for Period 2, see, Figure 3. We note that the model output shows an increase towards the end of the period, for both Period 1 and Period 2. Now, if we observe the model output for a few more weeks as in Figure 4, we find that it would stabilize at the average of the observed values in the latter half of the observed period. Unfortunately, in reality, the system suffered another serious perturbation (which our model is not meant to handle), which led to the next wave of the epidemic in South Africa.

6.3.   Optimal control

We run a numerical approximation using hypothetical values for the parameters that were introduced in the control problem. We choose $C_1 = 1$, equal to 1 as a reference value. Then, we argue that the economic costs of lockdown and distancing at the workplace, in businesses and schools are very costly, compared to the cost of vaccination. On that basis, we choose $C_2 = 2$. Finally, the balancing weight $C_3$ is chosen as $C_3 = 2\times 10^{-5}$. Our time horizon for the problem is 180 days. The optimal controls obtained for this particular instance of the control problem are shown in Figures 5 and 6.

These are quite informative. They tell us that the vaccination should be executed at maximum level for most of the time. They also tell us that the control on the contact rate parameter, $\theta$, must be run at maximum intensity over a period of some 60 days, after which it can be relaxed, at approximately 90% of the maximum intensity for another 60 days, and then it will dwindle away towards zero more or less linearly. The graphs of the $E$-class are shown in Figure 7. Here, we have a comparison between the system without control (the broken line) and the optimally controlled situation (the unbroken line). For these two scenarios we have computed the new cases that arise over the 180 days, and they are, respectively, 334970 and 288370. We observe a difference of 46600. We get a more impressive performance of the optimal control if we consider the same parameters, but starting at a different initial state. In this case the intervention averts more than 102 000 new cases. The graphs of the $E$-classes and $I$-classes are shown in Figures 8 and 9, respectively. An explanation for this is that the transmission rate is lower when the $I$-numbers are higher, and this is due to the form of the contact rate function, $\theta$.

7.   Concluding remarks

We have introduced a new disease model which allows for vaccination and inflow of infected. The model is mathematically sound with respect to existence of a unique global positive solution, and its equilibrium points have been shown to be stable. This model with general contact rate can be considered as additional to the class of SEI models that was presented in the paper ^[24] of Sigdel and McCluskey and subsequent contributions by ^[10,27] and others. We have also considered a control problem covering a two-fold intervention: vaccination and reduction of the contact rate. As a sample application we have considered the COVID-19 epidemic in South Africa. Our paper is one of a handful of optimal control works applicable to COVID-19 in South Africa, of which the closest to ours is ^[18]. These papers consider control on only non-pharmaceutical interventions. The current paper is the only one that includes vaccination in the model. Other papers have been more concerned with the initial rise of COVID-19 in the earlier part of 2020, whereas we compare our model output with observed data over two periods subsequent to peaks of COVID-19 incidence. The model was able to highlight an important phenomenon. It revealed how behavior change in terms of social distancing and contact shortly after a peak of COVID-19 incidence can be detrimental to the fight against the disease. More generally, it illustrates how a system can experience a significant perturbation from within itself. Of course with only four compartments, this SEIR model is limited in its accuracy. Furthermore, it may be necessary to introduce more classes if we wish to study the effect of asymptomatic infectives more thoroughly.

Acknowledgments

The authors wish to thank the reviewers for valuable comments which led to an improved manuscript.

Conflict of interest

The authors declare that there is no conflict of interest.