An inherently discrete–time $SIS$ model based on the mass action law for a heterogeneous population

1.   Introduction

In many diseases, one can indicate the group in the population with a higher susceptibility to infection compared to the other groups. For that reason, while constructing an epidemic model, it is natural to assume that the population consists of two subpopulations: with a low risk ($LS$) and a high risk ($HS$) of getting infected. We called such population a heterogeneous one – in other words, the heterogeneous population consists of two homogeneous populations. Moreover, while proposing an epidemiological model, we must indicate stages through which an individual passes. The most popular class of models is the $SIS$ (susceptible–infected–susceptible) class – in which a recovered individual does not gain immunity and can become infected again. The author in ^[1] described and investigated the exemplary $SIS$ model. $SIS$ models can be extended while including, for example, the $E$ (exposed – infected but firstly not infectious) or $R$ (recovered) classes. The analyses of the obtained sample $SEIS$, $SIR$, and $SEIRS$ models are presented in ^[2,3,4]. However, for many diseases, a lack of data often prevents us from doing such an extension.

A crucial thing in constructing the model is assuming whether biological processes are either continuous or discrete. Discrete–time epidemiological models are less analyzed than the continuous–time models. The general approach in constructing discrete–time systems is the discretization of continuous–time counterparts. Discretization methods range from simple, including the explicit Euler method ^[5], to complicated ones such as non–standard discretization schemes ^[6], which is an idea explained in ^[7]. However, applying any discretization method results in the occurrence of a step size of the discretization method, thus being an additional parameter in a system. The step size has no clear biological meaning. For that reason, it is better to construct a discrete–time epidemiological model on the grounds of its biological description without creating its continuous–time counterpart, followed by its discretization. The system obtained with this method is called an inherently discrete–time method. Although this approach is rare, there are several papers related to this issue. In ^[8], the author proposed discrete–time $SIS$ and $SEIR$ models and provided their analysis, thereby focusing on construction of a disease transmission function and a local stability analysis of the stationary states that appeared in the models. Paper ^[9] deals with an investigation of chaotic behavior, including bifurcation, in the $SIS$ model. In ^[10], the author studies an effect of seasonal diseases and trends on behavior $SIS$ of the $SIS$ model. An interesting approach to discrete–time epidemic modeling can be found in the recent paper ^[11]. The author constructed epidemic discrete–time models with two time scales. He considered two cases: when dynamics of the epidemic was either slower or faster compared to the dynamics of the other processes. Later, the author reduced the proposed systems to $SIS$ models and conducted a local stability analysis.

The mentioned papers concern the case of a homogeneous population. Our paper aims to construct and analyze a discrete–time $SIS$ model for a heterogeneous population without the discretization of its continuous–time counterpart. The case of population heterogeneity makes a model analysis complicated. However, the heterogeneity is crucial for the proper investigation of epidemic dynamics, which was previously proven medically ^[12]. What is important is that we assume that the values of every parameter are different for each subpopulation; this assumption makes the population purely heterogeneous.

Our paper continues our previous work ^[13], where we proposed a similar model. However, in our previous work, we proposed the functions reflect the probability of staying healthy of type $G \left(\beta \frac{I(t)}{N(t)} \right)$, where $\beta$ is a transmission coefficient, $I(t)$ is a number of infected people from a given subpopulation at time $t$, and $N(t)$ is a number of the whole given subpopulation. This type of function was proposed in ^[8]. In this paper, we introduced a simple form of function $G \Big(\beta I(t) \Big)$, which, in some sense, relied on the mass action law. Here and in ^[13], we assume that the disease is spread among each subpopulation separately and from $HS$ to $LS$. In many cases, people from $HS$ tend to infect themselves in their community, for example, highly sexually active people with sexually transmitted diseases or people living in poor conditions in the context of airborne infections. For that reason, we reasonably neglect the pathogen spread from $LS$ to $HS$. The components of models in ^[13] and in this paper, which describe particular biological phenomena, are analogical to those from our previous paper ^[14], where we proposed the continuous–time system and conducted a subsequent analysis.

We complement our results with a numerical simulation that concerns the case of a tuberculosis ($TB$) spread in the Warmian–Masurian Province in Poland between 2000–$2023$. As $LS$ and $HS$, we consider non–homeless and homeless people, respectively. The choice of the Warmian–Masurian Province is motivated by the implementation of Active Case Finding (ACF) programs among homeless people therein. These programs caused a decrease in the incidences of $TB$ in both the homeless and non–homeless communities in the region.

This paper is organized as follows. The next section describes the assumption of the model and presents its form. In Section 3, we investigate the existence of stationary states of the systems. The next two sections concern the local and global stability analyses of the states. Section 6 contains computations which yield the basic reproduction of our model. Then, we conduct a numerical simulation of the $TB$ case mentioned earlier. We conclude our results in Section 8.

2.   Mathematical model and its description

Now, we introduce the components of our model. If a variable or a parameter has a lower subscript equal to $1$ or $2$, it refers to $LS$ or $HS$, respectively. For the lower subscript expressed by $i$, we have $i \in \{1, 2 \}$. The variables $S_1$ and $S_2$ represent the size of $LS$ and $HS$, respectively. The number of infected people are denoted by $I_1$ and $I_2$, respectively. By $N_i : = S_i + I_i$, we clearly define the size of the $i$–th subpopulation. The new individuals move into healthy groups in each population with an inflow $C_i$. By $\gamma_i$ and $\alpha_i$, we denote the rates of recovery and disease–related deaths, respectively. Every individual survives with the probability $r_i$.

A healthy individual from class $S_i$ can move to class $I_i$ because of contact with an infected person. The efficiency of the diseases transmission in the sole $i$–th subpopulation is expressed by $\beta_i$ and $\beta$ reflects efficiency of the transmission from $HS$ to $LS$. As we stated in the introduction, we omit the case of the pathogen transmission from the $I_1$ to $S_2$ classes. Because of this omission, we must introduce two different functions which reflect the probability of remaining healthy in both LS and HS. We denote these functions by $G$ and $H$. We assume that $G(x)$ has the following properties:

Analogically, for $H$, we have the following:

After recovery, the individuals from class $I_i$ return to class $S_i$.

Our proposed model reads as follows:

where $S_n^{(i)}$ and $I_n^{(i)}$ are the values of variables at the $n$–th node of the discrete time scale. The parameters are fixed and positive. Because of their meaning, we assume that $r_i$, $\gamma_i$, $\beta_1$, $\beta \in (0, 1)$, and $\alpha_i \in (0, r_i)$.

Thanks to the inherent discretization, there is no need to include a step size of the discretization method in the system. Moreover, System (2.1) has no negative terms that occur during the discretization of a continuous–time one (cf. ^[5]).

In order to make the form of System (2.1) more transparent, we define the following:

and we rewrite System (2.1) as follows:

For the sake of simplification, if it does not yield ambiguity, we will use the following notation:

Observe that for a initial condition

we have $S_n^{(1)}, S_n^{(2)} > 0$, and $I_n^{(1)}, I_n^{(2)} \geq 0$ for any $n \in \mathbb{N}$.

3.   Existence of stationary states

Let us indicate the stationary states which appear in System (2.2). For every stationary state, the system reads as follows:

Adding either Eqs (3.1a) and (3.1b) or Eqs (3.1c) and (3.1d) yields the following:

giving

From Eq (3.2), we obtain the following:

where $\sigma_i : = 1-r_i + \alpha_i$. Clearly, $0 < \sigma_i < 2$. Taking $I_i = 0$ leads to the following disease–free stationary state:

This state always exists. If $I_i \neq 0$, then from (3.4), we obtain positivity of $S_1$ for

and we have positivity of $I_i > 0$ for

Observe that $\sigma_i$ is separated from $0$, hence $\frac{C_i}{\sigma_i}$ is finite. Moreover, thanks to Eq (3.5), we can restrict the domains of the functions $G$ and $H$ to accordingly

Additionally, we can assume the following:

Therefore, $G$ and $H$ are surjections. Because $G$ and $H$ are injections, we conclude that $G$ and $H$ are invertible.

For the sake of simplification, we will use the following parameter in the latter part of this work:

Clearly, $\kappa_i \in (0, 1).$

Observe that because of the lack of disease transmission from $LS$ to $HS$, the dynamics of growth of $HS$ is independent of the dynamics of growth of $LS$. Hence Eqs (2.2c) and (2.2d) constitute an autonomous system and can be solely investigated.

Let us analyse a case when $I_2 > 0$. We obtain the following theorem:

Theorem 1. System (2.2c) and (2.2d) has a positive stationary state $(S_2, I_2) = (\overline{S}_2, \overline{I}_2)$, where $\overline{S}_2$ reads

and $\overline{I}_2$ is a solution of an equation

under a condition

This state exists if

Proof. The form of $\overline{S}_2$ in Eq (3.9) results from Eq (3.4).

Now we substitute $S_2 = N_2 - I_2$ into (3.1d) and obtain the following:

From Eqs (3.4) and (3.5), we have

and Eq (3.11), accordingly. Considering Eq (3.14) in Eq (3.13) gives the following:

what can be transformed to Eq (3.10).

Now, let us analyze the left–hand side of Eq (3.10) as the following function:

This function is continuous in $\Big[0, \frac{C_2}{\sigma_2} \Big)$. We have $F(0) = 0$ and

Moreover, for $I_2 \in \Big[0, \frac{C_2}{\sigma_2} \Big)$, we have

and

Additionally, we define the following auxiliary function:

and investigate the intersection point of functions $F(I_2)$ and $\tilde{F}(I_2)$ for $I_2 \in \Big[0, \frac{C_2}{\sigma_2} \Big)$. $H$ is decreasing and convex, so $\tilde{F}$ is increasing and concave. Moreover, $\tilde{F}(0) = 0$. Clearly, $F$ and $\tilde{F}$ intersect at $0$. Furthermore, the properties of these functions yield the second unique intersection point $\overline{I}_2 > 0$ if and only if

From Eq (3.16), we have the following:

From Eq (3.17), we obtain the following:

Considering Eqs (3.19) and (3.20) in Eq (3.18) gives Eq (3.12).

Now, we focus on the case $I_2 = 0$, giving $S_2 = \tfrac{C_2}{1-r_2}$. Under these equalities, we expect to obtain a stationary state of (2.2) with $S_1, I_1 > 0$. If we take $(S_2, I_2) = \left(\tfrac{C_2}{1-r_2}, 0 \right)$ in Eqs (3.1a) and (3.1b), then the reasoning for the pair $(S_1, I_1)$ is analogical to this from the proof of Theorem 1. We obtain an equation, analogical to Eq (3.10), that reads as follows:

It has one positive unique solution $I_1 = \widehat{I}_1$. We formulate the following corollary analogical to Theorem 1:

Corollary 1. In System (2.2), there exists a stationary state

with

and $0 < \widehat{I}_1 < \tfrac{C_1}{\sigma_1}$ being a solution of Eq (3.21). This state exists if

Now, let us investigate the existence of a postulated positive stationary state.

Theorem 2. In System (2.2), there is a positive (endemic) stationary state $E_e : = \left(\overline{S}_1, \overline{I}_1, \overline{S}_2, \overline{I}_2 \right)$, where $\overline{I}_1$ is a solution of

held for

This state exists if (3.12) and

Proof. It is clear the pair $(\overline{S}_2, \overline{I}_2)$, being the stationary state of System (2.2c) and (2.2d), defined in Theorem 1, constitute the values of coordinates $(S_2, I_2)$ in the whole System (2.2). The proof of that theorem yields condition (3.12).

We repeat the approach from that proof. Considering $(\overline{S}_2, \overline{I}_2)$ in Eq (3.1b) yields Eq (3.24) under condition (3.25). We introduce the following functions:

We have $F', F'', \tilde{F}' > 0$, and $\tilde{F}'' < 0$. There is a unique positive intersection point $\overline{I_1}$ of the functions $F(I_1)$ and $\tilde{F}(I_1)$ if and only $F'(0) \leq \tilde{F}' \left(\beta \overline{I}_2 \right)$, which can be written as (3.26).

Let us focus on coordinates $\widehat{I}_1$ and $\overline{I}_1$ from states $E_1$ and $E_e$. Accordingly, they are positive unique solutions of (3.21) and (3.24). From the properties of function $G$, we have the following:

Hence, we get that $\widehat{I_1} < \overline{I_1}$.

Now, we compare Eqs (3.23) and (3.26) that appear in Corollary (1) and Theorem (2), which state the existence of $E_1$ and $E_e$, respectively. Because of the properties of function $G$, we have $G'(0) < G'\left(\beta \overline{I}_2 \right) < 0$. Clearly, Equation (3.26) is stricter than Eq (3.23). Thus, without considering Eq (3.12), we state that the $E_e$ existence occurs for smaller ranges of the parameter values in comparison to the existence of $E_1$. This fact is naturally eligible from the medical point.

4.   Stability of the stationary states

Now, let us investigate the local stability of the obtained stationary states. The Jacobian matrix for System (2.2) reads as follows:

This matrix can be written as a block matrix $J = \begin{pmatrix} J_1 & J^* \\ 0 & J_2 \end{pmatrix}$, where

In order to investigate the stability from the eigenvalues of $J$, it is sufficient to consider the eigenvalues of the matrices $J_1$ and $J_2$.

4.1.   Local stability of the disease–free state

We start from the local stability analysis for $E_{df}$.

Theorem 3. The state $E_{df}$ of System (2.2) is locally stable if

and

Proof. Since $G(0) = 1$ and $H(0) = 1$, we write the following:

From both matrices, we obtain the following eigenvalues:

Because of the meaning of $r_i$, we get $|\lambda_{1, 3}| < 1$. Condition $|\lambda_2| < 1$ can be expressed as follows

This inequality can be split into two separates ones:

and

The meaning of $\kappa_1$ yields negativity of the left–hand side of Eq (4.3). Clearly, its right–hand size is always positive. Therefore, Equation (4.3) is always true. The fulfillment of $|\lambda_2| < 1$ is determined by Eq (4.4), which can be expressed as Eq (4.1). For condition $|\lambda_4| < 1$, we conduct analogical reasoning as for $|\lambda_2| < 1$ to obtain Eq (4.2).

For Eq (3.23), we determine the condition for the $E_1$ existence. It stays on the contrary to Eq (4.1), which is one of the conditions for the $E_{df}$ stability. Hence, we conclude the following:

Corollary 2. If the state $E_1$ exists, the state $E_{df}$ loses its stability.

4.2.   Local stability of the state with the infection

Now, we focus on the local stability of $E_1$. In the next theorem and its proof, we will use the following notations:

We are interested in the long–time dynamics of epidemics, thus $\eta_1 > 0$. In the proof, $I$ represents an identity matrix.

Theorem 4. The existing state $E_1$ of System (2.2) is locally stable if (4.2),

and

Proof. The proof is conducted under the $E_1$ existence. Observe that for this state, we have $J_2 (E_1) = J_2 (E_{df})$, which yields condition (4.2). If we add the first row to the second one of the determinant $|J_1(E_1) - \lambda I|$, we obtain the following:

The characteristic polynomial of $J_1(E_1)$ reads as follows:

The discriminator of $P(\lambda)$ equals to the following:

If (4.5), then $\Delta > 0$ (we neglect the non-generic case $\Delta = 0$) and the polynomial's eigenvalues are always real and can be written as follows:

We investigate inequalities $\lambda_1 > -1$ and $\lambda_2 < 1$. The first one holds if

For $2-b > 0$, which can be written as (4.6), we raise both sides of Eq (4.8) to a square and eventually get $b < 1 + c$. This inequality can be transformed to Eq (4.7). The condition $\lambda_2 < 1$ provides $\sqrt{b^2-4c} < 2 + b$, which is a weaker condition than Eq (4.8).

Now, we determine the stability of $E_e$. We use the following notations:

The determinant $|J_1(E_e) - \lambda I|$ reads as follows:

what is similar to the determinant $|J_1(\widehat{E}_1) - \lambda I|$. Additionally, we have the following:

The analysis of this determinant is analogical to the analysis of $|J_1(\widehat{E}_1) - \lambda I|$. Hence, we can perform a reasoning analogical to this for the local stability of $E_1$. We conclude the following:

Corollary 3. The existing state $E_e$ of System (2.2) is locally stable if $\mathcal{G} < \gamma_1$, $\mathcal{H} < \gamma_2$,

and

4.3.   Case of modified functions $G$ and $H$

In the analysis of the stability of the stationary states of System (2.2), we refer to $G'(0)$ and $H'(0)$. Let us introduce the explicit forms of functions $G$ and $H$ to obtain the explicit results. We aim to provide straightforward results; hence, we want to apply as simple of functions as possible. For this reason, we treat $G$ and $H$ as linear functions defined on domains from (3.7). The linear character of the proposed functions implies the modification of the particular conditions for functions $G$ and $H$. Naturally, we must have $G''(x) = H''(x) = 0$. Moreover, we assume that $G(0) = H(0) = 1$ and (3.8).

Let us compute the value of slope $a_G$ of $G$. We obtain the following:

Analogically, we get that the value of the slope $a_H$ of $H$ is as follows:

Clearly, we have the following:

We formulate the lemma concerning the local stability of $E_{df}$ for System (2.2) with a linear $G$ and $H$. The conditions for the local stability of the other stationary states can be analogically obtained.

Lemma 1. Assume that (4.11). $E_{df}$ is locally stable if

and

Proof. Considering Eq (4.11) in Eqs (4.1) and (4.2), we accordingly obtain the following:

From (4.14), we get the following:

When we introduce the definition of $\kappa_1$, this can be transformed to the following:

Observe that it is enough to check the following two inequalities:

and

Equation (4.16) can be simplified to Eq (4.12). Equation (4.17) yields always true to the inequality $0 < \gamma_1 \sigma_1 + \alpha_1 \gamma_1$. Analogically, we obtain Eq (4.12).

5.   Global stability of stationary states

In this section, we determine conditions for the global stability of the stationary states.

5.1.   Global stability of the disease–free state

First, we focus on $E_{df}$. We formulate the following theorem:

Theorem 5. The state $E_{df}$ of System (2.2) is globally stable if (4.2) and

Proof. Substituting Eqs (3.3) and (3.4) in Eq (2.2d) gives the following:

Let us define a function $\Theta: \Big[0, \frac{C_2}{\sigma_2} \Big) \to \Big[0, 1 + \frac{C_2}{\sigma_2} \Big)$ such that

The function $\Theta (I_2)$ can be treated as a reproduction function for the infected individuals from the high–risk group. The set of iterates of $\Theta$ is equivalent to the set of the sequence, which is generated by Eq (2.2d). After differentiation $\Theta$ with respect to $I_2$, we obtain

and

Observe that

which can be simplified to

From Eq (4.2), we have the following:

Combining Eqs (5.4) and (5.5), we obtain $\Theta' (0) < 1$. Hence, a fixed point $I_2 = 0$ is locally stable under the $\Theta-$iteration. Observe that the signs of $\Theta'(I_2)$ and $\Theta''(I_2)$ do not explicitly depend on $I_2 \in \Big[0, \frac{C_2}{\sigma_2} \Big)$. Remind that $H''(\beta_2 I_2) > 0$. Hence, $\Theta''(I_2) < 0$. Combining it with the inequality $\Theta' (0) < 1$ gives $\Theta' (I_2) < 1$ and $\Theta (I_2) < I_2$. Therefore, a sequence $\left(I_n^{(2)} \right) _{n = 0}^{\infty}$ is strictly decreasing and bounded below by zero. In the interval $\Big[0, \frac{C_2}{\sigma_2} \Big)$, this sequence converges to the only fixed point, that is $I_2 = 0$. From Eq (3.4), we obtain the following:

Now, let us substitute Eqs (3.3) and (3.4) in Eq (2.2b) to obtain the following:

Because of its properties and Eq (3.5), function G can be assessed in a way

giving

Instead of Eq (5.6), we will consider the following inequality:

For the sake of simplification, we replace (5.7) with the following analogical equality:

where we introduce a new variable $Y_1$. Observe that the analysis of such an obtained equation is similar to the analysis of Eq (5.2). We obtain Eq (5.1) and get that the sequence $\left(Y_n^{(1)} \right) _{n = 0}^{\infty}$ is therefore strictly decreasing and bounded below by zero. Hence, $\lim\limits_{n \to \infty} Y_n^{(1)} = 0.$ Our aim is to analyze Eq (5.7), hence we conclude that the sequence $\left(I_n^{(1)} \right) _{n = 0}^{\infty} \leq \left(Y_n^{(1)} \right) _{n = 0}^{\infty}$ is also therefore strictly decreasing and, because of the definition of $I_1$, bounded below by zero. Naturally, we get $\lim\limits_{n \to \infty} I_n^{(1)} = 0$ and $\lim\limits_{n \to \infty} S_n^{(1)} = \frac{C_1}{1-r_1}$. To sum up, we obtain the following

which yields the global stability of $E_{df}$.

Let us look at Eqs (4.1) and (5.1) which appear in Theorems 3 and 5, respectively. From the properties of function $G$, we obtain that Eq (5.1) is stronger than Eq (4.1), which is expected.

5.2.   Global stability of the states with the infection

Now, we discuss the global stability of $E_1$.

Theorem 6. The existing state $E_1$ of System (2.2) is globally stable if (4.2) and

Proof. We assume that $E_1$ exists. For Eq (4.2), we repeat the reasoning from the proof of Theorem 5 and obtain the following:

Now, we focus on the convergence of sequence $I_n^{(1)}$. We again rely on the approach from the proof of the previous theorem. We define the function as follows:

In the above definition, we treat variable $I_2$ as any constant from the interval $\Big[0, \frac{C_2}{\sigma_2} \Big]$. Differentiating $\Psi$ with respect to $I_1$ at $I_1 = 0$ gives the following:

Assume that

Then, $\Psi' (0) > 1$ and point $I_1 = 0$ cannot be reached for any iteration of $\Psi$. Let us denote the smallest fixed point of $\Psi$ from interval $\Big[0, \tfrac{C_1}{\sigma_1} \Big]$ by $I_1^*$. See that $\Psi \left(\tfrac{C_1}{\sigma_1} \right) = (r_1 - \alpha_1) (1 - \gamma_1) \tfrac{C_1}{\sigma_1} < \tfrac{C_1}{\sigma_1}$. From the Intermediate Value Theorem, we get the existence of positive fixed point $I_1^* \in \Big(0, \tfrac{C_1}{\sigma_1} \Big)$ for which $\Psi(I_1^*) = I_1^*$, $\Theta(I_1) > I_1$ for $I_1 \in (0, I_1^*)$; as a consequence, $\Psi(I_1^*) \leq 1$. Similarly, as for $\Theta''(I_2)$ in the proof of Theorem 5, we have $\Psi''(I_1) < 0$, which implies that $\Psi'(I_1) < \Psi'(I_1^*) \leq 1$ for $I \in \Big(I_1^*, \tfrac{C_1}{\sigma_1} \Big)$. Clearly for such $I$, we have $\Psi(I_1) < I_1$. Hence, there exists a unique positive fixed point $I_1^* \in \left(0, \tfrac{C_1}{\sigma_1} \right)$ of $\Psi$. It is clear from the analysis of the existence of stationary states that appear in System (2.2) that $I_1^* = \widehat{I}_1$ from the $E_1$ state. Thus, we can write $\lim\limits_{n \to \infty} I_n^{(1)} = \widehat{I}_1$ and, using Eq (3.22), $\lim\limits_{n \to \infty} S_n^{(1)} = \widehat{S}_1$. Combining these convergences with Eq (5.9) yields the global stability of $E_e$ under the fulfillment of Eqs (4.2) and (5.12).

Observe that Eq (5.12) must hold for any $I_2 \in \Big[0, \frac{C_2}{\sigma_2} \Big]$. Therefore, from the properties of function $G$, we replace this inequality with Eq (5.8).

Observe that Eq (5.8) is stronger than Eq (3.23), being conditions for the $E_1$ existence. Hence, we do not provide any additional constrain Eq (5.8). Now, compare Theorems 4 and 6. The latter one provides less conditions for the global stability of $E_1$ than Theorem 4 stating the local stability of $E_1$.

The following theorem concerns global stability of $E_e$:

Theorem 7. The existing state $E_e$ of System (2.2) is globally stable if (5.8) and

Proof. Assuming the existence of $E_e$, we refer to the proof of Theorems 5 and 6. We use the function $\Theta: \Big[0, \frac{C_2}{\sigma_2} \Big] \to \Big[0, 1 + \frac{C_2}{\sigma_2} \Big]$ with the formula given by Eq (5.3) from the proof of Theorem 5. If (5.13), then $\Theta'(0) > 1$ and the fixed point $I_2 = 0$ is unstable. We repeat the reasoning conducted for function $\Psi$ from the proof of Theorem 6 and obtain the following:

Now, we again apply the function $\Psi$ defined by (5.10) from the previous proof. Following the analysis from that proof, we state that under Eq (5.8), being the strengthened condition (5.12) for any $I_2 \in \Big[0, \frac{C_2}{\sigma_2} \Big]$, we have the following:

Merging (5.14) and (5.15) provides the thesis of Theorem 7 for Eqs (5.13) and (5.8).

From properties of function $G$ and Eq (3.11), we get that Eq (5.8) is stronger than Eq (3.26) for the $E_e$ existence. Hence, there is no need to restrict Eq (5.8) in Theorem 7.

Furthermore, by comparing that Corollary 3 and Theorem 7 yield the global stability of $E_e$, we get less conditions than for the local stability of this state.

Now, let us compare the conditions for the local stability of the stationary states that appear in Theorems 5–7. Equations (5.1) and (5.8) from Theorems 5 and 6, respectively, are contradictory. The analogical conclusion holds for Eqs (4.2) and (5.13) from Theorems 6 and 7, respectively. Moreover, each inequality from Theorem 5 has its opposite counterpart from Theorems 7. It concludes that one stationary state is always globally stable at most, which is what we obviously expect. Moreover, if both inequalities

hold, then there is no global stability in the System (2.2).

6.   The basic reproduction number

Now, we compute the basic reproduction number $\mathcal{R}_0$ of System (2.2). As $\mathcal{R}_0$, we define a number of new infections incidences produced by one infectious individual in a population at the disease-free stationary state ^[15]. In order to compute $\mathcal{R}_0$, we will use of the next-generation approach, which was introduced in that paper. We consider general forms of functions $G$ and $H$ with properties from Section 2.

First, we arrange System (2.2) so that the first equations are for the infected variables. We obtain the following:

with the disease–free equilibrium having the following form:

The Jacobian matrix of System (6.1) at $\overline{E}_{df}$ reads as follows:

We express $\overline{J}(\overline{E}_{df})$ as the following block matrix:

where $F$, $T$, $\overline{J}_1$, and $\overline{J}_2$ are $2 \times 2$ matrices. Matrix $F$ concerns new infections and matrix $T$ reflects the other processes for the infective stages. We write them as follows:

The needed matrix is $F \cdot (I-T)^{-1}$ that reads as follows:

We define $\rho(M)$ as a spectral radius of a matrix $M$. In the next generation method, the inequalities $\rho(\overline{J}_2) < 1$ and $\rho(T) < 1$ must hold. They can be written as $\max(r_1, r_2) < 1$ and $\max(1 - \kappa_1, 1 - \kappa_2) < 1$, respectively. Clearly, both equalities are always true. We eventually obtain the following:

Let us express Theorem 3 in the context of $\mathcal{R}_0$. We formulate the following corollary:

Corollary 4. $E_{df}$ of System (2.2) is locally stable if $\mathcal{R}_0 < 1$.

7.   Numerical simulation

Here, we present a numerical simulation conducted for System (2.2). The simulation concerns the case of $TB$ described in Section 1. The values of the parameters besides those for the transmission coefficients are included in Table 1. Each value is the arithmetical mean of the corresponding ones for every year between the years 2000–2024. For computations, we used the annual data from ^[16]. The units of the parameters reflect the case of the corresponding continuous–time systems analyzed in our previous paper ^[14]. Model (2.2) is the inherently discrete–time system and does not concern an explicit change of values of the variables in time. Hence, the parameters in System (2.2) are dimensionless.

There is no explicit medical definition of the transmission coefficient. For this reason, there is a need to estimate the values of these coefficients. For this purpose, we used the built–in lsqcurvefit function in Matlab that is based on the Gauss–Newton algorithm ^[17]. The best–fitted values of the transmission coefficients were obtained thanks to a comparison of simulated values to the actual data. Table 2 shows the estimated values of these coefficients, which units reflect the analogical continuous–time system.

The simulated results and the associated actual data are shown in Figure 1. Observe that the simulations and the actual data for years 2022–2024 show opposed courses. In the last years, we observe the worldwide expansion of tuberculosis. However, the values of the epidemiological coefficients in our paper reflect over a 20-year period.

8.   Conclusions

In this paper, we introduced and analyzed the discrete–time system of an epidemic spread in a heterogeneous population. In this population, we indicated two subpopulations: a low ($LS$) risk and a high ($HS$) risk of getting infected. We constructed this system without the discretization of its continuous–time counterpart. The proposed system is analogical to those analyzed in our previous paper ^[13]. The difference considers a different function regarding the probability of staying healthy. In this paper, this function had a simpler form than previously used ^[13]. However, the chosen function, as in ^[13], was still not given with a specific formula, which rendered our model general.

We first indicated stationary states that appeared in the system with the condition of their existence. There are three states: disease–free ($E_{df}$), without infection in $LS$ ($E_1$), and endemic ($E_e$). We obtained the condition for their local and global stabilities. Moreover, the basic reproduction number $\mathcal{R}_0$ of the system was computed. In order to obtain explicit results, we also investigated the case when the functions that reflected the probability of remaining healthy were linear. We want to emphasize that we considered different values of the given parameter for each subpopulation in our analysis. This assumption made the population entirely heterogeneous.

The mathematical analysis provided analogical results to those from ^[13], which were expected from the epidemiological point. $E_{df}$ state was locally stable if $\mathcal{R}_0 < 1$ and lost its stability for $\mathcal{R}_0 > 1$. Additionally, the conditions for the local stability of $E_1$ and $E_e$ were provided. Moreover, we managed to provide conditions for the global stability of each stationary state.

Finally, we conclude that our model can be applied to the epidemiological modeling of a heterogeneous population for the discrete nature of epidemic dynamics. The structure of the model eliminates the problem of the step size of the discretization method, which has no clear biological meaning.

Use of AI tools declaration

The author declares he has not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The author declares there is no conflict of interest.