Reassessment of public awareness and prevention strategies for HIV and COVID-19 co-infections through epidemic modeling

TL;DR

This paper uses epidemic modeling to study how public awareness and interventions affect the spread of HIV and COVID-19 when they co-occur in a population.

Contribution

The paper introduces a co-infection model with optimal control strategies for HIV and COVID-19, revealing bistability and backward bifurcation phenomena.

Findings

The co-infection model's reproduction number is the maximum of HIV and COVID-19 reproduction numbers.

Bistability between disease-free and endemic states can occur in the single-infection COVID-19 model even when the reproduction number is less than one.

Optimal control strategies combining interventions like vaccination and public awareness are critical for reducing co-infection spread.

Abstract

A co–infection model between HIV and COVID-19 that takes into account COVID-19 vaccination and public awareness is discussed in this article. Rigorous analysis of the model is conducted to establish the existence and local stability conditions of the single-infection models. We discover that when the corresponding reproduction number for COVID-19 and HIV exceeds one, the disease continues to exist in both single-infection models. Furthermore, HIV will always be eradicated if its reproduction number is less than one. Nevertheless, this does not apply to the single-infection COVID-19 model. Even when the fundamental reproduction number is less than one, an endemic equilibrium point may exist due to the potential for a backward bifurcation phenomenon. Consequently, in the single-infection COVID-19 model, bistability between the endemic and disease-free equilibrium may arise when the basic reproduction number is less than one. From the co–infection model, we find that the reproduction number of the co–infection model is the maximum value between the reproduction number of HIV and COVID-19. Our numerical continuation experiments on the co–infection model reveal a threshold indicating that both HIV and COVID-19 may coexist within the population. The disease-free equilibrium for both HIV and COVID-19 is stable only if the reproduction numbers are less than one. Additionally, our two-parameter continuation analysis of the bifurcation diagram shows that the condition where both reproduction numbers equal one serves as an organizing center for the dynamic behavior of the co-infection model. An extended version of our model incorporates four different interventions: face mask usage, vaccination, and public awareness for COVID-19, as well as condom use for HIV, formulated as an optimal control problem. The Pontryagin’s Maximum Principle is employed to characterize the optimal control problem, which is solved using a forward-backward iterative method. Numerical investigations of the optimal control model highlight the critical role of a well-designed combination of interventions to achieve optimal reductions in the spread of both HIV and COVID-19.

Figures

Transmission diagram for the model.The Transmission diagram of Eq (3) is constructed to represent the key epidemiological processes and transitions considered in the study. The system of differential equations is subsequently derived based on the flow structure shown in this diagram.

HIV/AIDS-only transmission diagram.Transmission diagram of HIV/AIDS-only model in system (4).

COVID-19-only transmission diagram.Transmission diagram of COVID-19-only model in system (10).

Dynamical response of the HIV-COVID co-infection model (3).Solution was computed for the parameter values given in Table 1, with χ2=0.4, ω=1.1, η=0.5, ψ=1.2, θ=0.1, φ=0.1, β=0.4 and Λ1=Λ2=2.

Continuation of parameter-dependent equilibria of model (3).This experiment conducted with respect to Λ1, Λ2, θ, ν, χ1 and χ2, calculated for the parameter values used in Fig 4. Panels (a)–(f) show in the vertical axis on the left (blue) the variations of the HIV-related compartments of the model. The vertical axis on the right (red) displays the parameter-dependence of the COVID compartment C. In the bifurcation pictures, stable and unstable equilibria are represented by solid and dashed lines, respectively. The bifurcation analysis reveals the presence of several branching points listed as follows: BP1 (Λ2≈3.8405×10−4), BP2 (χ1≈0.01206), BP3 (χ1≈0.01297) and BP4 (χ2≈0.08761). Panel (g) shows time responses of system (3) computed at the test points P1 (χ1=0.0135), P2 (χ1=0.0125) and P3 (χ1=0.0116). In these plots, the red series corresponds to C(t) while the blue series to H(t). In what follows, this color code will be employed to present time plots of the H and C compartments in the same graph.

Codimension-2 of model (3).(a) Continuation of the bifurcation points found in Fig 5(e) with respect to the parameters χ1 and χ2. The diagram shows the continuation of the bifurcations BP3 (ℓ1) and BP2 (ℓ2), which intersect each other at the point BP5 (χ1≈0.01206, χ2≈0.5886). From this point, two further curves of branching points emanate, labeled ℓ3 and ℓ4. The resulting curves divides the parameter space into four regions as indicated in the diagram. Panel (b) shows time plots of system (3) calculated at the selected points P1 (χ1=0.011, χ2=0.56), P2 (χ1=0.011, χ2=0.62), P3 (χ1=0.0128, χ2=0.56) and P4 (χ1=0.0128, χ2=0.35).

Optimal control simulation for Strategy A.Simulation of (a) the effect of the control profiles u1(t), u3(t), and u4(t) on COVID-19 infected individuals (b) the effect of the control profiles u1(t), u3(t), and u4(t) on HIV infected individuals (c) the effect of the control profiles u1(t), u3(t), and u4(t) on AIDS infected individuals (d) the effect of the control profiles u1(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and HIV (e) the effect of the control profiles u1(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and AIDS (f) the control profiles u1(t), u3(t), and u4(t).

Optimal control simulation for Strategy B.Simulation of (a) the effect of the control profiles u2(t), u3(t), and u4(t) on COVID-19 infected individuals (b) the effect of the control profiles u2(t), u3(t), and u4(t) on HIV infected individuals (c) the effect of the control profiles u2(t), u3(t), and u4(t) on AIDS infected individuals (d) the effect of the control profiles u2(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and HIV (e) the effect of the control profiles u2(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and AIDS (f) the control profiles u2(t), u3(t), and u4(t).

Optimal control simulation for Strategy C.Simulation of (a) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on COVID-19 infected individuals (b) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on HIV infected individuals (c) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on AIDS infected individuals (d) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and HIV (e) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and AIDS (f) the control profiles u1(t), u2(t), u3(t), and u4(t).

Equations

Full text

1 Introduction

Human Immunodeficiency Virus (HIV) is a virus that attacks the human immune system, specifically targeting CD4 cells [1]. These cells are essential for the body’s defense ability against infections and illness. Consequently, individuals with HIV are more vulnerable to infections and illnesses compared to healthy individuals. The symptoms experienced by those infected with HIV can range from mild, such as fever, chills, and rash, to severe, including a progressively weakened immune system. If treatment is not received, HIV can progress into Acquired Immunodeficiency Syndrome (AIDS), a more severe and advanced condition where the immune system is severely weakened, making the body susceptible to potentially fatal opportunistic infections and cancers. Early diagnosis and treatment with antiretroviral therapy (ART) are vital in managing the progression of the disease and improving the quality of life for those affected [2]. According to the World Health Organization (WHO), approximately 39.9 million people were living with HIV at the end of 2023, with more than 600,000 deaths recorded in the same year [3]. These figures highlight the ongoing global burden of HIV/AIDS and the critical need for sustained efforts in prevention, early diagnosis, and access to antiretroviral therapy (ART).

In late 2019, the world was struck by a new variant of the coronavirus, known as Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), more commonly referred to as COVID-19. The virus was first detected in Wuhan, China, in late 2019 and rapidly spread worldwide, leading the World Health Organization (WHO) to declare a global pandemic in March 2020 [4]. The virus primarily spreads through respiratory droplets expelled when infected individuals sneeze, cough, talk, or even breathe. It can also spread when susceptible individuals touch contaminated surfaces and then touch their mouth, nose, or eyes [5]. The symptoms of COVID-19 might be moderate or severe. Fever, coughing, exhaustion, loss of taste or smell, and muscle soreness are examples of mild symptoms. Severe symptoms include shortness of breath, confusion, bluish lips or face, and, in some cases, death—particularly among individuals who are immunocompromised or older adults aged 65 years and above. Several methods are available for diagnosing COVID-19, such as Polymerase Chain Reaction (PCR) tests, rapid antigen tests, and imaging techniques like X-rays. Vaccination programs for COVID-19 began worldwide in 2020, offering high efficacy in preventing infection and significantly reducing the severity of symptoms in vaccinated individuals who contract the virus.

As previously explained, HIV can weaken the immune system of infected individuals, making them more susceptible to other diseases. Consequently, co-infection between HIV and COVID-19 is highly likely to occur during the COVID-19 pandemic, as reported in [6, 7]. COVID-19 infected individuals who are already living with HIV appear to have a more severe symptoms to COVID-19, and can even lead to death if not treated properly [8].

Since it was first introduced and formulated in a paper [9], a mathematical model has been used by many authors to model the spread of disease among the population, such as in dengue fever [10–12], malaria [13, 14], tuberculosis [15–19], HIV [20–23], COVID-19 [24–29], and many other types of diseases [30–35]. The concept of the basic reproduction number has been used by these authors as the endemic indicator of their mathematical models. In most cases, the basic reproduction number equal to one becomes the threshold. If the basic reproduction number is larger than one, then there will always exist an endemic equilibrium. On the other hand, there is a chance to eliminate the disease if the basic reproduction number is less than one. The basic reproduction number is defined as the estimated number of secondary cases caused by single primary case in a fully susceptible population. There are several methods that can be used to calculate the basic reproduction number, one of which is the next-generation matrix approach [36].

A mathematical model that describes co-infection between two or more diseases uses a set of mathematical equations to analyze interactions between diseases within the same population. These models often incorporate factors such as disease transmission (both internal and external interactions), recovery, and potential interactions between pathogens. Co-infection models can become complex and require a higher dimension of system due to intricate interactions between diseases. Mathematical models of co-infection involving HIV or COVID-19 with other diseases have been discussed extensively in the literature. Since HIV/AIDS causes immune deficiency, several studies have explored its potential co-infection with other diseases, such as tuberculosis [37–39] and malaria [40, 41]. On the other hand, because COVID-19 emerged only at the end of 2019, co-infection models involving COVID-19 and other diseases have been discussed in more recent studies, such as [42–45] for general co-infections, and specifically for COVID-19 with HIV in [46–50].

The best of our knowledge, the first epidemic model addressing co-infection between HIV and COVID-19 was introduced by the authors in [46], where they developed a six-dimensional system of ordinary differential equations. Their model includes a Susceptible-Infected (SI) compartment for the COVID-19 animal reservoir. Additionally, this article introduces an extended model using fractional-order differential equations. In 2021, the authors in [47] analyzed the model proposed by [46] using an ABC fractional derivative approach. The study explored the phenomenon of backward bifurcation, supported by numerical experiments for the fractional-order model. A more complex co-infection model between HIV and COVID-19 was introduced by the authors in [48]. This model incorporates a vaccinated compartment for COVID-19 and distinguishes between two variants of the virus, the Wild type and the Delta variant, based on their differing infection rates. The article also includes a data fitting process using COVID-19 incidence data and a sensitivity analysis of the model’s reproduction number, employing Partial Rank Correlation Coefficient (PRCC) combined with Latin Hypercube Sampling (LHS). Their findings highlight the significant impact of fractional derivatives on disease dynamics. An optimal control model for the co-infection between HIV and COVID-19 was also discussed by [48]. This eight-dimensional system of ordinary differential equations incorporates vaccination and treatment strategies for COVID-19 and prevention measures for HIV as control variables. The simulations demonstrate that COVID-19 prevention can significantly reduce the burden of co-infections with HIV, and vice versa. Recently, a co-infection model involving HIV, COVID-19, and Monkeypox was explored by the authors in [50]. The study analytically examines single infection models, co-infection between two diseases, and co-infection among three diseases, focusing on the existence of equilibrium points and reproduction numbers. Sensitivity analysis reveals that natural death rates and disease-induced death rates are key parameters influencing the spread of these diseases.

In the post-COVID-19 era, public attention to COVID-19 is no longer as heightened as it was during the pandemic period from 2020 to 2022. COVID-19 booster vaccines are no longer mandatory but have become optional for individuals seeking to enhance their protection against the virus. People who are more aware of COVID-19 are more likely to opt for a vaccine booster. Based on this context, this article presents a co-infection model between HIV and COVID-19 that incorporates population awareness of COVID-19. We conduct a mathematical analysis to examine the existence and stability of equilibrium points, along with the calculation of the reproduction number. Another novel contribution of this study lies in its use of numerical experiments with the continuation software COCO, which is employed to analyze the co-infection model—particularly focusing on the continuation and bifurcation analysis of parameter-dependent equilibria. COCO is chosen for its ability to efficiently track solution branches, detect bifurcation points, and handle large, structured systems, making it well-suited for exploring the rich dynamical behavior inherent in co-infection models. Furthermore, an optimal control extension of the model is analyzed to explore potential strategies for reducing the spread of HIV, COVID-19, and their co-infection. Our optimal control model involved four different interventions, namely the use of face masks, media campaign and vaccination for prevention efforts to reduce the spread COVID-19, and also the use of condom to prevent HIV infection.

The layout of this article is as follows. In Sect 2, we carefully formulate our model based on several key assumptions. The model is constructed as a system of nine-dimensional ordinary differential equations. The model analysis is presented in Sect 3, where we examine the HIV-only model, the COVID-19-only model, and the co-infection model. In Sect 4, we perform a numerical investigation of the co-infection model using COCO. The model extension as an optimal control model is discussed in Sect 5. Finally, Sect 6 provides conclusions and outlines future research directions.

2 Construction of the mathematical model

Let the human population be divided into nine compartments based on their health status and type of disease as follows: the susceptible unaware compartment, Su; the susceptible aware compartment, SA; the vaccinated with COVID-19 compartment, V; the COVID-19-infected compartment, C; the HIV-infected compartment, H; the AIDS-infected compartment, A; the compartment coinfected with COVID-19 and HIV, CH; the compartment coinfected with COVID-19 and AIDS, CA; and the recovered from COVID-19 compartment, R. Since HIV/AIDS is incurable, there is no recovery compartment for individuals with HIV/AIDS. Hence, the total human population, N, is given by:

$$[eqn]$$

The model construction follows the transmission diagram shown in Fig 1, with details described as follows. We assume that HIV/AIDS and COVID-19 are not transmitted vertically to newborns; transmission occurs only through direct or close contact between susceptible individuals and infected individuals. Hence, all newborns are assumed to be susceptible and enter the population through the SU compartment at a constant rate, $[eqn]$ . We incorporate population awareness of COVID-19 into our model. Awareness of HIV/AIDS is not included, as we assume HIV/AIDS is a long-established disease to which the population is already aware. On the other hand, awareness of COVID-19 still needs to be enhanced through media campaigns or other approaches. Susceptible unaware individuals are those who do not fully recognize the dangers of COVID-19 transmission. Consequently, they neither receive the COVID-19 vaccine nor take protective measures against COVID-19 transmission. There is a transition from unaware to aware susceptible individuals at a rate of $[eqn]$ . This awareness is assumed to be temporary; thus, there is a dropout rate from SA to SU at a constant rate, $[eqn]$ . Furthermore, only aware susceptible individuals receive the COVID-19 vaccine at a constant rate, $[eqn]$ .

Transmission diagram for the model.The Transmission diagram of Eq (3) is constructed to represent the key epidemiological processes and transitions considered in the study. The system of differential equations is subsequently derived based on the flow structure shown in this diagram.

HIV/AIDS infection occurs through direct contact between susceptible individuals ( $[eqn]$ , and R) and those infected with HIV/AIDS (H,A,CH and CA). We assume that the probability of a successful infection with HIV/AIDS is $[eqn]$ . Therefore, the transmission term for HIV/AIDS is given by

$$[eqn]$$

where $[eqn]$ is a correction parameter greater than one, reflecting that individuals in compartments A, CH, and CA are more infectious than those in H due to an increased viral load of HIV.

We assume that COVID-19 infection occurs through contact with individuals infected with COVID-19 (C) or those coinfected with both COVID-19 and HIV/AIDS, represented by CH or CA. The standard transmission rate is denoted by $[eqn]$ . Thus, the infection term for COVID-19 is given by

$$[eqn]$$

where $[eqn]$ is a correction parameter for COVID-19 infection due to contact with coinfected individuals. We assume that $[eqn]$ is less than one because individuals infected with HIV/AIDS are unable to engage in normal social activities like those who are not infected with HIV/AIDS. An important feature of our model is that individuals infected with HIV or AIDS have a higher likelihood of contracting COVID-19 due to their weakened immune systems. Therefore, we introduce correction parameters $[eqn]$ and $[eqn]$ , where $[eqn]$ , corresponding to increased susceptibility to COVID-19 for individuals in compartments H and A, respectively. Hence the co–infection term of H and A compartment is given by $[eqn]$ and $[eqn]$ , respectively.

Due to varying levels of awareness and vaccination within the susceptible population, the probability of successful COVID-19 transmission differs between compartments SU, SA, and V. The transmission rate for the vaccinated compartment, V, is the lowest, while SU has the highest transmission rate. Thus, new infections from the SU compartment are given by $[eqn]$ , where $[eqn]$ is the correction parameter for the infection term for SU individuals. For the SA compartment, it is given by $[eqn]$ , and for the V compartment, it is $[eqn]$ , where $[eqn]$ represents COVID-19 vaccine efficacy.

An additional key assumption in our model is that COVID-19 and AIDS may cause death at constant rates, denoted by $[eqn]$ , $[eqn]$ , and $[eqn]$ for individuals in compartments C/CH, A, and CA, respectively. Since COVID-19 is a recurrent disease, we define $[eqn]$ , $[eqn]$ , and $[eqn]$ as the recovery rates for individuals in compartments C, CH, and CA, respectively. Moreover, the transition rate from H to A, representing the progression of HIV infection, is denoted by $[eqn]$ , while the transition rate from CH to CA is represented by $[eqn]$ , with $[eqn]$ .

Based on above description, the mathematical model for co–infection between HIV/AIDS and COVID-19 is given as follows:

$$[eqn]$$

where f1 and f2 is the infection term for HIV/AIDS and COVID-19 given in Eqs (1) and (2), respectively. The description of model parameter is given in Table 1.

Table 1: Parameter definitions for model (3).

3 Model analysis

In this section, we discuss the basic reproduction number of our model, which is a critical threshold parameter determining whether an infectious disease can invade and persist in a population. Given the complexity of the co–infection model, which simultaneously incorporates both HIV/AIDS and COVID-19 dynamics, an analytical examination of becomes highly intricate. To facilitate a more manageable analysis, we begin by exploring the non-co–infection models—specifically, the HIV/AIDS-only model (Sect 3.1 and the COVID-19-only model (Sect 3.2). By deriving and analyzing for these simplified scenarios, we gain valuable insights into the individual disease dynamics. These insights serve as a foundation for anticipating the potential behavioral patterns and outcomes of the more complex co–infection model, thereby enhancing our understanding of its possible dynamic outputs.

3.1 HIV only model

Without COVID-19 co-infection, the total population is given by $[eqn]$ . The transmission diagram of the mathematical model for HIV-only infection, based on model (3), is shown in Fig 2, and the model is presented as follows:

$$[eqn]$$

HIV/AIDS-only transmission diagram.Transmission diagram of HIV/AIDS-only model in system (4).

Since the above model represents the human population, it is necessary that all variables in model (4) are non-negative for all t > 0. This properties is discussed in the following lemma.

Lemma 1. If the initial conditions of system (4) given by $[eqn]$ , then the solutions S**_U_(t),H(t), and A(t) of system (4) are non-negative for all time t > 0.

Proof: In the boundary of $[eqn]$ , we have:

$$[eqn]$$

From the above calculation, it can be seen that all directions on the boundary are always pushed back inward to $[eqn]$ . Hence, since the initial conditions are always non-negative in $[eqn]$ , the solutions SU(t), H(t), and A(t) remain positive for all time t > 0. Therefore, the proof is complete. $[eqn]$

Lemma 2. The feasible region defined by

$$[eqn]$$

is positively invariant and attracting with regard of system (4).

Proof: Since the initial conditions is non-negative, then it is sufficient to prove that SU(t),H(t),A(t) is bounded in $[eqn]$ . By summing all equation in (4) we have:

$$[eqn]$$

Hence, the solution will satisfy $[eqn]$ , implying $[eqn]$ . Furthermore, it is always satisfy $[eqn]$ . Therefore, the differential equation $[eqn]$ has the solution that is always bounded in the domain $[eqn]$ . Hence the proof is complete. $[eqn]$

3.1.1 HIV-free equilibrium point and the basic reproduction number ℛ0HIV.

Taking the right hand side of system (4) equal to zero, and solve it respect to susceptible S_U_by assuming H = A = 0, the disease-free equilibrium for HIV only model in Eq 4 is given by:

$$[eqn]$$

In the context of an HIV infection model, the disease-free equilibrium ( $[eqn]$ ) represents a critical theoretical state where no individuals are infected with HIV in the population — that is, all infection-related compartments (e.g., HIV and AIDS compartment) are at zero, and only susceptible compartment have nonzero populations. Furthermore, $[eqn]$ provides the baseline for assessing control strategies through the calculation of the respected basic reproduction number.

Next we calculate the basic reproduction number of the HIV-only model in (4). The basic reproduction number, representing the potential number of secondary cases caused by a single primary case in a completely susceptible population during the infectious period of the primary case, has been widely used by many researchers to determine the qualitative behavior of their epidemiological models, including whether the disease will disappear or persist. In our case, we calculate the basic reproduction number (denoted by $[eqn]$ ) using the next-generation matrix approach [36]. The transition (Fh) and transmission $[eqn]$ matrix from infected compartment H and A of model (4) is given by:

$$[eqn]$$

Since there is a zero row of $[eqn]$ , then the next-generation matrix formula is given by $[eqn]$ where $[eqn]$ . Hence, the next-generation matrix is given by:

$$[eqn]$$

Therefore, the basic reproduction of the HIV-only model is taken from the spectral radius of Kh, and given by:

$$[eqn]$$

3.1.2 The HIV-endemic equilibrium point.

The second equilibrium of the HIV-only model in (4) is the endemic equilibrium point, where all compartment in model (4) exist in the equilibrium state. This equilibrium is denoted by $[eqn]$ and is given by:

$$[eqn]$$

where:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Based on the above expression, we have the following theorem regarding the existence of HIV-endemic equilibrium point:

Theorem 1. The HIV-endemic equilibrium point of model (4), denoted by $[eqn]$ exist only if $[eqn]$ and not exist otherwise, where $[eqn]$ is the basic reproduction number of the HIV-only model in (4).

3.1.3 The stability of the HIV-free equilibrium point.

The linearization of model (4) in $[eqn]$ is given by:

$$[eqn]$$

We can see clearly that the first column of $[eqn]$ is $[eqn]$ . Hence, $[eqn]$ is the first eigenvalues of $[eqn]$ . The other eigenvalues is taken from the block matrix $[eqn]$ matrix that was left. The characteristic polynomial for the other two eigenvalues is given by:

$$[eqn]$$

Since

$$[eqn]$$

then all coefficient of $[eqn]$ will be positive if $[eqn]$ . Hence, we have the following theorem regarding the local stability criteria of $[eqn]$ .

Theorem 2. The HIV-free equilibrium $[eqn]$ is locally asymptotically stable if $[eqn]$ , and unstable if $[eqn]$ .

3.1.4 The stability of HIV-endemic equilibrium point.

We use Castillo-Song bifurcation theorem [53] to analyze the stability of equilibria around $[eqn]$ . First, we rewrite our model in (4) as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Let the bifurcation parameter is $[eqn]$ evaluated at $[eqn]$ , gives:

$$[eqn]$$

Evaluate the Jacobian matrix of 8a in $[eqn]$ and $[eqn]$ has three eigenvalues, i.e. $[eqn]$ , $[eqn]$ , and $[eqn]$ . Since we have a simple zero eigenvalue, while the other two eigenvalues are negative, then we can use the Castillo-Song bifurcation theorem.

In Castillo-Song bifurcation theorem, we need to calculate two bifurcation type indicators, namely $[eqn]$ and $[eqn]$ using the following recepies:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ is the right and left eigenvector respected to $[eqn]$ . By direct calculation, we have:

$$[eqn]$$ $$[eqn]$$

Substituting $[eqn]$ and wi to $[eqn]$ and $[eqn]$ , yields:

$$[eqn]$$ $$[eqn]$$

Since $[eqn]$ and $[eqn]$ , then we have that system 8a always exhibit a forward bifurcation at $[eqn]$ . Hence, we have the following theorem.

Theorem 3. The HIV-endemic equilibrium point $[eqn]$ of system (4) is locally asymptotically stable for $[eqn]$ , but close to one.

3.2 COVID-19 only model

The specific case of model (3), which considers a population affected solely by COVID-19, is represented by system (10) below. The corresponding transmission diagram is shown in Fig 3.

$$[eqn]$$

COVID-19-only transmission diagram.Transmission diagram of COVID-19-only model in system (10).

In a similar manner to that presented in Lemma 1, the following lemma can be similarly derived.

Lemma 3. If the initial conditions of system (10) satisfy $[eqn]$ and $[eqn]$ , then the solutions $[eqn]$ and R(t) remain non-negative for all t > 0.

The next lemma discusses the feasible region of the COVID-19 model in Eq 10. This theorem can be proved using a similar approach to that used in Lemma 2.

Lemma 4. The feasible region defined by

$$[eqn]$$

is positively invariant and attracting with regard of system (10).

3.2.1 COVID-19-free equilibrium and the basic reproduction number ℛ0COV.

In this section, we calculate the COVID-19-free equilibrium and the basic reproduction number corresponding to the COVID-19 model in Eq 10. By setting the right-hand side of system (10) to zero and assuming C = 0 and R = 0, the COVID-19-free equilibrium is obtained as follows:

$$[eqn]$$

with $[eqn]$ , $[eqn]$ $[eqn]$ , C^0^ = 0, and R^0^ = 0. We apply a similar method as in the previous section to determine the basic reproduction number of the COVID-19 model in (10), denoted by $[eqn]$ . Through direct calculation, we obtain:

$$[eqn]$$

3.2.2 The COVID-19 endemic equilibrium point.

Let $[eqn]$ represent the infection term in the steady-state condition. Hence, Eq 10 in its steady-state condition can be expressed as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Hence, the COVID-19 endemic equilibrium is given by:

$$[eqn]$$

where

$$[eqn]$$

with $[eqn]$ , $[eqn]$ , and $[eqn]$ . Substitute above equation in (??) to $[eqn]$ give:

$$[eqn]$$

Hence, substituting $[eqn]$ and $[eqn]$ to (??) yields the following polynomial for determining the existence of a solution for z:

$$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The polynomial P(z) in (??) is a cubic polynomial in z. According to the Fundamental Theorem of Algebra, it will always have a maximum of three roots. Furthermore, it is evident that since a3>0 and $[eqn]$ , the polynomial P(z) will always have at least one positive root for z if $[eqn]$ . Thus, we arrive at the following lemma.

Lemma 5. The COVID-19-only model in (10) will always have at least one endemic equilibrium if $[eqn]$ .

To determine the possible number of positive roots of P(z), we apply Descartes’ rule of signs, and the results are summarized in Table 2. It is clear that if $[eqn]$ , it is always possible to have an endemic equilibrium $[eqn]$ , either one or three. On the other hand, if $[eqn]$ , while it is possible to have no endemic equilibrium $[eqn]$ , it is also possible to have two $[eqn]$ for $[eqn]$ (refer to rows 3, 5, and 7 of Table 2).

Table 2: Descartes rules of sign for polynomial P(z).

3.2.3 The stability of COVID-19-free equilibrium point.

We apply the linearization technique to analyze the stability of system (10) around $[eqn]$ . By evaluating the Jacobian matrix of system (10) at $[eqn]$ , we obtain:

$$[eqn]$$

where $[eqn]$ .

Three explicit eigenvalues are obtained: $[eqn]$ and $[eqn]$ . It is evident that $[eqn]$ and $[eqn]$ are always negative, while $[eqn]$ is negative if and only if $[eqn]$ . The remaining two eigenvalues are determined from the roots of the following polynomial:

$$[eqn]$$

Since the coefficients of the above polynomial are always positive, the roots will always have negative real parts, regardless of any conditions. Thus, we can state the following theorem.

Theorem 4. The COVID-19 free equilibrium $[eqn]$ of COVID-19 only model in (10) is always locally asymptotically stable if $[eqn]$ and unstable if $[eqn]$ .

3.2.4 The stability of COVID-19 endemic equilibrium point.

We use the center-manifold approach around $[eqn]$ to analyze the local stability of the COVID-19 endemic equilibrium in this section. This method was first introduced by Castillo and Song in [53]. To apply this method, we rewrite the model variables as follows:

$$[eqn]$$

with $[eqn]$ . Hence, the COVID-19 model in (10) now reads as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Next, we choose $[eqn]$ as the bifurcation parameter. By solving $[eqn]$ with respect to $[eqn]$ , we obtain:

$$[eqn]$$

Linearizing system (18) at $[eqn]$ and at COVID-19 free equilibrium point yields:

$$[eqn]$$

with $[eqn]$ , $[eqn]$ , and $[eqn]$ .

From direct calculation, we find a simple zero eigenvalue of $[eqn]$ . Two explicit eigenvalues are negative, namely $[eqn]$ , while the remaining two eigenvalues are determined from the positive roots of the following polynomial:

$$[eqn]$$

Since the coefficients of the above polynomial are all positive, then $[eqn]$ and $[eqn]$ always have negative real parts. Given that there is a simple zero eigenvalue while the other four eigenvalues have a negative real part, we can apply the Castillo-Song bifurcation theorem to analyze the behavior of the system near the branching point at $[eqn]$ .

To determine the type of bifurcation near $[eqn]$ ,we need to evaluate the sign of the following indicator as specified by the Castillo-Song theorem:

$$[eqn]$$

where wk and $[eqn]$ are the right and left eigenvectors of $[eqn]$ correspond to the zero eigenvalue.

From direct calculation, the right eigenvector of $[eqn]$ correspond $[eqn]$ is given by:

$$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Furthermore, the left eigenvector corresponds to $[eqn]$ is given by

$$[eqn]$$

With the left and right eigenvectors at hand, we can directly compute the bifurcation indicators $[eqn]$ and $[eqn]$ . From direct calculation, we obtain:

$$[eqn]$$ $$[eqn]$$

where

$$[eqn]$$

Furthermore,

$$[eqn]$$ $$[eqn]$$

From the above calculation, it is evident that $[eqn]$ is always positive, while $[eqn]$ can be either positive or negative, depending on the sign of P2. Note that P2 can be rewritten in the following form:

$$[eqn]$$

where

$$[eqn]$$

Let $[eqn]$ . Then, we have $[eqn]$ and vice versa. Since $[eqn]$ is always positive, the following theorem describes the type of bifurcation of our COVID-19-only model in (10) at $[eqn]$ .

Theorem 5. The COVID-19 only model in (10) will:

undergoes a forward bifurcation at $[eqn]$ if $[eqn]$ , and

undergoes a backward bifurcation at $[eqn]$ if $[eqn]$ .

From Theorem 5 and the existence criteria of the COVID-19 equilibrium in Table 2, we can conclude that the COVID-19-only model in (10):

always have at least one stable endemic equilibrium for $[eqn]$ but close to one if $[eqn]$ ,always has two endemic equilibria for certain values of $[eqn]$ if $[eqn]$ , where the smaller equilibrium is ununstable nd the larger one is stable.

3.3 COVID-19 and HIV co–infection model

In this final subsection, we analyze the co–infection model of HIV and COVID-19 presented in Eq 3. The analysis includes the criteria for the positiveness of solutions and the calculation of the reproduction number.

Using a similar approach to that applied in the HIV-only and COVID-19-only models, the positiveness criteria for the solutions of the co–infection model between HIV and COVID-19 are provided in the following theorem.

Lemma 6. If the initial conditions of system (3) given by $[eqn]$ then the solutions $[eqn]$ , H(t),A(t),C_H_(t) and C_A_(t) are non-negative for all time t > 0.

The disease-free equilibrium of the co–infection model in (3) is given by:

$$[eqn]$$

with $[eqn]$ , $[eqn]$ and $[eqn]$ Using the next-generation method [53], the basic reproduction number of the co–infection model is given by:

$$[eqn]$$

where $[eqn]$ and $[eqn]$

More discussion of the dynamical properties of the COVID-19 and HIV co–infection model will be discussed in the following section.

4 Numerical study of the COVID-HIV co-infection model

This section will present a detailed numerical investigation of the dynamical behavior of the HIV-COVID co-infection model (3) introduced earlier, when key parameters of the model are perturbed. This study will be performed employing path-following (continuation) methods implemented using the continuation software COCO [54]. This is a numerical platform based on MATLAB designed to solve continuation problems covering, in a broad extent, the analysis and bifurcation detection routines available in classical continuation tools, such as AUTO [55] and MATCONT [56]. In particular, in this work we will use extensively the capabilities for continuation and bifurcation analysis of branches of equilibria in smooth systems of odes, with special emphasis on the detection of codimension-1 and -2 phenomena.

To start our study, we will set first the parameter values of the HIV-COVID co-infection model (3) according to Table 1, considering $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ , which fall within the defined realistic parameter ranges. Under this parameter scenario, the dynamical response of system (3) obtained via direct numerical integration is presented in Fig 4. This diagram shows the time response for selected compartments of the model: SA, V, C, H, CH and CA. The observed behavior reveals an evolution that after transients approaches a full endemic equilibrium, that is, a steady state where both diseases COVID and HIV are present in the system. This equilibrium will be the starting point of our numerical investigation, where one of the main questions concerns how this undesirable steady state is perturbed when key system parameters vary and how to bring this state to a disease-free scenario.

Dynamical response of the HIV-COVID co-infection model (3).Solution was computed for the parameter values given in Table 1, with χ2=0.4, ω=1.1, η=0.5, ψ=1.2, θ=0.1, φ=0.1, β=0.4 and Λ1=Λ2=2.

Next, we will use numerical continuation (path-following) methods to study how this full endemic steady state is affected when selected system parameters are perturbed. The result of this continuation analysis is presented in Fig 5, panels (a)–(f). In this picture, the horizontal axis displays the selected bifurcation parameter, while the vertical axis on the left (in blue) shows the variations of the Euclidean norm of the HIV-related compartments (H,A), which gives a measure of the intensity of the disease in the model. Moreover, the vertical axis on the right (in red) presents the response of the COVID compartment C.

Continuation of parameter-dependent equilibria of model (3).This experiment conducted with respect to Λ1, Λ2, θ, ν, χ1 and χ2, calculated for the parameter values used in Fig 4. Panels (a)–(f) show in the vertical axis on the left (blue) the variations of the HIV-related compartments of the model. The vertical axis on the right (red) displays the parameter-dependence of the COVID compartment C. In the bifurcation pictures, stable and unstable equilibria are represented by solid and dashed lines, respectively. The bifurcation analysis reveals the presence of several branching points listed as follows: BP1 (Λ2≈3.8405×10−4), BP2 (χ1≈0.01206), BP3 (χ1≈0.01297) and BP4 (χ2≈0.08761). Panel (g) shows time responses of system (3) computed at the test points P1 (χ1=0.0135), P2 (χ1=0.0125) and P3 (χ1=0.0116). In these plots, the red series corresponds to C(t) while the blue series to H(t). In what follows, this color code will be employed to present time plots of the H and C compartments in the same graph.

The numerical continuation of the full endemic equilibrium subject to variations to the COVID transmission factor for HIV patients $[eqn]$ is presented in Fig 5(a). As can be observed, an increment of this parameter produces a growth of the COVID compartment, while the HIV-related components decrease. This is due to the fact that larger $[eqn]$ induces a larger transmission rate of COVID, according to the second term of the CH-component of system (3). For this reason, the individuals with only HIV are transferred with a larger rate to the co-infection compartment CH, due to which the HIV-related compartments decrease. An analogous scenario is encountered when the COVID transmission factor for AIDS patients $[eqn]$ is varied, see Fig 5(b). Under the same reasoning, the COVID presence in the system increases, while the number of patients with HIV or AIDS only reduces as $[eqn]$ grows. In this case, however, a critical point BP1 ( $[eqn]$ ) is found, corresponding to a branching. Here, the stable endemic steady state turns into a stable COVID-free equilibrium with HIV in the system.

Fig 5(c) presents the continuation of the full endemic steady state considering variations of the COVID awareness campaign rate $[eqn]$ . Recall that higher $[eqn]$ means increased protection against COVID through vaccination and use of protection equipment such as face masks. Therefore, is is expected to observe a reduction in the COVID compartment as $[eqn]$ increases, see Fig 5(c). At the same time, the compartments related to HIV-COVID co-infection reduce, due to which the population with HIV only increases. An analogous situation is encountered when the COVID vaccine efficacy is perturbed, see panel (d). The COVID infection rate $[eqn]$ is studied in Fig 5(f), showing an increment of the COVID compartment as $[eqn]$ grows. This increases the HIV-COVID co-infections, hence producing a reduction of the population with HIV only. In this diagram, another branching bifurcation is found, this time at BP4 ( $[eqn]$ ). Here, a stable full endemic equilibrium becomes a stable endemic equilibrium with HIV and no COVID. A remarkably different situation is encountered when the HIV infection rate $[eqn]$ is perturbed, see Fig 5(e). A full endemic equilibrium is encountered for large $[eqn]$ , as expected. When this parameter decreases, a critical point BP3 ( $[eqn]$ ) is found, where the COVID compartment becomes zero and only HIV survives in the system. This means that, for the considered biological scenario, the presence of COVID requires a sufficiently high transmissibility of the human immunodeficiency virus. For a different parameter setting, of course, this needs not be the case, as will be seen later. After COVID disappears, a further reduction of $[eqn]$ produces a decrement of HIV until it vanishes from the system at the branching point BP2 ( $[eqn]$ ). Below this point, the disease-free equilibrium gains stability, where none of the diseases are present in the system, a situation that could not be observed by varying the previous parameters.

To conclude this section, we will now carry out a two-parameter study of the critical points encountered above. Specifically, we will perform a numerical continuation of the codimension-1 points BP3 and BP2 with respect to the HIV and COVID infection rates $[eqn]$ and $[eqn]$ , respectively. The outcome of this study is displayed in Fig 6(a). This picture presents the continuation of the bifurcations BP3 ( $[eqn]$ ) and BP2 ( $[eqn]$ ) with respect two parameters, where the resulting curves intersect each other at the point BP5 ( $[eqn]$ , $[eqn]$ ). Here, a degenerate codimension-2 branching phenomenon takes place, where the HIV and COVID endemic branches become zero simultaneously. This point is characterized by the condition $[eqn]$ (see (6) and (??)), which can be achieved, in a generic manner, by varying two parameters (hence codimension-2). The BP5 bifurcation serves as organizing center for the system’s biological behavior, where two further curves of branching points emanate, labeled $[eqn]$ and $[eqn]$ . The computed curves divides the parameter space into four zones as indicated in Fig 6(a). These regions characterize the asymptotic response of the HIV-COVID co-infection model (3) as described next. Region 1 (green): stable disease-free equilibrium; Region 2 (blue): COVID endemic equilibrium without HIV; Region 3 (yellow): HIV endemic equilibrium without COVID; Region 4 (red): full endemic equilibrium with both diseases. This behavior is numerically confirmed by the time plots computed at selected sample points in each region, see Fig 6(b).

Codimension-2 of model (3).(a) Continuation of the bifurcation points found in Fig 5(e) with respect to the parameters χ1 and χ2. The diagram shows the continuation of the bifurcations BP3 (ℓ1) and BP2 (ℓ2), which intersect each other at the point BP5 (χ1≈0.01206, χ2≈0.5886). From this point, two further curves of branching points emanate, labeled ℓ3 and ℓ4. The resulting curves divides the parameter space into four regions as indicated in the diagram. Panel (b) shows time plots of system (3) calculated at the selected points P1 (χ1=0.011, χ2=0.56), P2 (χ1=0.011, χ2=0.62), P3 (χ1=0.0128, χ2=0.56) and P4 (χ1=0.0128, χ2=0.35).

5 Model extension with optimal control

Let us introduce four new control variables into our co–infection model in (3). The first control variable is the use of face masks to reduce the probability of COVID-19 infection. Individuals in SU (the unaware compartment) are assumed to be unwilling to use face masks due to their lack of information about COVID-19. Let $[eqn]$ represent the proportion of individuals using face masks. Furthermore, let $[eqn]$ denote the efficacy of face masks. The first group using face masks consists of susceptible aware individuals SA. Under this assumption, the infection rate of COVID-19 for SA is given by:

$$[eqn]$$

The second group of individuals assumed to be using face masks consists of those in compartment V. Therefore, the infection term for V individuals is given by:

$$[eqn]$$

Individuals infected with HIV are assumed to be highly vulnerable to infections from other diseases. Hence, we assume that they are more aware of COVID-19 and voluntarily willing to use face masks. With this assumption, the infection term for individuals in compartment H with respect to COVID-19 is now expressed as:

$$[eqn]$$

while for individuals in compartment A, it is given by:

$$[eqn]$$

The second control variable in our model is the vaccination rate, which was previously denoted by the constant $[eqn]$ , and is now represented by u2(t). The third control variable in our model is the media campaign rate to increase public awareness of COVID-19, which was previously denoted by the constant $[eqn]$ , and is now represented by u3(t). The fourth control variable in our model, denoted by u4(t) (where $[eqn]$ ), represents the compliance rate to the usage of condom to prevent susceptible individuals from acquiring HIV infection. Furthermore, let $[eqn]$ denote the efficacy of condom usage. Hence, the force of infection for HIV is given as

$$[eqn]$$

With this assumption, the optimal control model of co–infection between HIV and COVID-19 under these control measures is given by:

$$[eqn]$$

For this, we consider the following objective functional

$$[eqn]$$

Here, the final time is denoted by tf, and the weight constants help to balance each term in the integrand (25). The term $[eqn]$ , represents the cost associated with monitoring the infected individuals at all stages. The term $[eqn]$ , denotes the cost associated with the public enlightenment campaign to educate the general public on the dynamics of COVID-19, implementing vaccination, and the cost of acquiring a face mask. The term $[eqn]$ , represents the cost associated with the public health education campaign to educate the general public on safer sex practices.

The goal is to reduce the number of infected individuals as well as the cost associated with the control strategies. The aim is to find an optimal control quadruple u1(t), u2(t), u3(t) and u4(t) such that

$$[eqn]$$

where the control set ( $[eqn]$ ) is defined as

$$[eqn]$$

is Lesbgue measurable.

5.1 Analysis of the optimal control model

An optimal control pair must satisfy the necessary conditions specified by Pontryagin’s Maximum Principle (PMP) [57]. This principle convert (24) and (25) into a problem of minimizing pointwise a Hamiltonian ( $[eqn]$ ) , with respect to the control quadruple u1(t), u2(t), u3(t) and u4(t). The Hamiltonian is given by;

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ are the adjoint functions associated with the state variables of the model (24).

Theorem 6. Given the optimal control sets $[eqn]$ and the solutions $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ of the corresponding state system (24) that minimizes $[eqn]$ over $[eqn]$ , then there exists adjoint functions $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ , such that

$$[eqn]$$

The expression for $[eqn]$ (for $[eqn]$ ) are given in S1 File. Furthermore, the transversality conditions is given by

$$[eqn]$$

The following characterization holds;

$$[eqn]$$

where

$$[eqn]$$

Preposition 1 Corollary 4.2 of (Fleming and Rishel [58] gives the existence of an optimal control sets (u**_1*(t), *u**2*(t), *u**3*(t) and *u**4*(t)) due to the convexity of the integrand of J with respect of (*u**1*(t), *u**2*(t), *u**3*(t) and *u**4_(t)), a prior boundedness of the state solutions, and the local Lipschitz property of the model (24) with respect to the variables.

Proof: Using the Pontryagin’s Maximum Principles, we obtained

$$[eqn]$$

and considering the optimality condition;

$$[eqn]$$

This optimal control sets(u1(t), u2(t), u3(t) and u4(t)) can be solved for subject to the state variables. Taking into account the bounds on the controls, the characterization can be solved as follows;

For the control u1(t), we have

$$[eqn]$$

so that

$$[eqn]$$

For the control u2(t), we have

$$[eqn]$$

Solve it respect to u2 gives

$$[eqn]$$

Next, taking the derivative of $[eqn]$ respect to u3 gives

$$[eqn]$$

Similarly, solve it respect to u3 yield

$$[eqn]$$

For the control u4(t), we have

$$[eqn]$$

Thus, we obtained

$$[eqn]$$

Clearly, the optimality conditions obtained by taking the derivatives of the Hamiltonian with respect to the controls on hold in the interior of the control set. This end the proof. $[eqn]$

5.2 Numerical illustration of the control model

The forward-backward sweep approach is often utilized to determine the optimal control solution. The method begins with an initial estimate of the control variables. Using this initial guess, the state equations are integrated forward in time with a fourth-order Runge-Kutta technique. After the forward integration is complete, the resulting state trajectories, coupled with the initial control estimate, are used to solve the adjoint equations. These adjoint equations are integrated backward in time, starting with the stated terminal conditions. The iterative technique is performed until convergence is reached. This backward time integration is carried out using the reverse fourth-order Runge-Kutta method. The control variables, denoted as u1(t), u2(t), u3(t), and u4(t), are then updated and used to resolve both the state and adjoint systems again. This iterative process continues until satisfactory convergence is achieved in the state, adjoint, and control variables, as described in previous studies [59–62]. The parameter values used in the numerical illustration are specifically those outlined in Table 1. The initial conditions used for the state variables are SU(0) = 9900, SA(0) = 0, $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , CH(0) = 0, and CA(0) = 0. The values of the weight factors used are b1 = 1, b2 = 1, b3 = 1, b4 = 1, b5 = 1, b6 = 1, b7 = 1, b8 = 1, and b9 = 1. The control strategies are combination of efforts involving the following

Strategy A: combination of u1(t), u3(t), and u4(t), while setting u2 to zero.Strategy B: combination of u2(t), u3(t), and u4(t), while setting u1 to zero.Strategy C: combination of u1(t), u2(t), u3(t), and u4(t).

The numerical simulations for Strategy A, B, and C are given in Figs 7, 8, and 9, respectively.

Optimal control simulation for Strategy A.Simulation of (a) the effect of the control profiles u1(t), u3(t), and u4(t) on COVID-19 infected individuals (b) the effect of the control profiles u1(t), u3(t), and u4(t) on HIV infected individuals (c) the effect of the control profiles u1(t), u3(t), and u4(t) on AIDS infected individuals (d) the effect of the control profiles u1(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and HIV (e) the effect of the control profiles u1(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and AIDS (f) the control profiles u1(t), u3(t), and u4(t).

Optimal control simulation for Strategy B.Simulation of (a) the effect of the control profiles u2(t), u3(t), and u4(t) on COVID-19 infected individuals (b) the effect of the control profiles u2(t), u3(t), and u4(t) on HIV infected individuals (c) the effect of the control profiles u2(t), u3(t), and u4(t) on AIDS infected individuals (d) the effect of the control profiles u2(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and HIV (e) the effect of the control profiles u2(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and AIDS (f) the control profiles u2(t), u3(t), and u4(t).

Optimal control simulation for Strategy C.Simulation of (a) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on COVID-19 infected individuals (b) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on HIV infected individuals (c) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on AIDS infected individuals (d) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and HIV (e) the effect of the control profiles u1(t), u2(t), u3(t), and u4(t) on individuals co-infected with both COVID-19 and AIDS (f) the control profiles u1(t), u2(t), u3(t), and u4(t).

5.3 Discussion of results

Fig 7 is the simulation of the effect of control strategy A, which comprises of the usage of face mask (u1(t)), media campaign rate (u3(t)), and the compliance rate to condom usage (u4(t)) on the COVID-19 infected compartment (C(t)), HIV infected compartment (H(t)), AIDS infected compartment (A(t)), co-infected compartment with COVID-19 and HIV (CH(t)), and co-infected compartment with COVID-19 and AIDS (CA(t)). When control strategy A is applied, the peak number of COVID-19 infected individuals decreases from 6800 to 5300 by week 20, as shown in Fig 7(a). For HIV infected individuals, the peak value drops from 1600 to 1350 by week 30, then rises from 200 to 300 by week 80, as depicted in Fig 7(b). Fig 7(c) shows that the number of individuals infected with AIDS increases from approximately 300 to 470 by week 80 under strategy A. Fig 7(d) illustrates that the peak number of individuals co-infected with both COVID-19 and HIV decreases from about 5800 to 5400 by week 30. Additionally, Fig 7(e) demonstrates that the number of individuals co-infected with both COVID-19 and AIDS reduces from around 2200 to 1800 by week 80 when control strategy A is implemented. Fig 7(f) is the corresponding controls u1(t), u3(t), and u4(t)

Fig 8 presents the simulation of the effect of control strategy B, which comprises of the COVID-19 vaccination rate (u2(t)), media campaign rate (u3(t)), and the compliance rate to condom usage (u4(t)) on the COVID-19 infected compartment (C(t)), HIV infected compartment (H(t)), AIDS infected compartment (A(t)), co-infected compartment with COVID-19 and HIV (CH(t)), and co-infected compartment with COVID-19 and AIDS (CA(t)). When control strategy B is implemented, it is observed that the peak number of COVID-19 infected individuals drops from 6800 to 2700 by week 20, as shown in Fig 8(a). The peak number of HIV infected individuals decreases from 1600 to 610 by week 20, as depicted in Fig 8(b). In Fig 8(c), it is seen that the number of AIDS infected individuals falls from 330 to below 150 between week 20 and week 100. Fig 8(d) shows that the peak number of co-infected individuals with both COVID-19 and HIV decreases from around 5800 to 2600 by week 30 under strategy B. Additionally, Fig 8(e) illustrates that the number of co-infected individuals with both COVID-19 and AIDS drops from approximately 2200 to 900 by week 80 when control strategy B is implemented. Fig 8(f) is the corresponding controls u2(t), u3(t), and u4(t).

Fig 9 is the simulation of the effect of control strategy C, which comprises of the usage of face mask (u1(t)), the COVID-19 vaccination rate (u2(t)), media campaign rate (u3(t)), and the compliance rate to condom usage (u4(t)) on the COVID-19 infected compartment (C(t)), HIV infected compartment (H(t)), AIDS infected compartment (A(t)), co-infected compartment with COVID-19 and HIV (CH(t)), and co-infected compartment with COVID-19 and AIDS (CA(t)). When control strategy C is implemented, it is observed that the peak number of COVID-19 infected individuals decreases from 6800 to 2100 by week 20, as shown in Fig 9(a). The peak number of HIV infected individuals reduces from approximately 1600 to 505 by week 20, as depicted in Fig 9(b). In Fig 9(c), it is observed that the number of AIDS infected individuals decreases from 330 to 200 by week 100 under strategy C. Fig 9(d) shows that the peak number of co-infected individuals with both COVID-19 and HIV drops from around 5800 to 2300 by week 30 when strategy C is applied. Additionally, Fig 9(e) illustrates that the number of co-infected individuals with both COVID-19 and AIDS decreases from around 2200 to 750 by week 80 under control strategy C. Fig 9(f) is the corresponding controls u1(t), u2(t), u3(t), and u4(t).

Our results indicate that the most effective strategy to mitigate and control the spread of COVID-19, HIV/AIDS, and their co-infection is the combination of the four control measures (Strategy C). This strategy includes the use of face masks, COVID-19 vaccination, media campaigns to raise public awareness about COVID-19, and compliance with condom usage.

Our optimal control simulations demonstrate that integrated intervention strategies can significantly reduce the morbidity associated with HIV and COVID-19 co-infections. For COVID-19, the combined implementation of face masks, vaccination, and public awareness campaigns substantially decreases transmission, lowers the peak number of infections, and shortens the duration of the epidemic. For HIV, condom use effectively reduces new infections over time. The simulations reveal that simultaneous implementation of all available interventions leads to the greatest impact—minimizing both the individual and combined disease burdens. These results emphasize that partial or isolated interventions are insufficient in controlling co-infections. Instead, long-term, coordinated public health strategies are required to achieve meaningful reductions in disease morbidity. Our findings provide quantitative support for prioritizing comprehensive intervention packages in policy planning, especially in regions affected by syndemic conditions.

6 Conclusion and future works

Mathematical models have been extensively utilized to understand the mechanisms of various disease transmissions within populations, including HIV, COVID-19, and their co-infection. In modeling HIV transmission, numerous studies have examined the impact of antiretroviral therapy [63–65], age-structured populations [21, 22, 66], community awareness [23, 67], co-infections with other diseases [37, 68, 69], and various other factors. For COVID-19, since its emergence as a world-scale pandemic in 2020, a substantial body of scientific articles has focused on mathematical models of its transmission. Topics of interest include vaccination strategies, the effects of population awareness, variations in infection rates, the dynamics of multiple viral strains, and more [70–75]. Co-infection models have also been explored by many authors, particularly regarding HIV or COVID-19 with other diseases. However, only a few studies have specifically addressed mathematical models for co-infection between HIV and COVID-19 [46, 46–48, 48, 50]. According to the explanation above, there are several reasons why studying HIV and COVID-19 co-infection using mathematical models is essential. First, the immune responses of individuals living with HIV may influence the course, spread, and outcomes of COVID-19. A mathematical model provides a framework for quantitatively investigating these interactions in order to gain a better understanding of the dynamics of co-infection. Additionally, mathematical models could evaluate the combined impacts of treatments such as the COVID-19 vaccine, offering insights on the most effective strategies to treat both illnesses simultaneously.

As mentioned before, mathematical models on the co–infection between HIV and COVID-19 are not so many discussed by mathematician. Furthermore, during post pandemic of COVID-19, the importance of community awareness becomes more essential since public awareness on COVID-19 are not so much high as during pandemic era. Hence, it is essential to discuss the impact of public awareness on the contagiousness of COVID-19 considering co–infections with HIV.

Hence, in this study we formulated and analyzed a deterministic epidemic co–infection model between HIV and COVID-19 in a population. The model was constructed as a system of ordinary differential equations. Our discussion focuses on a qualitative investigation of the HIV only model, COVID-19 only model and the co–infection model. Existence and local stability criteria of various types of equilibrium points appearing in the model were derived based on the bifurcation theorem in [53]. The main findings of this study are given as follows.

The model with HIV only. Our work shows that the HIV-free steady state is locally asymptotically stable if the basic reproduction number for HIV, denoted by $[eqn]$ (6), is less than one. Furthermore, $[eqn]$ serves as the threshold, where the stability of the HIV-free steady state changes from stable to unstable when $[eqn]$ , and, at the same time, the HIV-endemic equilibrium begins to emerge. This result emphasizes the importance of keeping $[eqn]$ below one through effective interventions to reduce the infection rate such as widespread of HIV awareness program and behavioral prevention strategies. Once $[eqn]$ exceeds one, the infection becomes self-sustaining in the population, indicating a potential for long-term persistence and increased disease burden. Therefore, identifying and targeting the factors that influence $[eqn]$ is critical for controlling and eventually eliminating HIV transmission.The model with COVID-19 only. Our dynamical analysis shows that the COVID-19-free steady state is locally stable whenever the basic reproduction number for COVID-19, denoted by $[eqn]$ (12), is smaller than one, and unstable if it is bigger than one. The COVID-19 endemic steady state always exists, is unique, and is stable if $[eqn]$ . Furthermore, a backward bifurcation at $[eqn]$ is possible due to the mortality rate induced by COVID-19. Consequently, a higher death rate induced by COVID-19 increases the likelihood of a stable COVID-19 endemic equilibrium even when $[eqn]$ . This implies that simply reducing $[eqn]$ below one may not be sufficient to eradicate the disease if mortality-induced backward bifurcation occurs. Therefore, public health strategies must focus not only on reducing transmission but also on minimizing COVID-19-related mortality through timely treatment, improved health care access, and vaccination. Failure to address this could result in persistent endemicity despite reproduction numbers being nominally below unity.The co–infection model. The basic reproduction number of the co-infection system is given by the maximum value between $[eqn]$ and $[eqn]$ (23). Our numerical experiments using the software COCO reveal that the disease-free equilibrium (free of HIV and COVID-19) is locally stable if both $[eqn]$ and $[eqn]$ are smaller than one. The single-disease endemic equilibrium may be stable only if $[eqn]$ and $[eqn]$ exceeds one. Furthermore, the condition where $[eqn]$ and $[eqn]$ both equal one serves as the organizing center for the dynamical behavior. Small changes around this condition can lead to significantly different dynamics in the co-infection model. This implies that even marginal increases in either reproduction number could trigger the persistence of one or both diseases in the population. As a result, public health strategies must be integrated and responsive to the co-infection context, as control measures for one disease may not be sufficient if the other remains above threshold. Coordinated intervention is essential to avoid overlapping epidemics and prevent a potential syndemic scenario.The optimal control model. The study underscores the significance of effective control measures for HIV and COVID-19 co-infections. Face masks, vaccination, public awareness, and condom use are all effective approaches for reducing infection rates. Combined actions have a greater impact than isolated efforts, with vaccination and mask use reducing COVID-19 spread and condom use limiting HIV transmission. Public awareness enhances the effectiveness of interventions, emphasizing the importance of coordinating initiatives. Numerical results show that long-term, well-planned efforts can result in infection elimination or stabilization of infections. These findings suggest that integrated, multi-faceted public health campaigns are essential to effectively manage co-epidemics. Neglecting one aspect of control—such as awareness or preventive behavior—could undermine the success of even well-funded medical interventions.

Although our model serves an important insight to the co–infection phenomena between HIV and COVID-19 under some important intervention scenarios, we believe that our model still can be extended on several different direction. First direction is by exploring the model by considering varying levels of awareness among different sub-populations, such as adults versus children or urban versus rural communities. The next possible direction is to investigate a non-standard infection term to account for possible saturation effects in infection rates to accommodate people response on the disease transmission in a community. Incorporating stochastic elements into the model to capture random fluctuations in disease transmission dynamics also can be another option for future research. Last, we believe that utilizing more comprehensive datasets, particularly for HIV, to enhance parameter estimation and model validation are necessary for future research direction.

Supporting information

S1 FileThe expression of the co-state system for the optimal control model.(PDF)