An SIS sex-structured influenza A model with positive case fatality in an open population with varying size

1.   Introduction

Influenza A is a highly contagious, respiratory, infectious, viral disease and is potentially lethal. The virus is transmitted from human to human through direct contacts (mainly by coughing, sneezing, or close contact). The genetic materials of influenza viruses are composed of single-stranded RNA, and frequent mistakes happen while copying themselves. Due to their high mutation rates, influenza viruses evolve rapidly. Some of their new generations quickly adapt with the new conditions, which helps them succeed in causing new epidemics and (sometimes) pandemics. Therefore, influenza is both seasonal and pandemic. Seasonal influenza infects more than one billion people annually ^[1]. Recently, it was estimated that seasonal influenza is associated with 294,000 to 518,000 annual respiratory deaths ^[2]; according to the World Health Organization (WHO) ^[1], the estimation is even larger. Moreover, Chaves et al. ^[3] found that influenza is associated with an increased risk of ischemic heart disease (IHD) mortality. The authors attributed the death of about 300,000 IHD adults of ages over 50 years old (every year and globally) to influenza. In their study, they reported a $4\%$ reduction in the global IHD deaths if no influenza was present ^[3]. Therefore, influenza A viruses do significantly impact human health and, consequently, the global economy ^[4].

It is evident that sex is a risk factor to influenza incidences and outcomes, as both males and females differ in their respones to the infection ^[5]. This disparity may reflect genetic and hormonal differences between the two sexes ^[5]. Based on the repeated influenza outbreaks and pandemics, the morbidity and mortality of females are significantly different from those of males ^[6]. The extent to which they are different is associated with other risk factors, including age and the chronic medical conditions from which they suffer ^[5,6]. However, for the purpose of this work, we will focus on mathematically studying the role of sex and gender risk factors on the transmission dynamics of influenza.

Mathematical models have been extensively used to study the dynamics of influenza A infection at both the population level ^{[7,8,9,10,11]} and at the cellular level ^[12,13]. This study is designed to help predict the spread of influenza A infection and to estimate its burden rather than studying the dynamics of cellular interaction and disease progression; therefore, we focus on models that treat the problem at the population level. In this respect, the literature shows the publication of various works that focus on exploring the dynamics in both the absence and in the presence of influenza interventions. For example, Casagrandi et al. ^[7] extended the classical SIR model by including a fourth compartment that represented an intermediate state between the fully susceptible state (S) and the fully protected state (R). This intermediate state is called cross-immune and is denoted by C. Therefore, an SIRC model was developed and mathematically analyzed. This SIRC model describes the dynamics of influenza A infection in a demographically stationary closed population in the absence of influenza-induced mortality, with the fact that exposing cross-immune individuals to a different strain (other than the running one) of influenza either boosts their immunity or reinfects them. The model assumes that the natural immunity acquired by experiencing the infection decays with time and, in consequence, R-compartment individuals move to the C-compartment. Additionally, the model assumes that a further decay of the immunity of C-compartment individuals moves them to the fully susceptible state. Also, motivated by spanish flu data in Australia, Samsuzzoha et al. ^[14] employed two deterministic models (precisely, a Susceptible-Infected-Recovered (SEIRS) and a Susceptible-Vaccinated-Exposed-Infected-Recovered-Susceptible (SVEIRS)) to capture the main characteristic of influenza transmission. The authors fitted the models to the data and estimated the parameters involved in both systems.

Other mathematical models were employed to include quarantine as a preventive intervention to control influenza. Vivas-Barber et al. ^[10] focused on exploring the role of asymptomatic (mild) infections on the long-term transmission dynamics of influenza, in the presence of a fully protective quarantine-based intervention, by introducing a compartment (A) to denote asymptomatically infected individuals to an extended SIQR model (here, Q denotes quarantine individuals). Therefore, an SAIQR model was obtained. The model assumes that the quarantine intervention is completely perfect, while the asymptomatic individuals are capable of transmitting (although less infectious than I-individuals) influenza to the totally susceptible individuals (S). Moreover, the model considers an open population, but ignores the infection lethality. The authors performed a standard analysis of the model and numerically showed (with parameter values typical to influenza) the existence of damped oscillations that described recurring epidemics. Moreover, Erdem et al. ^[11] extended the SIQR model (with perfect quarantine) to consider the fact that no quarantine program is completely perfect. Therefore, a so-called "$\sigma$-quarantine" SIQR model was introduced and mathematically analyzed, where $\sigma \in [0, 1]$ denotes the effectiveness of a quarantine-based intervention program. Other models that were used to describe the dynamics of influenza A infection include (but not limited to) those published by Nuño et al. ^[15], Alexander et al. ^[16], and Krishnapriya et al. ^[17]. However, the role of sex and gender disparities on the influenza dynamics is less explored mathematically.

The classical SIR model has been extended in a previous work to include sex and gender disparities to describe the dynamics of a single outbreak of influenza ^[4]. The model was mathematically analyzed and was employed to assess the impact of these disparities on the influenza A disease outcomes (i.e., the basic reproduction number and the prevalence of the infection at the endemic situation). The author further extended the model to include the impact of applying an imperfect vaccine-based intervention strategy to contain/eliminate influenza A infection. Basically, the author considered a demographically stationary closed population, with the assumption that both male and female sub-populations were equally recruited solely by births. Moreover, the author neglected the possible repeated exposure to influenza A infection as well as its lethality. However, in this work, we consider various biologically meaningful extensions. First, we extend the classical Susceptible-Infected-Susceptible (SIS) model to include the sex and gender disparities in an open population, where the individuals are recruited by births as well as immigrants, where the immigration rate is assumed to depend on the total population size and some kind of carrying capacities. Additionally, the influenza-induced mortality is taken into account, with the assumption that both the male and female sub-populations have different case fatalities. Therefore, the model will be formulated and proven to be well-posed in Section 2. The equilibrium and bifurcation analyses, where the model is shown to exhibit backward bifurcation, is presented in Section 3. Motivated by simulations, the asymptotic stability of the endemic equilibria is presented in Section 4. The endemic prevalence of the infection in the overall population and in each sex-structured sub-population are given in Section 5. The paper closes with a summary and conclusion in Section 6.

2.   Model formulation and well-posedness

2.1.   Model description

2.1.1.   Population demographic model

The total population is assumed to be structured based on the individuals' sex and gender into males (having subscript $1$) and females (having subscript $2$). Denoting the total population size at time $t$ with $N(t)$, the total male (female) population is denoted by $N_1(t)$ ($N_2(t)$), so that $N(t) = N_1(t) + N_2(t)$. It is assumed that new births occur at an average rate $\nu$ from the females; therefore, the total number of births at time $t$ is $\nu N_2(t)$, where a proportion ${ q_0 }$ of them are females. Moreover, the natural death rate is assumed to be $\mu$. Then, in the absence of infection, the demography of the closed population, is governed by the following system of equations:

If this population is assumed to be demographically stationary, in the sense that its size remains roughly constant over time, then $\nu N_2 (t) = \mu N (t)$. Furthermore, assume that more vacancies are created due to developmental reasons so that the population becomes open and new demographic recruitments through immigration are allowed at an average (per unit time) number of new recruitments $\Lambda (K - N (t)) / K$, where $K$ is a kind of carrying capacity, $K - N (t)$ represents the extended extra opportunities (more than needed to the original population), and $\Lambda$ is the maximum rate of immigration. If $q_2$ is the proportion of females amongst the immigrants, then the population dynamics, in the absence of infection, is described by the following dynamical system:

It is noteworthy that the human carrying capacity refers to the number of people a place can sustainably support. It mainly depends on the size of the population, the availability of resources, and how the people use these available resources. Knowing the carrying capacity is vital for sustainable growth without major setbacks due to environmental degradation. In practice, developed countries use the advancement of technology to increase the carrying capacity of a place, to expand economies of scale, and to use the natural resources in an efficient and effective way without degrading the natural environment ^[18]. Therefore, pull factors are created to pave the way for immigrations. The closer the population size is to the carrying capacity, the lower the immigration rate. This is clear in the right hand side of system (2.2).

2.1.2.   Endemic infection model

We are interested in modeling the transmission dynamics of a potentially lethal non-immunizing respiratory infection (with special reference to influenza A) in an open population, whose size changes over time. Therefore, both the male and female populations are assumed to split into two mutually exclusive categories (for each) according to the individuals' epidemiological status: susceptible (with time-dependent size $S(t)$) and infected (with time-dependent size $I(t)$). Precisely, the number of susceptible males and females at time $t$ are denoted by $S_1(t)$ and $S_2(t)$ respectivel, while the number of infected males and females are given by $I_1(t)$ and $I_2(t)$, respectively, so that $N_1(t) = S_1(t) + I_1(t)$ and $N_2 (t) = S_2 (t) + I_2(t)$. A schematic diagram for the transfer between the subpopulation model states is shown in Figure 1, and a brief description of the model state variables is presented in Table 1. Accordingly, new recruitments (either by births or immigrations) are susceptible.

Definitely, susceptible males are recruited due to births at the rate $(1 - q_0) \mu N (t)$ and due to immigration at the rate $(1 - q_2) \Lambda (1 - N (t) / K)$, while they die naturally at the rate $\mu$ and acquire the influenza A infection at a male force of infection $\lambda_1(t)$ and become infected. Infected males either die due to the infection at the rate $c_1 \gamma$ or recover from the infection (without acquiring immunity) and become susceptible again at the rate $(1 - c_1) \gamma$, where $c_1$ is the male case fatality. Hence, the spread of the infection in the male population is described by the following equations:

Similarly, susceptible females are recruited due to births at the rate ${ q_0 } \mu N (t)$ and due to immigration at the rate $q_2 \Lambda (1 - N (t) / K)$, while they decline either due to natural death at the rate $\mu$ or due to acquiring the infection and becoming infected and capable of transmitting the infection at an incidence rate $\lambda_2(t)$. Infected females either die naturally at the rate $\mu$ or die due to the infection at the rate $c_2 \gamma_2$, where $c_2$ is the infection case fatality in females and $\gamma_2$ is the removal rate of infected females. They recover without immunity at the rate $(1 - c_2) \gamma_2$. Therefore, the dynamics of the female population is described by the following two differential equations:

2.1.3.   Forces of infection and removal rates

The male and female forces of infections, $\lambda_1 (t)$ and $\lambda_2(t)$, respectively, are derived in a way similar to that shown in Safan ^[4]. More precisely, we denote the average number of contacts that an individual has with other individuals in the total population by $\tilde\beta$. Additionally, $g_1$ is the susceptibility of males (which is proportional to the probability of success that a susceptible male becomes infected due to contacts with infected males and females), while $g_2$ is the susceptibility of females (which is proportional to the probability that a susceptible female becomes successfully infected as a result of contacts with infected individuals). Assume that $r_{ij}$ accounts for the transmissibility of influenza from infected individuals in the population $i$ to susceptible individuals in the population $j$, for all $i, j \in \{1, 2\}$. Hence, the rate at which susceptible males acquire the infection $\lambda_1(t)$ and that at which susceptible females acquire infection $\lambda_2(t)$ read as follows:

Furthermore, assume that the infected males transmit the infection to either males or females at the same potential, in the sense that the effective rate at which infected males infect susceptible males/females is equal (i.e., $r_{11} = r_{12} = r_1$). Similarly, females transmit the infection to either males or females at the same potential, in the sense that the effective rate at which infected females infect susceptible males/females is equal (i.e., $r_{21} = r_{22} = r_2$). Hence, $\beta = g_1 r_{11} \tilde\beta = g_1 r_1 \tilde\beta : = g_1 \bar\beta$ accounts for the effective contact rate at which susceptible males acquire influenza and $g_2 r_{12} \tilde\beta = g_2 r_1 \tilde\beta = g_2 \bar\beta = (g_2/g_1) g_1 \bar\beta : = g \beta$ is the effective contact rate at which susceptible females acquire the influenza virus. Here, $g = g_2/g_1$ accounts for the relative susceptibility of females with respect to males. Hence, (2.5) reads as follows:

where $r = r_2 / r_1$ is the relative transmissibility of females with respect to male infections.

The removal (either by recovery or due to disease-induced death) rates $\gamma_1$ and $\gamma_2$ are rescaled by assuming that $\gamma_1 = \gamma$ and $\gamma_2 = a \gamma$, where $a$ is a dimensionless rescaling parameter which accounts for the relative recoverability of females with respect to males. A brief description of the model parameters is shown in Table 2.

2.2.   Overall mathematical model

Motivated by the above-detailed assumptions, the populations' overall dynamics is described by the following system of ordinary differential equations:

where

is the male force of infection and the model is defined on the set

It is worth assuring that the letter ${}^T$ in (2.13) denotes a vector transpose.

2.3.   Well-posedness of the model

The model is well-posed in the sense that solutions starting with the initial conditions $(S_1(0), I_1 (0), S_2(0), I_2(0), N(0))^T \in \Omega$ remain in $\Omega$ for all positive times. The following proposition, whose proof is deferred to Appendix 10.1, summarizes the results on the existence and uniqueness of the above model time-dependent solutions.

Proposition 1. The set $\Omega$ is positively invariant and attracts all solutions in $\mathbb{R}_+^5$. In particular, any solution $(S_1(t), I_1(t), S_2(t), I_2(t), N(t))^T$ of the dynamical system (2.7)–(2.11), starting with non-negative initial values $(S_1(0), I_1 (0), S_2(0), I_2(0), N(0))^T \in \Omega$, remains in $\Omega$ for all $t \ge 0$ and is unique.

3.   Equilibrium and bifurcation analysis

To analyze the equilibrium for models (2.7)–(2.11), we put the derivatives in its left hand side equal zero and solve the resulting nonlinear algebraic system of equations, together with (2.12), in the model-state variables. Particularly, the Eqs (2.8) and (2.10) imply that the number of infected males and females, at equilibrium, are given respectively by the following:

Unless otherwise stated, in (3.1) and throughout the rest of this work, the quantities $S_1, I_1, S_2, I_2, N_1, N_2, N,$ and $\lambda$ denote their equilibrium status. By using (3.1) in the equilibrium equations of (2.7), (2.9), (2.11) and performing some rearrangements, we botain the following:

Now, we use (3.1) in (2.12) at equilibrium to obtain the following:

3.1.   Infection-free equilibrium and the basic reproduction number

Equation (3.5) implies that two cases arise: either $\lambda = 0$ or $\lambda \ne 0$. By (3.1), the first case implies that $I_1 = 0$ and $I_2 = 0$. Thus, by (3.2)–(3.4), we get $S_1 = (1 - { q_0 }) K, \; S_2 = { q_0 } K$ and $N = K$. Therefore, the model has an influenza-free equilibrium $E_0$, where

The basic reproduction number of the models (2.7)–(2.11) is computed by following the approach shown in ^[20]. To this end, we consider the Eqs (2.8) and (2.10) in computing the non-negative matrix for the new-infection term $\bf{T}$ and the non-singular matrix for the remaining transfer term ${\bf \Sigma}$ as follows

Consequently, the basic reproduction number $\mathcal{R}_0$ is the spectral radius of the matrix ${\bf T \: \Sigma^{-1} }$. Therefore,

It is noteworthy that the basic reproduction number $\mathcal{R}_0$ is not influenced by the immigration parameters $q_2, \Lambda,$ and $K$. However, it is affected by the entire population parameters $q_0, r, g, a, \beta, \gamma,$ and $\mu$. Definitely, $\mathcal{R}_0$ increases with the increase of the proportion of newborn females $q_0$ if and only if

where

● $D_m = \dfrac{1}{\gamma + \mu}$ is the length of males' infectious period, and

● $D_f = \dfrac{1}{a \gamma + \mu}$ is the length of females' infectious period.

Equation (3.9) says that the contagiousness and transmissibility of influenza A increases with the increase of the newborn females proportion if and only if the product of the relative susceptibility and the relative transmissibility (of females with respect to males) is higher than the ratio between the male and female infectious periods.

The local stability analysis of the infection-free equilibrium $E_0$ has been established based on linearization. The analysis revealed that $E_0$ is locally asymptotically stable if and only if $\mathcal{R}_0 < 1$ and the proof of this result is deferred to in Appendix 10.2.

Based on the above results, we show the following proposition.

Proposition 2. Models (2.7)–(2.11) have an infection-free equilibrium, given by $E_0 = \Big((1- q_0) K, 0, q_0 K, 0, K\Big)^{T}$, which is locally asymptotically stable if and only if $\mathcal{R}_0 < 1$, where $\mathcal{R}_0$ is the basic reproduction number.

3.2.   Endemic infection equation

If $\lambda \ne 0$, then Eq (3.5) implies that

On solving the algebraic system of Eqs (3.2)–(3.4) in terms of $\dfrac{S_1}{N}, \dfrac{S_2}{N}$ and $\dfrac{\Lambda}{N}$, we obtain the following:

where

and

A complete derivation of (3.11)–(3.15) is deferred to in Appendix 10.3. It is noteworthy that the formulas in (3.11) imply that

By using (3.16) into (3.1), we obtain the following:

Now, we substitute (3.11) into (3.10) to obtain the following:

By using (3.12)–(3.14) in (3.18) and rearranging the terms, we arrive at a second-degree polynomial equation (in $\lambda$) in the following form:

where

Equation (3.19) is the endemic force of the infection equation. It is quadratic in $\lambda$ and may have up to two feasible solutions. Here, the feasibility means that the solution values of $\lambda$ satisfy $\lambda \in [0, \infty)$. Once a solution of (3.19) is obtained, we substitute it in (3.13)–(3.15) and then in (3.16) and (3.17) to obtain the corresponding equilibrium components, and consequently the corresponding equilibrium point.

3.3.   Direction of bifurcation and the existence of multiple endemic equilibria

The polynomial in the left hand side of Eq (3.19) could be considered a function in the variable $\lambda$ and the parameter $\beta$, given that the other model parameters are kept fixed. Therefore, Eq (3.19) could be seen as a bifurcation equation, with $\beta$ being the bifurcation parameter. Hence, at $\lambda = 0$, there is a bifurcation point $(\beta_0, 0)$ in the plane $(\beta, \lambda)$, where

To investigate the direction of bifurcation at the bifurcation point $(\beta_0, 0)$, we make use of the implicit function theorem by following the same approach shown in ^[21]. Consequently, we compute (and study the sign of) the following expression:

where

Thus, the model exhibits the existence of a backward bifurcation if and only if $A_1{\Big\vert}_{ \beta = \beta_0} < 0$. The following two propositions, whose proofs are deterred to in Appendixes 10.4 and 10.5, summarize the conditions required for the existence of a backward bifurcation.

Proposition 3. Models (2.7)–(2.11) exhibit a backward bifurcation if and only if the following condition holds:

where

The following proposition presents a more specific equivalent set of inequalities to inequality (3.24).

Proposition 4. The inequality (3.24) is equivalent to either of the following two sets of inequalities (see Figure 2):

Figure 2 shows the regions in the plane $(g, r)$ for which a backward bifurcation occurs. The figure is drawn with parameter values as shown in Table 2. The figure shows that the backward bifurcation phenomenon (i.e., the existence of subcritical endemic states) does possibly exist if a high enough $r$ is chosen (i.e., $r > r_1$), while a small enough $g$ is chosen (i.e., $g < \ell_1$) and vice versa (i.e., $r < r_2 \& g > \ell_2$), given that the remaining parameters are kept fixed.

It is worth mentioning that if the inequality (3.24) does not hold, then models (2.7)–(2.11) show the existence of a forward bifurcation (i.e., supercritical endemic states), which, in sense, means that the model has a unique endemic equilibrium that exists if and only if $\beta > \beta_0$, or equivalently, if and only if $\mathcal{R}_0 > 1$.

3.4.   The effective threshold

By the effective contact rate threshold ^[22], we refer to a critical value of the successful contact rate $\beta$ below which the infection dies out without any effort. In case the model parameters are selected such that the model shows only a forward bifurcation, the effective threshold is the value of $\beta$ at which the basic reproduction number $\mathcal{R}_0 = 1$ (i.e., $\beta^\star = \beta_0$). However, if the backward bifurcation condition (3.24) holds, then two positive solutions of the quadratic equation (3.19) exist for values of $\mathcal{R}_0 < 1$. In this case, both solutions closely approach each other with the decrease of the value of $\beta$ until they coalesce at the turning point, see the schematic diagram in Figure 3.

The value of $\beta$ at the turning point, say $\beta_1$, is the effective threshold. To derive the formula of $\beta_1$, we proceed as follows. From the implicit function theorem, the conditions for the turning point are as follows:

The two conditions (3.28) are equivalent to the following:

Hence, the critical rate $\beta_1$ is the solution of the Eq (3.29) with respect to the contact rate $\beta$. By performing some computations and rearranging the terms, we obtain the following:

where

In summary, we show the following proposition.

Proposition 5. The effective contact rate threshold $\beta^\star$ is given by the following:

where $\beta_0$ and $\beta_1$ are defined in (3.20) and (3.30), respectively.

As a function of the relative transmissibility parameter $r$, the effective contact rate threshold $\beta^\star$ is depicted in Figure 4. The left Subfigure 4(a) corresponds to the backward bifurcation set of inequalities (3.26), while the right Subfigure 4(b) corresponds to the inequality set (3.27). In each subfigure, the solid curve represents $\beta = \beta_0$ (i.e., $\mathcal{R}_0 = 1$), while the broken curve represents $\beta = \beta_1$. In the region above the solid curve, the model has a unique endemic equilibrium. However, in the region that lies in between the solid and broken curves, the model has two endemic equilibria. Otherwise, the model has no endemic equilibrium.

In terms of the basic reproduction number $\mathcal{R}_0$, the effective basic reproduction threshold $\mathcal{R}_0^\star$ is the value below which the infection disappears and does not persist. In case the model shows a forward bifurcation, this effective threshold is $\mathcal{R}_0^\star = 1$. However, if the model exhibits a backward bifurcation, the effective threshold is as follows:

In summary, we have the following proposition.

Proposition 6. The effective basic reproduction number threshold is as follows:

The effective basic reproduction threshold $\mathcal{R}_0^\star$ is drawn in the plane $(r, \mathcal{R}_0)$ and shown in Figure 5. Figure 5(a) shows that, for small enough values of the relative susceptibility $g$ (of females with respect to males), multiple endemic equilibria do exist for values of $\mathcal{R}_0^1 \le \mathcal{R}_0 < 1$ in a very narrow region on the right of the plane $(r, \mathcal{R}_0)$. In this case, the multiple equilibria region becomes wide with the increase of the relative transmissibility parameter $r$. However, for high enough values of the relative susceptibility parameter $g$, the multiple equilibria region exists on the left of the plane $(r, \mathcal{R}_0)$ and diminishes with the increase of $r$, see Figure 5(b).

In case the model undergoes a backward bifurcation, it should be underlined that the effective contact rate threshold $\beta_1$ depends on the proportion of immigrating females $q_2$. Consequently, the effective basic reproduction number threshold $\mathcal{R}_0^1$ becomes influenced by any increase or decrease in $q_2$.

3.5.   The endemic equilibria: solutions and the bifurcation curve

Motivated by the results shown in Sections 3.3 and 3.4, Eq (3.19) has two feasible solutions if the condition (3.24) holds. The first solution is given by

and exists if and only if the condition (3.24) holds together with the inequality $\beta_1 \le \beta \le \beta_0$. It is worth mentioning that the inequality $\beta_1 \le \beta \le \beta_0$ is equivalent to $\mathcal{R}_0^1 \le \mathcal{R}_0 \le 1$, where $\mathcal{R}_0 = \beta/\beta_0$ and $\mathcal{R}_0^1 = \beta_1 / \beta_0$. The other solution of Eq (3.19) is given by

and exists if and only if $\beta > \beta_1$ (i.e., if and only if $\mathcal{R}_0 \ge \mathcal{R}_0^1$).

Both solutions are depicted in the plane $(\mathcal{R}_0, \lambda)$ and shown in Figure 6(a). The solution $\lambda^{-}$ is represented by the broken curve, while the solution $\lambda^{+}$ is represented by the solid curve. They collide at the turning point.

However, if the backward bifurcation condition (3.24) does not hold, then the Eq (3.19) has a unique feasible solution that definitely exists if and only if $\beta > \beta_0$ (or equivalently, $\mathcal{R}_0 > 1$). This unique solution is given by the formula (3.35) and is depicted in Subfigure 6(b).

By using (3.34) and (3.35) within (3.13)–(3.15) and then with (3.16) and (3.17), we obtain the corresponding equilibrium points whose formulas are determined in the following proposition.

Proposition 7. The equilibrium analysis of models (2.7)–(2.11) reveals the following results on the existence of endemic equilibria:

● If the condition (3.24) does not hold, then the model has a unique endemic equilibrium that exists if and only if $\mathcal{R}_0 > 1$ and is given by

where

● If the backward bifurcation condition (3.24) holds, then the model has two endemic equilibria. One of them is in the form (3.36) and exists if $\beta \ge \beta_1$ (i.e., $\mathcal{R}_0 \ge \mathcal{R}_0^1$). The other one exists if $\mathcal{R}_0^1 \le \mathcal{R}_0 < 1$ and is given by

where

● Otherwise, the model has no endemic equilibrium.

In case the model undergoes a backward bifurcation, it is worth noting that the solution $\lambda^-$ lies between $\lambda = 0$ and $\lambda = \lambda^+$. Moreover, the equilibrium point $E^-$ is unstable, while the equilibrium point $E^+$ is locally asymptotically stable, whenever it exists. Therefore, the curve representing the solution $\lambda = \lambda^-$ in the plane $(\mathcal{R}_0, \lambda)$ of Figure 6(a) is drawn as a broken line to distinguish it from the solid curve representing the solution $\lambda = \lambda^+$ that corresponds to the stable endemic equilibrium.

4.   Asymptotic stability of the endemic equilibria

Due to the complicated terms in the formulas of the endemic equilibria, simulations have been performed to study the asymptotic stability of the model's equilibrium points. To this end, the function ode45 in Matlab has been used to numerically solve the models (2.7)–(2.10) with various randomly selected initial conditions chosen so that they lie in the set of definition $\Omega$. Although the state variables in the model represent the number of individuals, we draw the solutions in the form of proportions to better present the results, as shown in Figures 7–9. These figures are produced with parameter values as shown in Table 2; the parameters $r$ and $g$ are given the values $r = 2$ and $g = 0.8$, while the contact rate $\beta$ is chosen so that the basic reproduction number $\mathcal{R}_0$ takes three different values to explain different scenarios for the fate of the trajectory solution (i.e., the attracting equilibrium point(s)). It is worth confirming that the state variable notations in the legend of the vertical axis in Figures 7–9 denote proportions rather than numbers.

Based on the chosen parameter values, the model undergoes a backward bifurcation at $\mathcal{R}_0 = 1$ and the turning point occurs at the point $(\mathcal{R}_0^1, \lambda^1) = (0.9973, 0.0163)$ in the plane $(\mathcal{R}_0, \lambda)$. Therefore, Figure 7 has been produced with a value of $\mathcal{R}_0 = 0.8 < \mathcal{R}_0^1$. In this case, no endemic equilibrium exists, and, in consequence, the figure shows that all trajectory solutions are attracted by the infection-free equilibrium ($I_1 = I_2 = 0$). However, the solutions are shown in Figure 8 with a value of $\mathcal{R}_0\; = \; 0.998\; \in\; (\mathcal{R}_0^1, 1)$. The figure shows that the solutions are attracted either by the infection-free equilibrium ($I_1 = I_2 = 0$) or by an endemic equilibrium (i.e., $I_1 > 0$ and $I_2 > 0$). Finally, the model has been solved with a value of $\mathcal{R}_0 = 1.3 > 1$ and the solutions are shown in Figure 9. The figure shows that all solutions are attracted by a unique non-trivial endemic equilibrium (i.e., $I_1 > 0$ and $I_2 > 0$).

Motivated by the results shown in Figures 7–9, our in silico simulations show three scenarios for the evolution of influenza infection if the model undergoes backward bifurcation at $\mathcal{R}_0 = 1$. The first scenario is that both of the time-dependent proportions of infected males and infected females eventually approach zero (i.e., their values at the influenza-free equilibrium as shown in Figure 7(b), (d)); therefore, the infection washes out without any further efforts. This scenario is ensured if the combination of the model parameters is chosen so that $\mathcal{R}_0 < \mathcal{R}_0^1$. Another scenario is that both of the time-dependent proportions of infected males and infected females eventually approach a positive level (i.e., their correspondents at the endemic equilibrium as shown in Figure 9(b), (d)); therefore, the infection persists in the population. This scenario happens if the combination of the model parameters used in the simulations satisfies $\mathcal{R}_0 > 1$. The third scenario mixes the above-mentioned two scenarios, where some of the solutions eventually approach the influenza-free equilibrium, while the others eventually approach an influenza endemic equilibrium (as shown in Figure 8(b), (d)). This scenario occurs if the combination of the model parameters used to simulate the model are chosen to satisfy $\mathcal{R}_0^1 \le \mathcal{R}_0 \le 1$. In the third scenario, the infection elimination depends on the initial conditions. It is worth noting that in case the model only undergoes a forward bifurcation, then there are only the first two scenarios.

5.   The endemic prevalence of infection

In the absence of an endemic infection (i.e., $\lambda = 0$), the formulas (3.12)–(3.15) imply that

Therefore, by (3.11) and (3.17), we obtain the subpopulation proportions

and the total population size $N = K$.

However, in the presence of an endemic infection (i.e., $\lambda \ne 0$), the formulas (3.11)–(3.15) and (3.17) help compute the following expressions.

● The total population size at equilibrium in the endemic situation is given by the following:

● The male equilibrium proportion in the endemic situation, say $p_m$, is given by the following:

● The female equilibrium proportion in the endemic situation, say $p_f$, is given by the following:

● The endemic prevalence of infection in the male population, say $p_{I_1}$, is given by the following:

● The endemic prevalence of infection in the female population, say $p_{I_2}$, is given by the following:

● The endemic prevalence of infection in the overall population, say $p_{I}$, is given by the following:

It is worth mentioning that $\lambda$ is the males' equilibrium force of infection and is(are) the feasible solution(s) of the Eq (3.19).

Based on various values of the relative transmissibility parameter $r$ and the relative susceptibility parameter $g$ selected from the different regions in the plane $(g, r)$, the expressions in the formulas (5.3)–(5.8) have been drawn as functions of the basic reproduction number $\mathcal{R}_0$ and presented in Figures 10–13, while keeping the other parameters fixed as shown in Table 2. The values of the pair $(g, r)$ are shown in the head of each subfigure.

Figures 10 and 11 have been produced with values of $r$ and $g$ such that the model undergoes a forward bifurcation at $\mathcal{R}_0 = 1$. The curves in both figures correspond to the solution $\lambda = \lambda^+$ (i.e., they are connected to the endemic equilibrium $E^{+}$). Each figure is produced with three different values of the pair $(g, r)$, where $r$ is kept fixed for the same figure, though $g$ is allowed to change. The solid curves are drawn with higher values of the relative susceptibility ($g = 1.5 > \ell_2 = 1.2932$) than in the cases of the dotted curve ($g = 0.8 < \ell_1 = 1.1888$) and the broken curve ($g = 1.10 < \ell_1 = 1.1888$).

From a demographic perspective, it is clear that the Figure 10(a), (d) shows that higher values of the relative susceptibility level $g$ implies a reduction in the equilibrium total population size and the equilibrium proportion of female subpopulation, while the contrast is remarkable for the proportion of males at equilibrium $p_m$ as shown in the Figure 10(c). Figure 10(c), (d) shows that, in the case of a high enough relative susceptibility level, the proportion of males $p_m$ (females $p_f$) at equilibrium strictly initially increases (decreases) with an increase of $\mathcal{R}_0$ till reaching a maximum (minimum) and then decreases (increases) with an increase of $\mathcal{R}_0$. However, contrasting qualitative behaviors of both $p_m$ and $p_f$ are remarkable in the case of low enough relative susceptibility levels ($g = 0.8 < \ell_1 = 1.1888$ and $g = 1.1 < \ell_1 = 1.1888$).

From an epidemiological perspective, Figure 10(b), (e), (f) shows that the high relative susceptibility of females with respect to males increases the endemic prevalence of infection in the overall population $p_I$, the endemic prevalence of infection within males $p_{I_1}$, and the endemic prevalence of infection within females $p_{I_2}$, respectively. These prevalences strictly increase with an increase of the basic reproduction number $\mathcal{R}_0$.

It is worth explaining that even if the relative transmissibility parameter $r$ is increased (but the model still undergoes only forward bifurcation at $\mathcal{R}_0 = 1$), then the qualitative behavior doesn't differ from that shown in Figure 10, see Figure 11.

The above-described qualitative behavior is ensured as long as the males case fatality is higher than that of the females (i.e., if $c_1 > c_2$). However, in the opposite case (i.e., if $c_1 < c_2$), this behavior is different as shown in Figure 14, where two contrasting behaviors are remarkable for values of $\mathcal{R}_0 > 1$. Definitely, for values of $\mathcal{R}_0$ in the right-neighbourhood of $\mathcal{R}_0 = 1$, the equilibrium total population size $N$ and the equilibrium proportion of females $p_f$ decrease with an increase of the relative susceptibility $g$. However, for high enough values of $\mathcal{R}_0$, they increase with an increase of the relative susceptibility parameter $g$. On the contrary, the endemic prevalence levels $p_I, p_{I_1},$ and $p_{I_2}$ decrease with an increase of the relative susceptibility parameter $g$ for values of $\mathcal{R}_0$ slightly above one, while they increase for high enough values of $\mathcal{R}_0$. Based on values of $r$ and $g$ that ensure the model undergoes a backward bifurcation at $\mathcal{R}_0 = 1$, the expressions in (5.3)–(5.8) are drawn as functions of the basic reproduction number $\mathcal{R}_0$ and presented in Figures 12 and 13. There are two curves for values of $\mathcal{R}_0$ slightly less than one. The broken curve corresponds to the unstable endemic equilibrium (computed based on a value of $\lambda = \lambda^-$), while the solid one corresponds to the stable endemic equilibrium (computed based on a value of $\lambda = \lambda^+$), as implicated by the occurrence of a backward bifurcation. Only the solid curves are of interest, as they represent the stable endemic equilibrium computed based on a value of the endemic force of the infection solution $\lambda^+$. Both figures show that the qualitative behavior of the functions based on the stable endemic equilibrium remains the same as the behavior shown in Figure 10 and is detailed in the above description. It is worth mentioning that several simulations have been performed to explore the qualitative behavior if the females' case fatality is higher than the males' one. The simulations show that the qualitative behavior of the demographic expressions $N, p_m,$ and $p_f$ is similar to their correspondences in Figure 14, while the endemic prevalences $p_I, p_{I_1},$ and $p_{I_2}$ keep the same behavior as in the case $c_1 > c_2$.

Motivated by the above analysis, we come up with the following.

1) If the male's case fatality is higher than the female's one, then we have the following results:

● Reducing the relative susceptibility of females with respect to males reduces the endemic prevalence of the infection in the total population and in each sex-structured subpopulation, though it increases for both the female proportion $p_f$ and the total population size at equilibrium.

● For high enough levels of the relative susceptibility ($g > \ell_2$), reducing the basic reproduction number $\mathcal{R}_0$ reduces the equilibrium proportion of females until a minimum close to the right-neighbourhood of $\mathcal{R}_0 = 1$ is reached, and then increases again, while the converse is true for small levels of the relative susceptibility parameter $g < \ell_1$.

2) However, if the male's case fatality is lower than the female's one, then the two contrasting scenarios are remarkable. For small values of $\mathcal{R}_0$ which lie in the right neighbourhood of $\mathcal{R}_0 = 1$, the above-mentioned implications are reversed, while those implications do still work for higher values for $\mathcal{R}_0$.

6.   Summary and conclusions

The dynamics of influenza A infection has been extensively explored based on mathematical models. However, mathematical models that take the inequalities due to differences in sex and gender into account are less investigated. Motivated by a previous work ^[4], a SIS model was extended and adapted to describe the dynamics of influenza A in an open population with varying size, where the infection lethality was taken into account. From an epidemiological perspective, the model was adapted to consider differences in susceptibility, infectivity, infection-induced mortality (i.e., lethality of the infection), and recoverability between males and females. However, from a demographic point of view, the model was adapted to consider an open population with a population-size-structured immigration rate. Additionally, the inequality of the birth ratio of females and males was included.

The model has been mathematically analyzed. Definitely, the well-posedness of the model was shown, where the existence and uniqueness of time-dependent solutions and the positive invariance of the model's definition set was proven. The model's equilibrium analysis revealed that the model has an influenza-free equilibrium that was proven to be locally asymptotically stable if and only if $\mathcal{R}_0 < 1$, where $\mathcal{R}_0$ is the model's basic reproduction number. Moreover, the bifurcation analysis showed that the model underwent a backward bifurcation at $\mathcal{R}_0 = 1$ for a certain space-set of the model parameters. The conditions for the occurrence of a backward bifurcation was determined and presented in the form of either of the two inequality sets ($r > r_1, g < \ell_1$) or ($r < r_2, g > \ell_2$), where $r_1, r_2, \ell_1$ and $\ell_2$ were defined within the text (formulas (3.25)–(3.27)). The occurrence of a backward bifurcation made the model's behavior more complicated than in case of a forward bifurcation, especially when discussing the possibility to eliminate the infection.

The mathematical implication of the backward bifurcation phenomenon in epidemic models is that a two influenza-endemic (i.e., with positive levels of the infection's state variables) equilibria co-exists with the influenza-free equilibrium for values of $\mathcal{R}_0^1 \le \mathcal{R}_0 < 1$, where the endemic equilibrium with a higher infection level (i.e., with $\lambda = \lambda^+$) is locally asymptotically stable, while the endemic equilibrium with a lower level of the infection (i.e., with $\lambda = \lambda^-$) is unstable. The asymptotic local stability of the model's equilibrium solutions was numerically investigated.

The epidemiological implication of the existence of a backward bifurcation is that reducing the basic reproduction number $\mathcal{R}_0$ to values slightly less than one is a necessary but no longer a sufficient condition to eliminate the infection. In other words, in the case of a backward bifurcation, the value $\mathcal{R}_0 = 1$ is no longer a threshold, while $\mathcal{R}_0 = \mathcal{R}_0^1$ is the effective threshold (whose formula is given in (3.32)), and strategies aiming to eliminate the infection would be based on reducing $\mathcal{R}_0$ to slightly below the effective threshold value $\mathcal{R}_0^1$. Therefore, the minimum effort required to eliminate the infection becomes increased ^[21,23].

Our analysis showed that the proportion of female new-immigrants didn't affect the value of the basic reproduction number, but it affected the value of the effective basic reproduction number threshold and its influence was associated with the value of the relative susceptibility parameter $g$. However, the basic reproduction number was affected by the value of the female newborns proportion $q_0$ and the impact's type was determined by the inequality (3.9).

Some demographic quantities (precisely, the total population size and proportion of male and female subpopulations at equilibrium) and some epidemiological expressions (precisely, the endemic prevalence of the infection in the total population, as well as in both male and female subpopulations) were computed and numerically investigated. The effect of changes in the relative transmissibility and susceptibility parameters ($r$ and $g$, respectively), as well as in the infection's case fatality in both the male and female subpopulations on these demographical and epidemiological expressions, were numerically investigated.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors declare that this research is not supported by funding. And, the authors are grateful to the anonymous reviewers for their constructive comments.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.1.   Proof of Proposition 1

To show that the model is well-posed, we first consider the $N$-equation in (2.11), where

Hence, by the comparison theorem, we obtain the following:

where $N(0)$ is the total population size at $t = 0$. Additionally, we have

and using the comparison theorem, we obtain the following:

Hence,

i.e., $N(t)$ is upper-bounded.

In the same manner, we may use (2.7)–(2.11) to write

and thereby we get

Now, since $0 \le S_1(t) + I_1 (t) + S_2(t) + I_2 (t) = N(t)$, and since both $N(t)$ is upper-bounded, then $S_1(t), I_1 (t), S_2(t), I_2 (t)$ are all upper-bounded and non-negative for all $t \ge 0$. Thus, the set $\Omega$ is positively invariant.

To show the uniqueness of the time-dependent solutions of models (2.7)–(2.11), we note that the right hand side of these equations are all continuous in the models state variables $S_1, S_2, I_1,$ and $I_2$. Moreover, it is easy to check that all partial derivatives of the functions in the right-hand side of models (2.7)–(2.11) are also continuous in the model state variables. Therefore, they are locally Lipschitz, and, hence, any time-depending solution starting with initial conditions in $\Omega$ is unique.

A.2.   Proof of the local stability of the influenza-free equilibrium $E_0$

We apply the linearization principle to establish the local stability analysis of the influenza-free equilibrium $E_0$. Rather than traditionally using the $S_1, I_1, S_2,$ and $I_2$ equations in models (2.7)–(2.11), we consider the $N_1, N_2, I_1,$ and $I_2$ equations, as $S_1 = N_1 - I_1$ and $S_2 = N_2 - I_2$. Therefore, we consider the following model:

which has the following influenza-free equilibrium $\widehat{E}_0 = \Big((1- q_0) K, q_0 K, 0, 0 \Big)^{T}$ that corresponds to $E_0$. The local stability of $\widehat{E}_0$ implies the local stability of $E_0$. The Jacobean matrix evaluated at $\widehat{E}_0$ is a block matrix in the following triangular form:

where

It is noteworthy that the characteristic polynomial of the block triangular matrix (A2.1) is the product of the characteristic polynomials of the matrices $\mathbb{A}$ and $\mathbb{D}$. Therefore, the local stability of the influenza-free equilibrium $\widehat{E}_0$ is ensured if and only if the following conditions on the matrices $\mathbb{A}$ and $\mathbb{D}$ hold:

To this end, we compute the following:

It is clear that $\det (\mathbb{D}) > 0$ if and only if $\mathcal{R}_0 < 1$. However, if $\mathcal{R}_0 < 1$, then we have $(1 - q_0) \beta < (\gamma + \mu)$ and $q_0 r g \beta < (\mu + a \gamma)$, which implies that $\text{Tr} (\mathbb{D}) < 0$. Thus, the condition (A2.2) holds if and only if $\mathcal{R}_0 < 1$.

A.3.   Derivation of the persistence-infection polynomial equation (3.19)

To derive Eq (3.19), we rewrite the algebraic system of equations (3.2)–(3.4) in the following simple form:

where

Now, we solve the systems (A3.1)–(A3.3) in terms of $\dfrac{S_1}{N}, \dfrac{S_2}{N}$, and $\dfrac{\Lambda}{N}$ to obtain the following:

where

where

A.4.   Proof of Proposition 3

To prove Proposition 3, we make use of the formulas (3.21)–(3.23). They imply that our model exhibits a backward bifurcation if and only if $A_1 \Big\vert_{\beta = \beta_0} < 0$. Now, we have the following:

Hence, $A_1 \Big\vert_{\beta = \beta_0} < 0$ if and only if

i.e.,

i.e.,

i.e.,

where

A.5.   Proof of Proposition 4

To prove Proposition 4, we notice that the numerator of $\ell_1$ is less than the numerator of $\ell_2$, while the denominator of $\ell_1$ is bigger than the denominator of $\ell_2$. Therefore, $\ell_2 > \ell_1$. Hence, $(\ell_1 - g)$ and $(\ell_2 - g)$ have the same sign (either negative or positive) if and only if either $g < \min(\ell_1, \ell_2) = \ell_1$ or $g > \max(\ell_1, \ell_2) = \ell_2$. Therefore, the inequality (A4.1) holds if either $g < \ell_1$ or $g > \ell_2$. In the first case, we obtain the following condition:

while in the second case, we obtain the following condition:

This completes the proof.