Global dynamics of SARS-CoV-2/malaria model with antibody immune response

1.   Introduction

Coronavirus disease 2019 (COVID-19) is a respiratory disease that emerged in China in late 2019. It is attributed to a virus called severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). It continues to spread and causes new cases and deaths daily around the world. According to the update published on 5 April 2022 by the World Health Organization (WHO)^[1], over 489 million confirmed cases and over 6 million deaths have been reported. The highest numbers of new deaths were registered in the United States of America, the Russian Federation, the Republic of Korea, Germany, and Brazil ^[1]. The coinfection of SARS-CoV-2 with other diseases represents an extra challenge to health care systems. One of the concerns is coinfection with malaria and COVID-19. In fact, many SARS-CoV-2/malaria coinfections have been reported in malaria-endemic countries ^[2]. This condition has increased the need to understand the dynamics of coinfection between malaria and SARS-CoV-2.

SARS-CoV-2 is a single-stranded positive-sense RNA virus ^[3]. It is a part of the family Coronaviridae ^[3]. It enters the host cell through the angiotensin-converting enzyme 2 (ACE2) receptor ^[4,5]. ACE2 is expressed in the heart, gastrointestinal tract, kidney, blood vessels and other organs ^[6]. However, it is expressed at high levels in the alveolar epithelial cells of the lungs ^[6]. The transmission of SARS-CoV-2 occurs through respiratory droplets that carry viral particles or by human-to-human contact ^[7]. Ten COVID-19 vaccines have been approved for use by the WHO, including Novavax, Serum Institute of India (Novavax formulation), Moderna, Pfizer/BioNTech, Janssen (Johnson & Johnson), Oxford/AstraZeneca, Serum Institute of India (Oxford/AstraZeneca formulation), Bharat Biotech, Sinopharm, and Sinovac ^[8]. The progress of mRNA-based vaccines is elaborately discussed in a previous study ^[9]. On October 22, 2020, the U.S. Food and Drug Administration (FDA) approved the antiviral drug Veklury (remdesivir) for the treatment of patients with COVID-19 requiring hospitalization ^[10]. It is used to treat adults and pediatric patients aged 12 years and older (with weight$\geq 40$ kg)^[10].

Malaria is a parasitic disease caused by Plasmodium parasites ^[7,11]. Plasmodium falciparum is the deadliest parasite causing malaria. The other parasites causing malaria are P. vivax, P. ovale, P. malariae and P. knowlesi. In 2020, approximately 241 million confirmed cases of malaria and 627000 deaths were reported worldwide ^[12]. The African Region had 95$\%$ malaria cases and 96$\%$ malaria-related deaths ^[12]. Malaria is transmitted through the bites of infected Anopheles mosquitoes ^[11]. Malaria infection in humans consists of two stages: the liver stage and the blood stage ^[11]. Most of the clinical symptoms occur in the blood stage. During the blood stage, the parasites in the form of merozoites infect red blood cells and replicate within them ^[11]. After cells rupture, 8-32 daughter merozoites are released to infect healthy red blood cells ^[11]. Preventive chemotherapies are used to treat malaria infection and its consequences ^[12]. The WHO also recommended the use of RTS and the S/AS01 malaria vaccine for children living in regions with medium to high P. falciparum malaria transmission ^[12]. In this paper, we concentrate on malaria infection in the blood stage.

SARS-CoV-2 and malaria infections share many symptoms, such as fever, fatigue, difficulty breathing, headache, and myalgia ^[13,14,15]. This similarly may lead to problems in the clinical diagnosis of these diseases or ignoring the possibility of coinfection ^[6,7]. SARS-CoV-2 and malaria also share a close incubation period of 7-14 days for Plasmodium falciparum malaria and 2-17 days for SARS-CoV-2, which increases the chance of coinfection ^[7,16]. Indeed, SARS-CoV-2/malaria coinfection has been recorded in many countries ^[16,17,18]. Some studies mentioned that coinfection potentially increases the severity of COVID-19 ^[2,16,19]. However, a large number of studies have stated that neutralizing antibodies directed against Plasmodium falciparum are effective against SARS-CoV-2 particles, reducing the severity of SARS-CoV-2 infection in coinfected patients ^[6,20,21,22]. Therefore, an understanding of the dynamics of coinfection is crucial to develop functional treatments that minimize the risk of death.

Mathematical models have been considered a robust tool to support biological and medical studies of most epidemics. Malaria models have been widely developed and investigated (see, for example, ^{[23,24,25,26,27,28]}). Similarly, SARS-CoV-2 models have been developed in many studies (see, for example, ^{[29,30,31,32,33,34,35,36,37,38,39]}). Nevertheless, between-host models have attracted more attention than within-host models that study infection within a human body ^[40]. To the best of our knowledge, a SARS-CoV-2/malaria coinfection model has not yet been developed. In this paper, we develop a within-host model of SARS-CoV-2/malaria coinfection. This model inspects the interactions between seven components: healthy red blood cells, infected red blood cells, free merozoites, healthy epithelial cells, infected epithelial cells, free SARS-CoV-2 particles, and antibodies. Using the developed model, we (i) confirm the biological acceptance of solutions by proving the boundedness and nonnegativity, (ii) compute all steady-state points with the corresponding conditions of their existence, (iii) prove the global stability of all steady-state points, and (iv) validate the theoretical results by performing numerical simulations.

The article is arranged as described below. Section 2 provides a description of the developed model. Section 3 proves that all solutions are bounded and nonnegative. Moreover, it computes all steady state points. Section 4 uses Lyapunov functions to prove the global stability of these solutions. Section 5 presents some numerical simulations. Finally, Section 6 discusses the results with some future directions.

2.   SARS-CoV-2/malaria model with an antibody immune response

This section provides a description of the model under consideration. The model consists of seven ordinary differential equations and takes the form

where $X(t)$, $I(t)$, $M(t)$, $Y(t)$, $N(t)$, $V(t)$, and $Z(t)$ denote the concentrations of uninfected red blood cells, infected red blood cells, free merozoites, uninfected epithelial cells, infected epithelial cells, free SARS-CoV-2 particles, and antibodies, respectively. Uninfected red blood cells are recruited at a constant rate $\sigma_{1}$. They are infected by free merozoites at rate $\beta_{m}XM$ and die at rate $d_{1}X$. Infected red blood cells burst to produce $\eta$ merozoites per infected cell and die at rate $d_{2}I$. Free merozoites are eliminated by antibodies at rate $q_{1}MZ$ and die at rate $d_{3}M$. Healthy epithelial cells are generated at a constant rate $\sigma_{2}$, infected by SARS-CoV-2 at rate $\beta_{v}YV$ and die at rate $d_{4}Y$. Infected epithelial cells produce SARS-CoV-2 particles at rate $eN$ and die at rate $d_{5}N$. SARS-CoV-2 particles are eliminated by antibodies at rate $q_{2}VZ$ and die at rate $d_{6}V$. B cells are stimulated to produce antibodies to clear malaria merozoites and SARS-CoV-2 particles at rates $p_{1}MZ$ and $p_{2}VZ$, respectively. Antibodies die at a natural rate of $d_{7}Z$. The meanings of the different parameters are summarized in Table 2.

3.   Basic properties

This section verifies the basic properties of the model (2.1), including the existence, nonnegativity and boundedness. Additionally, it lists all possible steady state points of the model with the corresponding existence conditions.

Theorem 1. Assume that $\phi_{i} > 0$ for $i = 1, 2, 3, 4, 5$. Define the compact set $\Phi = \bigg \lbrace (X, I, M, Y, N, V, Z)\in \mathbb{R}_{+}^{7}:0\leq X(t),$ $I(t)\leq \phi_{1}, 0\leq Y(t),$ $N(t)\leq \phi_{2}, 0\leq M(t)\leq \phi_{3}, 0\leq V(t)\leq \phi_{4},$ $0\leq Z(t)\leq\phi_{5} \bigg \rbrace$. Then, $\Phi$ is a positively invariant set for system (2.1).

Proof. For system (2.1), we obtain

This equation ensures that $\big(X(t), I(t), M(t), Y(t), N(t), V(t), Z(t)\big)\in \mathbb{R}_{+}^{7}$ for all $t\geq 0$ when the initial conditions $\big(X(0), I(0), M(0), Y(0), N(0), V(0), Z(0)\big)\in \mathbb{R}_{+}^{7}$.

We define the following equation to prove the boundedness of solutions:

By taking the derivative with respect to $t$, we obtain the following formula:

where $\gamma_{1} = \min \left \lbrace d_{1}, d_{2} \right \rbrace$. Hence, we obtain

where $\phi_{1} = \dfrac{\sigma_{1}}{\gamma_{1}}$. Thus, $X(t)\leq \phi_{1}$ and $I(t)\leq \phi_{1}$. We define the following equation to prove the boundedness of $Y(t)$ and $N(t)$:

Then, we obtain

where $\gamma_{2} = \min \left \lbrace d_{4}, d_{5} \right \rbrace$. Therefore, the formula becomes

where $\phi_{2} = \dfrac{\sigma_{2}}{\gamma_{2}}$. This equation implies that $Y(t)\leq \phi_{2}$ and $N(t)\leq \phi_{2}$. Finally, we define the following equation to prove the boundedness of $M(t)$, $V(t)$, and $Z(t)$:

Then, we obtain

where $\gamma_{3} = \min \left \lbrace d_{3}, d_{6}, d_{7} \right \rbrace$. Therefore, we have

Thus, $M(t)\leq \phi_{3}$, $V(t)\leq \phi_{4}$, and $Z(t)\leq \phi_{5}$, where $\phi_{3} = \dfrac{\eta d_{2}\phi_{1}}{\gamma_{3}}+\dfrac{eq_{1}p_{2}\phi_{2} }{p_{1}q_{2}\gamma_{3}}$, $\phi_{4} = \dfrac{\eta p_{1}q_{2}d_{2}\phi_{1}}{q_{1}p_{2}\gamma_{3}}+\dfrac{e\phi_{2} }{\gamma_{3}}$, and $\phi_{5} = \dfrac{\eta p_{1}d_{2}\phi_{1}}{q_{1}\gamma_{3}}+\dfrac{ep_{2}\phi_{2} }{q_{2}\gamma_{3}}$.

From the equations listed above, the set $\Phi$ is positively invariant.

Theorem 2. The conditions $\mathcal{R}_{0m} > 0$, $\mathcal{R}_{1m} > 0$, $\mathcal{R}_{p} > 0$ $\mathcal{R}_{0v} > 0$, and $\mathcal{R}_{1v} > 0$ exist such that the model (2.1) has seven steady states when the following conditions are met:

(1) The uninfected steady state $E_{0}$ always exists;

(2) The SARS-CoV-2-free steady state without an immune response $E_{1}$ exists if $\mathcal{R}_{0m} > 1$;

(3) The SARS-CoV-2-free steady state $E_{2}$ exists if $\mathcal{R}_{1m} > 1$;

(4) The malaria-free steady state without an immune response $E_{3}$ exists if $\mathcal{R}_{0v} > 1$;

(5) The malaria-free steady state $E_{4}$ exists if $\mathcal{R}_{1v} > 1$;

(6) The SARS-CoV-2/malaria coinfection immune-free steady state $E_{5}$ exists if $\mathcal{R}_{0m} > 1$ and $\mathcal{R}_{0v} > 1$; and

(7) The SARS-CoV-2/malaria coinfection steady state $E_{6}$ exists if $\mathcal{R}_{p} > 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})}$, $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}} > 1+\dfrac{e\beta_{v}q_{1}\sigma_{2}p_{2}}{q_{2}d_{3}d_{5}(p_{2}d_{4}+\beta_{v}d_{7})}$, and $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}} > 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}}$.

Proof. The steady states of system (2.1) satisfy the following algebraic system:

By solving (3.1), we obtain the following steady states:

(1) The uninfected steady state $E_{0} = (X_{0}, 0, 0, Y_{0}, 0, 0, 0)$, where

Therefore, $E_{0}$ always exists.

(2) The SARS-CoV-2-free steady state without an immune response is designated $E_{1} = (X_{1}, I_{1}, M_{1}, Y_{1}, 0, 0, 0)$, where

where $\mathcal{R}_{0m} = \dfrac{\eta \beta_{m}\sigma_{1}}{d_{1}d_{3}}$. Notably, $X_{1}$ and $Y_{1}$ are positive, while $I_{1}$ and $M_{1}$ are positive for $\mathcal{R}_{0m} > 1$. Thus, $E_{1}$ exists when $\mathcal{R}_{0m} > 1$. The threshold parameter $\mathcal{R}_{0m}$ marks the initiation of malaria infection in the body.

(3) The SARS-CoV-2-free steady state with an immune response $E_{2} = \left(X_{2}, I_{2}, M_{2}, Y_{2}, 0, 0, Z_{2} \right)$. The components are defined as follows:

where $\mathcal{R}_{1m} = \dfrac{\eta \beta_{m}\sigma_{1}p_{1}}{d_{3}(p_{1}d_{1}+\beta_{m}d_{7})}$. We note that $X_{2}$, $I_{2}$, $M_{2}$ and $Y_{2}$ are always positive, while $Z_{2} > 0$ if $\mathcal{R}_{1m} > 1$. Hence, $E_{2}$ exists if $\mathcal{R}_{1m} > 1$. The threshold parameter $\mathcal{R}_{1m}$ determines the initiation of the antibody immune response against malaria merozoites.

(4) The malaria-free steady state without an immune response is defined as $E_{3} = (X_{3}, 0, 0, Y_{3}, N_{3}, V_{3}, 0)$. The components are calculated using the following equation:

where $\mathcal{R}_{0v} = \dfrac{e\beta_{v}\sigma_{2}}{d_{4}d_{5}d_{6}}$. Clearly, $X_{3}$ and $Y_{3}$ are always positive, while $N_{3}$ and $V_{3}$ are positive if $\mathcal{R}_{0v} > 1$. The threshold parameter $\mathcal{R}_{0v}$ marks the initiation of SARS-CoV-2 infection in the body.

(5) The malaria-free steady state is $E_{4} = (X_{4}, 0, 0, Y_{4}, N_{4}, V_{4}, Z_{4})$, where

where $\mathcal{R}_{1v} = \dfrac{e \beta_{v}\sigma_{2}p_{2}}{d_{5}d_{6}(p_{2}d_{4}+\beta_{v}d_{7})}$. Notably, $X_{4}$, $Y_{4}$, $N_{4}$ and $V_{4}$ are always positive, while $Z_{4} > 0$ if $\mathcal{R}_{1v} > 1$. Hence, $E_{4}$ exists if $\mathcal{R}_{1v} > 1$. The threshold parameter $\mathcal{R}_{1v}$ establishes an antibody-mediated immune response against SARS-CoV-2 particles.

(6) The SARS-CoV-2/malaria coinfection immune-free steady state is defined as $E_{5} = (X_{5}, I_{5}, M_{5}, Y_{5}, N_{5}, V_{5}, 0)$, where

The components $X_{5}$ and $Y_{5}$ are always positive. $I_{5}$ and $M_{5}$ are positive if $\mathcal{R}_{0m} > 1$, while $N_{5}$ and $V_{5}$ are positive if $\mathcal{R}_{0v} > 1$. Hence, $E_{5}$ exists if $\mathcal{R}_{0m} > 1$ and $\mathcal{R}_{0v} > 1$. These two conditions are needed to establish a SARS-CoV-2/malaria coinfection.

(7) The SARS-CoV-2/malaria coinfection steady state is defined as $E_{6} = (X_{6}, I_{6}, M_{6}, Y_{6}, N_{6}, V_{6}, Z_{6})$, where

By substituting $Y_{6}$ in the fourth equation of the model (2.1), we obtain

Thus, $V_{6}$ satisfies the following equation:

We define a function $H(V)$ as follows:

where

By evaluating the value of $H(V)$ at $V = 0$, we obtain

where $\mathcal{R}_{p} = \dfrac{e\beta_{v}q_{1}\sigma_{2}+q_{2}d_{3}d_{4}d_{5}}{q_{1}d_{4}d_{5}d_{6}}$. Notably, $H(0) > 0$ if

In addition, we find

Thus, $H\left(\dfrac{d_{7}}{p_{2}}\right) < 0$ if

This result implies that a root $0 < V^{*} < \dfrac{d_{7}}{p_{2}}$ exists such that $H(V^{*}) = 0$. Let $V_{6} = V^{*}$ and note that for $0 < V_{6} < \dfrac{d_{7}}{p_{2}}$ and $\mathcal{R}_{1m} > 1$ ($\mathcal{R}_{1m} > 1$ is naturally satisfied at $E_{6}$ because $E_{2}$ coexists with $E_{6}$ when $\mathcal{R}_{1m} > 1$. However, it will not be stable as concluded from Theorem 5, where $X_{6} > 0$, $I_{6} > 0$, $M_{6} > 0$, $Y_{6} > 0$, $N_{6} > 0$ and $Z_{6} > 0$. Similarly, we determine the third existence condition of $E_{6}$ by formulating a function of $Z$ and extracting the conditions under which a positive root is obtained as follows:

Therefore, $E_{6}$ exists if conditions (3.2), (3.3), and (3.4) are satisfied.

4.   Global properties

This section proves the global stability of all steady states of model (2.1) by selecting proper Lyapunov functions.

We define a function $\Pi_{i}(X, I, M, Y, N, V, Z)$ and let $S^{'}_{i}$ be the largest invariant subset of $S = \bigg \lbrace(X, I, M, Y, N, V, Z)\ | \ \dfrac{d \Pi_{i}}{dt} = 0\bigg \rbrace$.

Theorem 3. The steady state $E_{0}$ is globally asymptoticallystable (GAS) if $\mathcal{R}_{0m}\leq 1$ and $\mathcal{R}_{0v}\leq 1$.

Proof. We define a Lyapunov function as follows:

By taking the time derivative of $\Pi_{0}(t)$, we obtain

In this case, $\dfrac{d \Pi_{0}}{d t}\leq 0$ if $\mathcal{R}_{0m}\leq 1$ and $\mathcal{R}_{0v}\leq 1$. Additionally, $\dfrac{d \Pi_{0}}{d t} = 0$ when $X = X_{0}$, $Y = Y_{0}$, and $M = V = Z = 0$. The solutions tend to $S^{'}_{0}$ which contains elements with $M = V = 0$ and then $\dfrac{d M}{d t} = 0$ and $\dfrac{d V}{d t} = 0$. From the third and sixth equations of the model (2.1), we obtain $I = N = 0$. Therefore, $S^{'}_{0} = \lbrace E_{0}\rbrace$, which guarantees the global asymptotic stability of $E_{0}$ when $\mathcal{R}_{0m}\leq 1$ and $\mathcal{R}_{0v}\leq 1$ according to LaSalle's invariance principle ^[41].

Theorem 4. Suppose that $\mathcal{R}_{0m} > 1$. Then, the steady state $E_{1}$ is GAS if $\mathcal{R}_{0v}\leq 1$ and $\mathcal{R}_{1m}\leq 1$.

Proof. We establish the following Lyapunov function:

Then, we obtain

By applying the steady state conditions at $E_{1}$

and collecting the terms of Eq. (4.1), we obtain

Thus, $\dfrac{d \Pi_{1}}{d t}\leq 0$ if $\mathcal{R}_{0v}\leq 1$ and $\mathcal{R}_{1m}\leq 1$. Furthermore, $\dfrac{d \Pi_{1}}{d t} = 0$ when $X = X_{1}$, $I = I_{1}$, $M = M_{1}$, $Y = Y_{1}$, and $V = Z = 0$. The solutions tend to $S^{'}_{1}$ which has $V = 0$ and then $\dfrac{d V}{d t} = 0$. From the sixth equation of (2.1), we obtain $N = 0$. Hence, $S^{'}_{1} = \lbrace E_{1}\rbrace$. Accordingly, LaSalle's invariance principle ^[41] ensures the global asymptotic stability of $E_{1}$ if $\mathcal{R}_{0m} > 1$, $\mathcal{R}_{0v}\leq 1$ and $\mathcal{R}_{1m}\leq 1$.

Theorem 5. Assume that $\mathcal{R}_{1m} > 1$. Then, the steady state $E_{2}$ is GAS if $\mathcal{R}_{p}\leq 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})}$.

Proof. We define a Lyapunov function

Then, we obtain

After collecting the terms of Eq. (4.2) and using the following steady state conditions at $E_{2}$

we obtain

Now, we must calculate the last term of Eq. (4.3). From Theorem 2, we obtain

Thus, the derivative in (4.3) is rewritten as follows:

Importantly, $\dfrac{d \Pi_{2}}{d t}\leq 0$ if $\mathcal{R}_{p}\leq 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})}$. In addition, $\dfrac{d \Pi_{2}}{d t} = 0$ when $X = X_{2}$, $I = I_{2}$, $M = M_{2}$, $Y = Y_{2}$, and $V = 0$. We easily show that the elements of $S^{'}_{2}$ satisfy $N = 0$ and $Z = Z_{2}$. This result implies that $S^{'}_{2} = \lbrace E_{2}\rbrace$. Therefore, the global asymptotic stability of $E_{2}$ follows LaSalle's invariance principle ^[41] when $\mathcal{R}_{1m} > 1$ and $\mathcal{R}_{p}\leq 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})}$.

Theorem 6. Suppose that $\mathcal{R}_{0v} > 1$. Then, the steady state $E_{3}$ is GAS if $\mathcal{R}_{0m}\leq 1$ and $\mathcal{R}_{1v}\leq 1$.

Proof. We establish a Lyapunov function as follows:

Then, we obtain

Using the steady state conditions at $E_{3}$

the derivative in (4.4) is transformed to

This result implies that $\dfrac{d \Pi_{3}}{d t}\leq 0$ if $\mathcal{R}_{0m}\leq 1$ and $\mathcal{R}_{1v}\leq 1$. Additionally, $\dfrac{d \Pi_{3}}{d t} = 0$ when $X = X_{3}$, $I = 0$, $M = 0$, $Y = Y_{3}$, $N = N_{3}$, $V = V_{3}$, and $Z = 0$. Thus, $S^{'}_{3} = \lbrace E_{3}\rbrace$. Accordingly, LaSalle's invariance principle ^[41] ensures the global asymptotic stability of $E_{3}$ when $\mathcal{R}_{0v} > 1$, $\mathcal{R}_{0m}\leq 1$ and $\mathcal{R}_{1v}\leq 1$.

Theorem 7. Assume that $\mathcal{R}_{1v} > 1$. Then, the steady state $E_{4}$ is GAS if $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\leq 1+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\mathcal{R}_{1v}$.

Proof. We consider the following Lyapunov function:

Then, we obtain

The steady state conditions at $E_{4}$ are written as follows:

Using Eq. (4.6) to collect the terms of Eq. (4.5), we obtain

Hence, $\dfrac{d \Pi_{4}}{d t}\leq 0$ if $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\leq 1+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\mathcal{R}_{1v}$. In addition, $\dfrac{d \Pi_{4}}{d t} = 0$ when $X = X_{4}$, $M = 0$, $Y = Y_{4}$, $N = N_{4}$, and $V = V_{4}$. In this case, $S^{'}_{4} = \lbrace E_{4}\rbrace$. According to LaSalle's invariance principle ^[41], the steady state $E_{4}$ is GAS if $\mathcal{R}_{1v} > 1$ and $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\leq 1+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\mathcal{R}_{1v}$.

Theorem 8. Suppose that $\mathcal{R}_{0m} > 1$ and $\mathcal{R}_{0v} > 1$. Then, the steady state $E_{5}$ is GAS if $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}}\leq 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}}$.

Proof. We consider the following Lyapunov function:

Then, we obtain

The steady state conditions for $E_{5}$ are written as follows:

Using the aforementioned conditions, the derivative in (4.7) is transformed into

We evaluate the last term in (4.8) by computing

Thus, we obtain

Hence, $\dfrac{d \Pi_{5}}{d t}\leq 0$ if $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}}\leq 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}}$. Additionally, $\dfrac{d \Pi_{5}}{d t} = 0$ when $X = X_{5}$, $I = I_{5}$, $M = M_{5}$, $Y = Y_{5}$, $N = N_{5}$, $V = V_{5}$, and $Z = 0$. Thus, $S^{'}_{5} = \lbrace E_{5}\rbrace$. According to LaSalle's invariance principle ^[41], $E_{5}$ is GAS if $\mathcal{R}_{0m} > 1$, $\mathcal{R}_{0v} > 1$, and $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}}\leq 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}}$.

Theorem 9. Suppose that $\mathcal{R}_{p} > 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})}$, $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}} > 1+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\mathcal{R}_{1v}$, and $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}} > 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}}$. Then, the steady state $E_{6}$ is GAS.

Proof. We consider the following Lyapunov function:

By taking the time derivative, we obtain

The equilibrium conditions at $E_{6}$ are written as follows:

Utilizing the aforementioned conditions, Eq. (4.9) becomes

Hence, we have $\dfrac{d \Pi_{6}}{d t}\leq 0$. Furthermore, $\dfrac{d \Pi_{6}}{d t} = 0$ when $X = X_{6}$, $I = I_{6}$, $M = M_{6}$, $Y = Y_{6}$, $N = N_{6}$, $V = V_{6}$, and $Z = Z_{6}$. Thus, $S^{'}_{6} = \lbrace E_{6}\rbrace$. Based on LaSalle's invariance principle ^[41], $E_{6}$ is GAS when $\mathcal{R}_{p} > 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})}$, $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}} > 1+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\mathcal{R}_{1v}$, and $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}} > 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}}$.

All steady states with their existence and stability conditions are summarized in Table 1.

5.   Numerical simulations

In this section, we present some numerical simulations to verify the results of the previous sections. These numerical simulations are performed using MATLAB solver ode45 with the following sets of initial conditions:

(i) $X(0) = 5$, $I(0) = 0.0001$, $M(0) = 0.0002$, $Y(0) = 10$, $N(0) = 0.02$, $V(0) = 0.01$, $Z(0) = 0.1\times 10^{10}$,

(ii) $X(0) = 10$, $I(0) = 0.0001$, $M(0) = 0.002$, $Y(0) = 5$, $N(0) = 0.03$, $V(0) = 0.001$, $Z(0) = 0.2\times 10^{10}$,

(iii) $X(0) = 15$, $I(0) = 0.002$, $M(0) = 0.003$, $Y(0) = 1$, $N(0) = 0.04$, $V(0) = 0.002$, $Z(0) = 0.3\times 10^{10}$.

The numerical simulations of model (2.1) are divided into seven cases according to the global stability of each of the steady states calculated in Theorem 2. In each case, the values of $\beta_{m}$, $\beta_{v}$, $p_{1}$, $p_{2}$, and $d_{7}$ vary, while the remaining parameters are predetermined as described in Table 2. The seven cases are as follows:

(1) We choose $\beta_{m} = 2\times 10^{-10}$, $\beta_{v} = 0.1$, $p_{1} = 3\times 10^{-8}$, $p_{2} = 0.96$, and $d_{7} = 2\times 10^{-1}$ to yield $\mathcal{R}_{0m} = 0.6667 < 1$ and $\mathcal{R}_{0v} = 0.9122 < 1$. This equation implies that the steady state $E_{0} = (10\times 10^9, 0, 0, 22.41, 0, 0, 0)$ is GAS (see Figure 1), which coincides with Theorem 3. This equation represents the case of a healthy individual who does not suffer from either malaria infection or SARS-CoV-2 infection.

(2) We select $\beta_{m} = 2\times 10^{-9}$, $\beta_{v} = 0.1$, $p_{1} = 2\times 10^{-9}$, $p_{2} = 0.96$, and $d_{7} = 2\times 10^{-1}$. Hence, we obtain $\mathcal{R}_{0m} = 6.6667 > 1$, $\mathcal{R}_{0v} = 0.9122 < 1$, and $\mathcal{R}_{1m} = 0.7407 < 1$. Consistent with Theorem 4, the steady state $E_{1} = (1.5\times 10^9, 4.25\times 10^{8}, 7.08\times 10^7, 22.41, 0, 0, 0)$ is GAS (see Figure 2) and represents the case of a patient who only has malaria in the absence of an antibody-mediated immune response.

(3) We take $\beta_{m} = 2\times 10^{-9}$, $\beta_{v} = 0.1$, $p_{1} = 3\times 10^{-8}$, $p_{2} = 0.96$, and $d_{7} = 2\times 10^{-1}$. The threshold parameters are obtained as $\mathcal{R}_{1m} = 4.3478 > 1$ and $\mathcal{R}_{p} = 44.6137 < 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})} = 191.0065$. This calculation implies that the steady state $E_{2} = (6.522\times 10^9, 1.739\times 10^{8}, 6.667\times 10^6, 22.41, 0, 0, 1.607\times 10^{10})$ is GAS (see Figure 3), which supports Theorem 5. At this point, the antibody-mediated immune response is activated to eliminate free merozoites from the blood of patients with malaria.

(4) We set $\beta_{m} = 2\times 10^{-10}$, $\beta_{v} = 0.9$, $p_{1} = 3\times 10^{-8}$, $p_{2} = 0.96$, and $d_{7} = 2\times 10^{-1}$. This gives $\mathcal{R}_{0v} = 8.2099 > 1$, $\mathcal{R}_{0m} = 0.6667 < 1$, and $\mathcal{R}_{1v} = 0.0436 < 1$. This calculation leads to the global asymptotic stability of the steady state $E_{3} = (10\times 10^9, 0, 0, 2.73, 0.1789, 0.008, 0)$ as approved by Theorem 6 (see Figure 4). The patient in this case suffered only from SARS-CoV-2 infection, while the antibody-mediated immune response had not yet been activated.

(5) We choose $\beta_{m} = 2\times 10^{-10}$, $\beta_{v} = 0.9$, $p_{1} = 3\times 10^{-8}$, $p_{2} = 3.9$, and $d_{7} = 1\times 10^{-2}$. These values produce $\mathcal{R}_{1v} = 2.4821 > 1$, $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}} = 0.6895 < 1+\dfrac{e\beta_{v}q_{1}\sigma_{2}p_{2}}{q_{2}d_{3}d_{5}(p_{2}d_{4}+\beta_{v}d_{7})} = 1.0568$. Accordingly, the steady state $E_{4} = (10\times 10^9, 0, 0, 6.775, 0.1421, 0.0026, 1.628\times 10^{8})$ is GAS (see Figure 5). This result is consistent with Theorem 7. At this point, the antibody-mediated immune response is activated to clear SARS-CoV-2 particles from the body of patients with COVID-19, reducing the concentration of SARS-CoV-2 particles.

(6) We consider $\beta_{m} = 4\times 10^{-10}$, $\beta_{v} = 0.9$, $p_{1} = 3\times 10^{-8}$, $p_{2} = 0.96$, and $d_{7} = 0.8$. This gives $\mathcal{R}_{0m} = 1.3333 > 1$, $\mathcal{R}_{0v} = 8.2099 > 1$, and $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}} = 1.338 < 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}} = 1.4272$. Thus, the steady state $E_{5} = (7.5\times 10^9, 1.25\times 10^8, 2.08\times 10^7, 2.73, 0.1789, 0.008, 0)$ is GAS (see Figure 6), consistent with Theorem 8. In this situation, SARS-CoV-2/malaria coinfection occurs in the absence of an antibody-mediated immune response. Coinfection increases the concentrations of malaria merozoites and SARS-CoV-2 particles and worsens the medical condition of the patient.

(7) We take $\beta_{m} = 4\times 10^{-10}$, $\beta_{v} = 3.9$, $p_{1} = 3\times 10^{-8}$, $p_{2} = 0.5$, and $d_{7} = 0.3$. In this case, the threshold parameters are $\mathcal{R}_{p} = 79.2777 > 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})} = 51.2316$, $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}} = 1.3562 > 1+\dfrac{e\beta_{v}q_{1}\sigma_{2}p_{2}}{q_{2}d_{3}d_{5}(p_{2}d_{4}+\beta_{v}d_{7})} = 1.0003$, and $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}} = 1.3358 > 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}} = 1.1601$. Consistent with Theorem 9, the steady state $E_{6} = (8.623\times 10^9, 6.887\times 10^7, 9.994\times 10^6, 4.75, 0.1605, 0.0015, 7.185\times 10^8)$ is GAS (see Figure 7). At this point, the antibody-mediated immune response is directed against both malaria and SARS-CoV-2 infections. This immune response works to reduce the concentrations of both malaria merozoites and viral particles.

6.   Discussion

SARS-CoV-2/malaria coinfection represents a real concern, especially in malaria-endemic areas. Thus, a crucial goal is to address the dynamics of this coinfection. Here, we develop a within-host SARS-CoV-2/malaria coinfection model. This model considers the interactions between uninfected red blood cells, infected red blood cells, free merozoites, uninfected epithelial cells, infected epithelial cells, free SARS-CoV-2 particles, and antibodies. The model has seven steady states as follows:

(1) The uninfected steady state $E_{0}$ always exists. It is GAS when $\mathcal{R}_{0m}\leq 1$ and $\mathcal{R}_{0v}\leq 1$. This point simulates the condition of a person without SARS-CoV-2 or malaria infections.

(2) The SARS-CoV-2-free steady state without an antibody-mediated immune response $E_{1}$ exists if $\mathcal{R}_{0m} > 1$. It is GAS when $\mathcal{R}_{0v}\leq 1$ and $\mathcal{R}_{1m}\leq 1$. This state represents the condition of a malaria monoinfected patient with an inactive immune response.

(3) The SARS-CoV-2-free steady state $E_{2}$ exists if $\mathcal{R}_{1m} > 1$. It is GAS when $\mathcal{R}_{p}\leq 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})}$. Here, the antibody-mediated immune response is activated to eliminate malaria merozoites.

(4) The malaria-free steady state without an immune response $E_{3}$ exists if $\mathcal{R}_{0v} > 1$. It is GAS when $\mathcal{R}_{0m}\leq 1$ and $\mathcal{R}_{1v}\leq 1$. This point simulates the situation of a SARS-CoV-2 monoinfected patient with an inactive immune response.

(5) The malaria-free steady state $E_{4}$ exists if $\mathcal{R}_{1v} > 1$. It is GAS when $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\leq 1+\dfrac{e\beta_{v}q_{1}\sigma_{2}p_{2}}{q_{2}d_{3}d_{5}(p_{2}d_{4}+\beta_{v}d_{7})}$. The immune response is activated in SARS-CoV-2 monoinfected patients.

(6) The SARS-CoV-2/malaria coinfection immune-free steady state $E_{5}$ exists if $\mathcal{R}_{0m} > 1$ and $\mathcal{R}_{0v} > 1$. It is GAS when $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}}\leq 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}}$. Here, coinfection occurs in the absence of an immune response.

(7) The SARS-CoV-2/malaria coinfection steady state $E_{6}$ exists, and it is GAS if $\mathcal{R}_{p} > 1+\dfrac{\eta \beta_{m}\sigma_{1}p_{1}q_{2}}{q_{1}d_{6}(p_{1}d_{1}+\beta_{m}d_{7})}$, $\mathcal{R}_{0m}+\dfrac{q_{1}d_{6}}{q_{2}d_{3}} > 1+\dfrac{q_{1}d_{6}}{q_{2}d_{3}}\mathcal{R}_{1v}$, and $\mathcal{R}_{0m}+\dfrac{e\beta_{m}\sigma_{2}p_{2}}{p_{1}d_{1}d_{5}d_{6}} > 1+\dfrac{\beta_{m}(p_{2}d_{4}+\beta_{v}d_{7})}{\beta_{v}p_{1}d_{1}}$. This state represents the occurrence of SARS-CoV-2/malaria coinfection in the presence of an antibody-mediated immune response.

The numerical results are fully consistent with the theoretical results. We found that SARS-CoV-2/malaria coinfection may be protective, as the shared antibody-mediated immune response works by eliminating SARS-CoV-2 particles from the body. This immune response may cause less severe SARS-CoV-2 infection. This result is consistent with many studies that confirmed the positive effect of the shared antibody immune response ^[6,20,21,22]. However, other studies mentioned an increased risk of death in patients with malaria who are infected with SARS-CoV-2 ^[2,16,19] or that malaria infection does not elicit an antibody immune response against COVID-19 ^[45]. Therefore, further studies are needed to evaluate the effect of coinfection of malaria and SARS-CoV-2, to inspect the role of the immune system during coinfection and to develop better methods to treat coinfected patients. The model investigated in the present study can be developed (i) using real data to estimate the values for the parameters and examine the accuracy of the model by analyzing other Plasmodium parasites responsible for human infection, (ii) testing the effect of time delays that may occur during infection or the production of malaria merozoites and SARS-CoV-2 particles, (iii) expanding to a multiscale model to obtain a deeper understanding of coinfection dynamics ^[40,46,47], (iv) considering a coinfection model with more aggressive variants of SARS-CoV-2 (variants of concern), such as alpha, beta, gamma, delta, lambda, and omicron ^[9,46,47], and (v) incorporating the role of CTLs in killing infected cells.

Acknowledgment

The authors acknowledge that this project was funded by the Deanship of Scientific Research (DSR), University of Business and Technology, Jeddah 21361, Saudi Arabia. The authors, therefore, gratefully acknowledge the DSR technical and financial support.

Conflict of interest

The authors declare there is no conflict of interest.