Backward bifurcation and optimal control in a co-infection model for SARS-CoV-2 and ZIKV.

In co-infection models for two diseases, it is mostly claimed that, the dynamical behavior of the sub-models usually predict or drive the behavior of the complete models. However, under a certain assumption such as, allowing incident co-infection with both diseases, we have a different observation. In this paper, a new mathematical model for SARS-CoV-2 and Zika co-dynamics is presented which incorporates incident co-infection by susceptible individuals. It is worth mentioning that the assumption is missing in many existing co-infection models. We shall discuss the impact of this assumption on the dynamics of a co-infection model. The model also captures sexual transmission of Zika virus. The positivity and boundedness of solution of the proposed model are studied, in addition to the local asymptotic stability analysis. The model is shown to exhibit backward bifurcation caused by the disease-induced death rates and parameters associated with susceptibility to a second infection by those singly infected. Using Lyapunov functions, the disease free and endemic equilibria are shown to be globally asymptotically stable for R 0 1 , respectively. To manage the co-circulation of both infections effectively, under an endemic setting, time dependent controls in the form of SARS-CoV-2, Zika and co-infection prevention strategies are incorporated into the model. The simulations show that SARS-CoV-2 prevention could greatly reduce the burden of co-infections with Zika. Furthermore, it is also shown that prevention controls for Zika can significantly decrease the burden of co-infections with SARS-CoV-2.

Introduction

Arbovirus diseases (ARBOD) transmitted by Aedes aegypti, such as zika, dengue and chikungunya and the concurrent circulation of these diseases are of major public health concerns in tropical and subtropical regions. The Coronavirus pandemic caused by the “severe acute respiratory syndrome coronavirus 2” (SARS-CoV-2) has posed serious health challenges in countries with overlapping epidemics, consequently increasing the burden on public health system [1]. This is why, SARS-CoV-2 and arboviruses (ARBOD) epidemics co-occurrence has become a matter of great concern to government and health agencies in tropical regions of the world. Indeed, the resemblance in clinical symptoms of Zika and SARS-CoV-2, especially at the early stages of infection is a great challenge which makes appropriate diagnosis very difficult. Hence, the delay in the administration of an appropriate treatment leads to increase in the spread of infection. [2], [3]. Wrong diagnosis can result in lack of the proper care of the right disease and leads to worst health conditions [1], [4]. Rosario and Siqueira [5], in a recent study, also observed that arboviral infections could have life-threatening implications, such as Guillain-Barré syndrome (GBS), encephalitis, myelitis and others.

Mathematical modeling has become an important tool for studying the dynamics of infectious diseases [6], [7], [8], [9], [10], [11], [12]. Several models have been developed to study the dynamics of SARS-CoV-2 [13], [14], [15], [16], [17]. Atangana [13] developed and analyzed a fractal-fractional model for SARS-CoV-2 to assess the impact of lockdown prior to the advent of vaccination, and showed that effective lockdown strategy was very appropriate to contain the spread of the disease at the onset of the pandemic. Also, Khan and Atangana [14] modeled the dynamics of SARS-CoV-2 with quarantine and isolation. They analyzed the dynamical behavior of the disease by describing the interactions among the bats and unknown hosts. Kolebaje and co-authors [15] modeled the dynamics of COVID-19 in some African countries using a real data. They estimated the basic reproduction numbers for some countries and also presented how the disease could be controlled. In another study, Bonyah and co-workers [16] investigated a fractional optimal control model for COVID-19. They highlighted the importance of different control strategies in mitigating the spread of the disease. Moreover, the stability and optimal analysis for a COVID-19 model with quarantine and media awareness were discussed by the authors in [17].

Numerous mathematical studies have investigated the dynamics of SARS-CoV-2 and its co-infection with other diseases such as dengue [26], HIV [27], diabetes [28], [29], [30], [31], tuberculosis [32], [33] and malaria [34], [35], [36]. Most of the co-infection models in the literature do not include the assumption that susceptible individuals can get incident co-infection with the two diseases (an assumption which is possible for some diseases, and yet always ignored). Since there has not been any model yet to study the co-infection between SARS-CoV-2 and Zika virus, we therefore consider a robust novel mathematical model for the co-interactions between these two diseases, capturing incident co-infection by susceptible individuals. We shall also examine how this assumption could influence the dynamics of a co-infection model.

The major contributions of the paper are highlighted as follows:

The positivity and boundedness of solution of the model are discussed.

The model presented herein is qualitatively analyzed for the occurrence of backward bifurcation.

Using Lyapunov functions, the stability of both the disease free and endemic equilibria are examined, when and , respectively.

Time dependent controls are incorporated into the model and analyzed via the Pontryagin’s principle.

The entire model is simulated to examine the impact of various optimal control strategies on the dynamics of SARS-CoV-2, Zika virus and their co-infections.

Model formulation

At any time , the total human population consists of the following epidemiological states: Susceptible humans , infectious humans with SARS-CoV-2 , infectious humans with Zika virus , humans co-infected with SARS-CoV-2 and zika virus , with denotes infected population recovered from SARS-CoV-2 and zika virus, respectively. The total vector population, at any time , consists of the following states: , denoting susceptible vectors and vectors infected with Zika virus, respectively. Susceptible humans catch SARS-CoV-2 at the rate . Individuals in this state may catch zika virus either from infected humans or vectors at the rate , respectively. Furthermore, since concurrent infection with both diseases is possible, we have assumed that susceptible individuals can get co-infected with SARS-CoV-2 and zika virus at the rate . Human–human-transmission of Zika has been investigated in the literature (see, for example [22]). It is also assumed that the natural death rate for each epidemiological group is . Infected individuals with SARS-CoV-2 can get infected with zika virus at the rate, . Likewise, those infected with zika virus can get infected with SARS-CoV-2 at the rate . The death rates due to SARS-CoV-2, zika or co-infection are given by and , respectively. Moreover, recovery rates from SARS-CoV-2, zika virus or co-infection are denoted by and , respectively. In this model, we have assumed that no reversion or re-infection after recovery either from single or dual infections. Incorporating re-infection could be an extension to the proposed model. Parameters in the model are well defined in Table Table 1.

Table 1

Description of parameters in the model (1).

Parameter	Description	Value	References
η_C	SARS-CoV-2 disease-induced death rate	0.015/day	[18]
η_Z	Zika disease-induced death rate, respectively	0.001	[19]
ζ_C	SARS-CoV-2 recovery rate	115/day	[20], [21]
ζ_Z	Zika recovery rates	0.09−0.15	[22]
η_CZ	Co-infected disease-induced death rate	0.015/day	Assumed
ζ_C	Co-infected recovery rate	115/day	Assumed
Λ_h	Human recruitment rate	4,108,50874.9×365 per day	[23]
Λ_v	Vector recruitment rate	20,000 per day	[19]
β_1	Contact rate for SARS-CoV-2 infection	0.5944	[24]
β_2	Contact rate for zika infection (human to human)	0.0100	[22]
β2h	Contact rate for zika infection (vector to human)	0.43	[22]
β2v	Contact rate for zika infection (human to vector)	0.60−0.75	[22]
β_3	Co-infection contact rate (human to human)	0.200	Assumed
μ_h	Human natural death rate	174.9×365 per day	[25]
μ_v	Vector removal rate	121 per day	[19]
α_1,α_2	Modification parameters	1.0	Assumed

Analysis of the model

We shall now analyze the model qualitatively (1) without considering the controls. We begin with the following:

Positivity of solutions

For the model (1) to be epidemiologically meaningful, it is appropriate to show that all its state variables are non- negative over time. We prove the results below:

Let the initial data be .

Then the solutions, , of the model (1) are non-negative for all time .

See Appendix A “Proof of Theorem 1”.

Boundedness

We claim the following result:

The closed set , with is positively invariant with respect to the model (1) .

Adding all the equations corresponding to the human components of the system (1), we have It follows from (2) that where .

which can be re-written as By applying the comparison theorem [37] and simplifying, we obtain that Therefore, the total human population, as . Following the arguments similar to those given above, the total vector population, . Hence, the system (1) has the solution in . Thus, the given system is positively invariant.

The basic reproduction number of the model

By setting the right-hand sides of the equations in the model (1) to zero, we obtain the disease free equilibrium (DFE) as follows: The stability of the DFE is studied by applying the next generation operator method [38] to the system (1). The transfer matrices are given by where, The basic reproduction number of the model (1) is given by

, where , and are the associated reproduction numbers for SARS-CoV-2, Zika and co-infection of both diseases, respectively and are given by For the sake of simplicity, reproduction number associated with the human-to-human Zika transmission is denoted by , and the reproduction number associated with the vector-to-human-to-vector Zika transmission denotes . Thus, the Zika associated reproduction number can be re-written as

Local asymptotic stability of the disease free equilibrium (DFE) of the model

The DFE, , of the model (1) is locally asymptotically stable (LAS) if , and unstable if .

The local stability of the model (1) is analyzed by the Jacobian matrix of the system (1) evaluated at the disease-free equilibrium, and is given by: The eigenvalues are given by whereas, the remaining eigenvalues are the solutions of the equations given by Applying the Routh Hurwitz criterion, the roots of equations in (8) have negative real parts if and only if and . Thus, the DFE, is locally asymptotically stable if .

Endemic equilibrium points of the model

Suppose the reproduction number .

Also, let the factors playing a significant role in the disease transmission for the proposed model at steady state be denoted by: Then the model (1) will have multiple endemic equilibria , where

Backward bifurcation analysis of the model to assess the impact of incident co-infection

In this section, we shall examine the impact of the assumption of allowing incident co-infection with both diseases

with the help of backward bifurcation analysis. The phenomenon of backward bifurcation, which has been observed in several disease models, is usually characterized by the co-existence of a stable disease free equilibrium and a stable endemic equilibrium when the associated reproduction number of the model is less than unity. The public health implication of the backward bifurcation phenomenon of model (1) is that, the classical epidemiological requirement of having the reproduction number less than unity, although necessary, is no longer sufficient for the effective control of the diseases. The following result is obtained using the approach in [39].

To carry out the analysis, we set the disease induced death rates equal to zero because the disease-induce death rates give rise to backward bifurcation in a vector-host model (see for example, [19]). Note that, the sub-model for SARS-CoV-2 does not undergo backward bifurcation, as already shown by the authors in [36]. Furthermore, some researchers on co-infection models (without the assumption of incident co-infection by susceptible individuals) [36], [40], [41] are of the view that the dynamical behavior of the sub-models usually predict the behavior of the complete co-infection models. This motivates the analysis in this paper, under the assumption of incident co-infection with both diseases.

Setting the disease induced death rates , results in a constant population model, where, with . The resulting model with bilinear incidence rates is thus presented in (11). We are interested to investigate the question that if backward bifurcation occurs or not under this situation. If it does not occur, then the disease induced death will be the actual cause and the dynamics of the sub-model surely influences the complete model. However, if it occurs, we need to determine the parameter(s) which cause(s) such occurrence in addition to disease induced death rate in the co-infection model.

The model (1) with negligible induced death rates , exhibits backward bifurcation if the coefficient given below is positive.

Suppose denote an arbitrary endemic equilibrium of the model. By the following change of variables, the model (1) can be re-presented in the following form Consider the case when . If the contact rate (say) is chosen as a bifurcation parameter, then solving for from we have Similarly, for , we have . Evaluating the Jacobian of the system (11) at the DFE, , we obtain: Using the approach in [39], the matrix has a right eigenvector associated with the zero eigenvalue of given by , where the components are: The non-zero components of the left eigenvector of , satisfying are Using Theorem 4.1 in [39] and computing the non-zero partial derivatives of at the disease free equilibrium, ()), the associated bifurcation coefficients are defined below where, It clearly shows that even under the assumption that the disease induced death rates are zero, the bifurcation coefficients and are positive which indicates the occurrence of backward bifurcation. However, if we set the parameters associated with the susceptibility to infection with a second disease by individuals infected with single infection , then the co-efficient satisfies which rules out the possibility of backward bifurcation in the co-infection model. Thus, it is concluded that, human disease-induced death rates and terms associated with susceptibility to additional infection by singly infected individuals can induce backward bifurcation in the co-infection model for two diseases. Therefore, it is observed that, in addition to disease induced death rates (which induced bifurcation in the sub-model), parameters associated with infection with a second disease cause backward bifurcation in the complete model. Hence, by allowing incident co-infection with both diseases, the dynamics of the sub-model does not always drive or influence the dynamics of the complete co-infection model.

Global asymptotic stability of the disease-free equilibrium of the model (1) for a special case

In the absence of infection with a second disease by singly infected individuals (that is, ), the DFE of the model (1) given by , is GAS in provided that .

See Appendix B “Proof of Theorem 5”

The epidemiological implication of Theorem 5 is; if those already infected with a single disease do not get infected with a second disease, then both SARS-CoV-2 and Zika can be eliminated from the population provided that the threshold quantity, , regardless of the initial sizes of the sub-populations.

Global asymptotic stability of endemic equilibrium of the model (1) for a special case

In the absence of infection with a second disease by singly infected individuals (that is, ), the endemic equilibrium , of the model (1) with is globally asymptotically stable (GAS) in with whenever , where

See Appendix C “Proof of Theorem 6”

The epidemiological significance of Theorem 6 is; if those already infected with a single infection do not get infected with a second disease, and if diseases induced death is negligible, then both SARS-CoV-2 and HBV will persist in the population provided that the threshold quantity, .

Optimal control analysis

It was observed in the preceding sections that the occurrence of backward bifurcation in the model (1) makes the effective control of both diseases difficult in the population. The aim of this section is to incorporate the time dependent controls into the model (1) to obtain the optimal interventions for the elimination of the co-infections. They are defined as follows: : SARS-CoV-2 prevention control, : Zika prevention control, : Control against incident co-infection, and : Control against infection with a second disease, and the optimal control model is given by: subject to the initial conditions Let us consider the following objective function where is the final time. The total cost includes the cost of SARS-CoV-2 and arboviruses preventive measures. We need to find an optimal control such that where is the control set such that are measurable with for . The Hamiltonian is given by:

Existence

We now establish the existence of solution for the optimal control that minimizes the objective functional .

Supposeis defined on the control setsubject to system(15)with non-negative initial conditions at, then there exists an optimal controlsuch that, if the following conditions given in [42] hold:

The admissible control set U is convex and closed.

The state system is bounded by a linear function in the state and control variables.

The integrand of the objective functional in (16) is convex with respect to the controls.

The Lagrangian is bounded below by , where, .

Let be the control set consisting of , and the right hand of (15), that is (Eq. (19) is given in Box I),

To prove Theorem 7, we proceed as follows:

It is obvious to note that the is closed. In addition, consider any two arbitrary elements , where . Then, That is, , and hence is a convex set.

The control system (15) can be expressed as a linear function of control variables with the coefficients as functions of time and state variables: with ( is given in Box II) As the parameters and variables of the model are positive, we have

The optimal control problem’s Lagrangian is given by Let us consider two arbitrary elements , with . and . We now show that, it follows from (20) that, Taking the difference of the two equations given above, we have Thus, gives and hence convexity of .

There exists constants and such that, , , , We now establish the bound on . Note that . As , . Now, Hence,

Suppose the set minimizes over , then adjoint variables , satisfy the adjoint equations with (where the adjoint functions given in the Appendix D “Adjoint functions”) Furthermore,

Considering and the associated solutions , Pontryagin’s Maximum Principle [43] is applied to obtain the following: On the interior of the set, where    (), we have Therefore,

Numerical simulations

In this section, simulations of the control system (15) are carried out. This is done with MATLAB using the forward backward sweep by Runge Kutta method. The quadratic cost functions and are used. The weight constants are assumed as: and .

In the subsequent sections, we investigate the impact of various control strategies on the co-dynamics of both diseases and their co-infection. For the demographic parameters, and , we obtained their values based on the total population and life expectancy of Espirito Santo state, Brazil (a country with co-endemicity of both SARS-CoV-2 and zika virus), which are approximately given by 4,108,508 and 74.9 respectively [23], [25]. The initial conditions used for the simulations are set as follows: .

Strategy A: Impact of SARS-CoV-2 prevention control ()

The simulations of the optimal control system (15) when the strategy that prevents SARS-CoV-2 () is implemented, are depicted in Fig. 1, Fig. 1, respectively. On implementation of this intervention strategy, for , so that , we observe a significant decrease in the total number of individuals infected with SARS-CoV-2, as expected (shown in Fig. 1(a)). Interestingly, this strategy has significant impact on co-infected cases. It is observed that, significant co-infected cases of SARS-CoV-2 and Zika are averted by this control strategy (as depicted in Fig. 1(b)). The control profile for this intervention strategy is depicted in Fig. 1(c), showing very high impact against incident infections either with SARS-CoV-2 or co-infections.

Fig. 1

Impact of SARS-CoV-2 prevention () on individuals in and epidemiological classes. Here, , so that .

Strategy B: Impact of Zika prevention controls ()

The simulations of the optimal control system (15) when the strategy that prevents Zika transmission (from both human and vectors) () is implemented, are depicted in Figs. 2(a), 2(b), 2(c), respectively. Implementation of this intervention strategy, for , so that shows a significant decrease in the total number of individuals infected with Zika virus, as expected (shown in Fig. 2(a)). Also, this strategy has high positive population level impact on co-infected cases, as observed in Fig. 2(b). It is interesting to note that zika prevention strategy averts more co-infection cases than the SARS-prevention strategy. Equally, this strategy caused great reduction in infected vector populations as observed in Fig. 2(c). The control profiles for this strategy are depicted by Fig. 2(d). It is observed that the control has very high effect against Zika transmission.

Fig. 2

Impact of Zika prevention () on individuals in and , and vectors in epidemiological class. Here, , so that .

Strategy C: Impact of control against co-infections with both diseases ()

The simulations of system (15) when the strategy that prevents incident co-infection with SARS-CoV-2 and Zika virus () is implemented, are depicted in Fig. 3, Fig. 3. Implementation of this intervention strategy, for , so that indicates a significant decrease in the total number of individuals co-infected with SARS-CoV-2 and Zika virus as depicted by Fig. 3(a). Thus, to avert co-infection cases in the population, it is not enough to prevent new single cases. Efforts must also be made to prevent individuals from getting infected with concurrent infections. This control strategy also caused a great reduction of infected vector population, as shown in Fig. 3(b). Also, the simulation of the system (15) when the strategy that prevents infection with a second disease by those already infected, is depicted by Fig. 3(c). It is observed that a good number of co-infected cases is averted when this strategy is strictly put in place. The control profiles for prevention of incident co-infections and prevention of an additional infection by singly infected individuals, are presented in Fig. 3, Fig. 3. It is observed that, the control strategy against incident co-infections attains its peak value nearly after 50 days into the simulation period and throughout maintains this value. Moreover, the control against infection with a second disease, was at its peak throughout the simulation period.

Fig. 3

Impact of control against incident co-infections and on individuals in and , and vectors in epidemiological class. Here, , so that .

Conclusion

In this paper, a new mathematical model for SARS-CoV-2 and Zika co-dynamics was presented, incorporating incident co-infection by susceptible individuals which is not common in the existing models. We discussed the impact of this assumption on the dynamics of a co-infection model. The model also incorporated sexual transmission of Zika. The positivity and boundedness of solution of the developed model was investigated, in addition to local asymptotic stability analysis. The model was shown to exhibit backward bifurcation, caused by disease induced death rates and parameters associated with susceptibility to a second infection by those already infected. Employing the Lyapunov functions, the disease free and endemic equilibria were shown to be globally asymptotically stable for and , respectively. To manage the co-circulation of both infections in an effective manner, under an endemic setting, time dependent controls in the form of SARS-CoV-2, Zika and co-infection prevention strategies were incorporated into the model. The simulations showed that SARS-CoV-2 prevention could greatly reduce the burden of co-infections with Zika. Furthermore, it was shown that prevention controls for Zika can significantly reduce the burden of co-infections with SARS-CoV-2. We expect that the findings of this paper will open new avenues of research in this direction.

The model proposed in this paper focused only on SARS-CoV-2 and Zika co-infection, without vaccination or re-infection either with single or both diseases. Incorporating vaccination for SARS-CoV-2 and re-infection for both diseases, one could obtain an extension of the model. Also, the emergence of different variants of SARS-CoV-2 attracts further studies on their co-infections with other diseases, such as Zika, dengue, TB, influenza, Malaria and others. One could therefore, consider multi-strains of SARS-CoV-2 and co-infection with zika virus. Also, more investigations could be carried out to the mathematical (stochastic, agent based modeling, within/intra-host) and epidemiological dynamics of SARS-CoV-2 and Zika co-infection. The current research was not a case study because of insufficient data and information about the co-infection of both diseases. For instance, little is known about infection acquired cross-immunity between both diseases. Not much information is available to answer the question: whether the current available vaccines against SARS-CoV-2 could have any impact on the dynamics of zika virus. Thus, further studies with more reliable data and detailed information about the co-infection of both diseases is viable.

CRediT authorship contribution statement

Andrew Omame: Conceptualization, Formal analysis, Methodology, Writing – original draft, Software. Mujahid Abbas: Writing – original draft, Validation, Writing – review & editing, Supervision. Chibueze P. Onyenegecha: Writing – original draft, Review & editing.

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.