SIS CRISS-CROSS MODEL OF TUBERCULOSIS IN HETEROGENEOUS POPULATION

. In this paper we propose a model of tuberculosis (TB) transmission in a heterogeneous population consisting of two diﬀerent subpopulations, like homeless and non-homeless people. We use the criss-cross model to describe the illness dynamics. This criss-cross model is based on the simple SIS model with constant inﬂow into both subpopulations and bilinear transmission function. We ﬁnd conditions for the existence and local stability of stationary states (disease-free and endemic) and ﬁt the model to epidemic data from Warmian-Masurian Province of Poland. Basic reproduction number R 0 is considered as a threshold parameter for the general model. Applying local center manifold theory we show that when R 0 = 1 a supercritical bifurcation occurs, and with R 0 increasing above this threshold the disease-free stationary state loses stability and locally asymptotically stable endemic stationary state appears. Our analysis conﬁrms the hypothesis that homeless individuals may be a speciﬁc reservoir of the pathogen and the disease may be transmitted from this subpopulation to the general population.


1.
Introduction. Tuberculosis (TB) is an airborne disease caused by Mycobacterium tuberculosis and is usually transmitted by individuals with active (smearpositive) TB [4]. Generally, the probability of becoming infected is low, that is why epidemics of TB should be analyzed over long windows in time [1]. The multiple drug treatment can be applied for treating most forms of TB. However, antibiotic resistant strains evolve fast, and therapy costs are high.
Parts of populations living in poverty are considered as the groups of high risk of infection. The correlation between incidences of TB and the poverty is thoroughly investigated and presented in the literature [16]. The standard of living of homeless people is generally low, that is why we consider them as the group with high probability of being infected. The homelessness escalates the risk of getting TB several times [2]. According to biological research, homeless people can be a pathogen reservoir and transmit the infection into the whole population [25]. What is interesting, the homeless shelters are good places for TB transmission since both susceptible and infectious people meet there [8]. Since TB can be transmitted from homeless individuals to non-homeless ones and conversely, we use criss-cross model.
Every criss-cross model is built up using component models. In our case the component model describes transmission of the disease within a single subpopulation. We would like to exploit simple SIS model as our component, because even for simple component models crossing them leads to complex models for which analysis could be difficult. The concept of using component models can be applied also for modeling TB epidemics dynamics. In [20] the authors based their model on the appearance of HLA-DR2 (human leukocyte antigen -DR isotype) allele in the human genome. This allele is strongly connected to pulmonary TB and its frequency is high in India, where the occurrence of TB is around 50% [24]. The model consists of two components where the first one is related to the people having susceptible genotype and the second one reflects individuals without this genotype. Component models can be naturally applied for modeling other kinds of TB, like in the models of cattle TB presented in [21], where each component corresponds to the single species (in [21] the populations of cows and badgers are considered).
Let us introduce the following component SIS model describing epidemic dynamics in a homogeneous population, which we divide into the subpopulations of susceptible (S) and infected (I) people. We would like to focus on the model based on the mass-action law. Hence, we consider the system of two differential equations (cf. [10] and [15]):Ṡ = C − βSI + γI − µS, where C is a constant inflow into the population, β is a transmission coefficient, γ stands for recovery rate, µ is a natural death coefficient and α reflects diseaserelated death. We assume that the coefficients are fixed and positive. However, in the following we consider subpopulation of homeless people and in such a case it would be good to reach a situation when C = 0. Note that in [23], C = 0 in the component model, however µ is not necessarily positive, as it covers both birth and death processes. It means that the model is based on the Malthusian law, and therefore it allows unbounded growth of the population. Let us make some scaling for system (1) to decrease the number of parameters appearing in the model. First, let us introduce the following new independent variable τ = γt for which we have 1 γ and the rest of the coefficients are scaled by γ. Next, we multiply the first and second equations of system (1) by β. As a result, we obtain where x = βS, y = βI, k = α + µ + 1, and C is scaled accordingly. It is clear that k > 1. Note that the variables x, y are functions of τ now. Furthermore, let us denote First, we analyze the basic properties of the model (cf. [15]). The form of the right-hand side of system (2) implies existence, uniqueness and positiveness of the solutions for positive initial conditions. If we add by sides the equations of system (2), we get which is a typical form of the equation used in epidemiological modeling. Note that due to the form of equation (3), the Malthusian properties do not appear in system (1) unlike the component model from [23]. Now, we estimate the size of the population. Since y 0, we have yielding Hence, if w(0) > C µ , then the population decreases, and therefore x(τ ) and y(τ ) are bounded above, as If w(0) C µ , then the population can increase, but the growth of the variables is limited. Clearly, where c = const > 0. Thus, x(τ ) and y(τ ) are defined for every τ > 0. Let us consider again equation (3), which together with the inequality y w gives w C − (µ + α)w.

Equation (3) with inequality (4) gives
which yields that Ω = w : w ∈ C α + µ , C µ is the invariant set. System (2) has two stationary states: disease-free E 0 = (x 0 , y 0 ) = C µ , 0 and endemic We see that the condition C µ > k implies that E e exists. It is obvious that E 0 exists independently of the values of the model parameters.
The Jacobian matrix for system (1), denoted by M, reads From the form of M we get that E 0 is locally stable if C µ < k. For C µ > k (i.e. when E e exists) E 0 is a saddle point. We also get that E e is always locally stable under its existence. Phase portraits with sample trajectories are presented in Fig. 1. Notice that in this system local stability implies global stability in R 2 + . Clearly, we  Figure 1. Nullclines C − xy + y − µx = 0, x − k = 0 and phase portraits for system (2), for C µ < k (A) and C µ > k (B).
can apply the Dulac-Bendixson Criterion to show that there is no close orbit for system (1): Hence, the Poincaré Theorem implies global stability -either E 0 is globally stable whenever E e does not exist, or E e is globally stable whenever exists.
The threshold condition C µk = 1 defines so-called basic reproduction number corresponding to system (2), that is R SIS 2. Simple criss-cross model of TB. Now, we turn to the main goal of this article, that is to criss-cross model of TB. We assume that the population can be divided into two homogeneous subpopulations based on individual behavior (cf., e.g. [13]), such that individuals in a given subpopulation are indistinguishable from one another. It means that parameters are identical for all individuals within a given subpopulation, but may vary significantly from subpopulation to subpopulation. Each group is then subdivided into epidemiological compartments. Thus, basing on the model described in the previous section, the dynamics of the disease is described by the following system of four differential equations: where β ij , i, j = 1, 2, are transmission coefficients between specific groups in heterogeneous population. In this paper, we assume that these coefficients are fixed and positive. Parameters α 1 and α 2 stand for disease-related death rates for the non-homeless and homeless individuals, accordingly; γ 1 , γ 2 are recovery coefficients for the non-homeless and homeless people, respectively; parameters µ 1 and µ 2 are natural death coefficients for the non-homeless and homeless individuals. Additionally, C i (i = 1, 2) reflects a constant inflow of humans into the given subpopulation. Here, they are numbers of newborn and net migrating individuals. We assume that µ i , C i > 0 for i = 1, 2. Notice that for simplicity we do not assume explicit exchange of people between the groups of homeless and non-homeless individuals.
Similarly to system (1), let us make some scaling of system (5) to reduce the number of parameters appearing in the model. First, let us introduce the following new independent variable τ = γ 1 t for which we have 1 γ 1 and the rest of coefficients are scaled by γ 1 . Let us denote β 1 = β 12 and β 2 = β 21 . Next, we multiply the first and second equations of system (5) by β 11 and the third and fourth ones by β 22 . As a result, we obtain where x 1 = a 1 S 1 , y 1 = a 1 I 1 , x 2 = a 2 S 2 , y 2 = a 2 I 2 , a i = β ii and β i and C i are scaled accordingly. We will write η instead of η 2 , as η 1 will not appear in the presented formulas. Note that the variables S i , I i are functions of τ now. Furthermore, let us denote Note that κ i > 1 and k i > η i > 0 for i = 1, 2.
Similarly to system (2), local existence and uniqueness of solutions of system (6) are a simple consequence of the form of the right-hand side of this system. Assume now that x i (0) and y i (0) are positive. If there exists a pointτ > 0 such that x i (τ ) = 0, then for the first such point we have This means that x i is repelled from 0, and moreover, for τ > 0 we have x i (τ ) > 0. Next, assuming that there is a pointτ > 0 such that y 1 (τ ) = 0 and the rest of the variables are positive, we obtain yielding non-negativity of y 1 , and similar inequality is valid for y 2 . Moreover, from equation (6b) or (6d) we obtain y i ≥ −k i y i and, by the Chaplygin-Perron differential inequality, we get Adding equations (6a) and (6b), as well as equations (6c) and (6d), we get Repeating the same arguments as for system (2) we obtain the invariant set and moreover, this set attracts all solutions of system (6). Thus, x i (τ ) and y i (τ ) are defined for every τ > 0. What is more, for populations described by system (6), the Malthusian properties do not appear as well.
If C 2 ≤ µ 2 k 2 , then ξ 2 = 0 and y 1 ( Figure 2. Graphic representation of solutions of system (10) when Red curves represent the graph of the function y 2 (y 1 ), blue ones represent the graph of y 1 (y 2 ). Dotted lines bound the region defined y j (0) = 0 and y j (0) = 0, for i, j = 1, 2, i = j. Hence, we can conclude that if at least one of the conditions: Fig. 2), then there exists the positive solution of system (10), which satisfies the following inequalities On the other hand, if both C i < µ i k i , then the positive solution of system (10) exists only if y 2 (0) ≤ 1 y 1 (0) . Thus, the positive stationary state exists if and does not exist otherwise (cf. Fig. 3). Summing up the analysis presented above we obtain the following corollary.
Corollary 1. The positive stationary state exists if at least one of the following conditions holds: Figure 3. Graphic representation of solutions of system (10) for (11) holds (A) and does not hold (B). Red curves represent the graph of the function y 2 (y 1 ), blue ones represent the graph of y 1 (y 2 ). Red line is the tangent line to the function y 2 (y 1 ) at zero, blue one is the tangent line to y 1 (y 2 ) at zero. Dotted lines bound the region defined by 3.2. Basic reproduction number and stability of the disease-free stationary state. In this subsection we calculate the basic reproduction number of system (6). The basic reproduction number, denoted by R 0 , is defined as the average number of secondary infections which a single infected individual would produce in a completely susceptible population [9]. Developed for the study of demographics, it was also studied for vector-borne diseases such as malaria, and directly transmitted human infections. It is now widely used in the study of infectious diseases. It determines whether or not a disease will spread through a population, and can be treated as a threshold (bifurcation) parameter. If R 0 < 1, then each infected individual will produce on average less than one new infected individual in its entire infectious period, and therefore the disease cannot invade the population. If R 0 > 1, then any infected individual will produce on average more than one new infected individuals, and the invasion is always possible. When considering a single infected compartment, the basic reproduction number is simply the product of the mean duration of the infection and the infection rate. In such a case, the condition R 0 < 1 also means that the disease-free stationary state is locally asymptotically stable and unstable otherwise [14]. For more complex models, in which several infected compartments are included, however, this simple heuristic definition of R 0 is insufficient. A more general definition of the basic reproduction number can be formulated as the number of new infections produced by a typical infective individual in a population at a disease-free stationary state [11]. We follow by the definition given by Diekmann et al. [9], where R 0 is taken as the spectral radius of the next generation matrix. In this approach, R 0 depends on the infected compartments of the given model.
According to the concept of the next generation matrix presented in [11], we sort the compartments in such a way that the first compartments correspond to infected individuals. Therefore, we introduce z = (y 1 , y 2 , x 1 , x 2 ). Using this notation, the disease-free stationary state reads z 0 = (ỹ 1 ,ỹ 2 ,x 1 ,x 2 ) = 0, 0, C1 µ1 , C2 µ2 and system (6) can be written as where We define the following 2 × 2 matrices Thus, the next generation matrix of system (12) has the following form Therefore, the basic reproduction number reads where denotes the spectral radius.
Let us now determine the stability conditions for the disease-free stationary state E 1 .
It is easy to verify that hypotheses (A1) − (A4) in [11, p. 31] are fulfilled. The hypothesis (A5) is satisfied if all eigenvalues of the matrix D(F − V)(z 0 ) for F(z) set to zero (the absence of new infections) have negative real parts, which is also easy to verify. Thus, by Theorem 2 [11] we have that the disease-free stationary state is locally asymptotically stable if R 0 < 1, but unstable if R 0 > 1, where R 0 is defined by formula (13). Thus, the following inequalities guarantee local stability of E 1 : Therefore, we have the following stability conditions: Proposition 1. The disease-free stationary state E 1 exists for C i > 0, i = 1, 2, and is locally stable for and is unstable when at least one of these inequalities does not hold.
Proof. We only need to show that inequalities (14) and (15) are equivalent. It is easy to notice that the second inequality in (14) or (15) is satisfied only when both C i < µ i k i , i = 1, 2, or when both C i > µ i k i . However, in the latter case neither the first inequality in (14) nor the first inequality in (15) is satisfied. Hence, E 1 is locally stable under the assumptions of the proposition.
Notice that if both C i ≥ µ i k i , then the first inequality in Proposition 1 does not hold, while if C 1 ≥ µ 1 k 1 and C 2 < µ 2 k 2 or vice versa, or both C i < µ i k i and (µ 1 k 1 − C 1 ) (µ 2 k 2 − C 2 ) ≤ β 1 β 2 C 1 C 2 , then the second inequality in Proposition 1 does not hold. This means that for R 0 ≥ 1 the conditions of Corollary 1 are satisfied and the positive endemic stationary state E 2 exists. Moreover, if this state exists, then the state E 1 loses stability.

3.3.
Local stability of the endemic stationary state. Now, we will check the conditions under which the state E 2 = (x 1 ,ȳ 1 ,x 2 ,ȳ 2 ) is locally stable.
We assume that E 2 exists. Forx i = 0, from (6a) and (6c) we have and from (6b) and (6d) we havē Hence, the Jacobian matrix for E 2 reads Calculating the characteristic polynomial we obtain where We will use the Routh-Hurwitz criterion. We build an auxiliary matrix If this matrix has all entries and main minors positive, then the state E 2 is stable according to the criterion. Let us assume that the value of C i are large enough that Then from equation (8) we havex 1 > 1 andx 2 > η. Let us assume that Then Conditions (17) are sufficient to ensure that all a i , i = 1, 2, 3, 4, are positive. Therefore, we can calculate subsequent main minors: 3 − a 4 a 2 1 . Note that for a 4 > 0 the condition ∆ 4 > 0 is equivalent to ∆ 3 > 0. Moreover, ∆ 3 > 0 reads and it is a stronger condition than ∆ 2 > 0. Thus, under assumptions (16) and (17), we get the stability condition Finally, we obtain the following Proposition 2. If β 1 and β 2 are sufficiently small such that inequalities (17) hold, as well as C 1 and C 2 are sufficiently large such that inequalities (16) hold, then the positive stationary state E 2 exists and is locally stable if inequality (18) holds.

The case without inflow into the subpopulation of homeless people.
In this section we consider equations (6) with C 2 = 0. This assumption means that there is no inflow of new individuals into the subpopulation of homeless people, and it is somehow an ideal situation to which local government should try to lead to. It is very desirable in developed societies, may be realistic but difficult to achieve. Under this assumption, the population of homeless people will go to extinction. Let us consider stationary states for this case. The stationary state E 1 changes to E 0 = (x 1 , 0, 0, 0), withx 1 = C1 µ1 . Let (x * 1 , y * 1 , x * 2 , y * 2 ) denote a positive stationary state. Adding the two last equations of Eqs. (6) with C 2 = 0 yields ηy * 2 − k 2 y * 2 − µ 2 x * 2 = 0. Because k 2 > η, we obtain that which means that the positive stationary state does not exist if C 2 = 0. Let us now consider the state E 3 = (x 1 ,ȳ 1 , 0, 0). From (6b) we have thatx 1 = k 1 , and then, by (6a), it follows thatȳ 1 = C1−µ1k1 k1−1 . Thus, the semi-positive stationary state E 3 exists only if C 1 > µ 1 k 1 . Now, we turn to the analysis of local stability of the states E 0 and E 3 . Notice that the form of Jacobian matrix J does not change for C 2 = 0 comparing to the case C 2 > 0. Hence, for E 0 we have three negative eigenvalues −µ 1 , −µ 2 and −k 2 , while the eigenvaluex − k 1 is negative under the assumption and the characteristic polynomial reads where P 1 (λ) = λ 2 + λ (ȳ 1 + µ 1 ) +ȳ 1 (k 1 − 1) , For the polynomial P 1 , there are either two real negative zeros or two complex zeros with negative real part. Because k 2 > η, we can estimate k 2 (β 2ȳ1 + µ 2 ) − β 2ȳ1 η = β 2ȳ1 (k 2 − η) + µ 2 k 2 > 0, and therefore P 2 has the same properties as P 1 . This means that there are four eigenvalues with negative real parts for J(E 3 ). Thus, the stationary state is always locally stable.
Proposition 3. If C 2 = 0 and C 1 < µ 1 k 1 , then there exists only one stationary state E 0 which is locally stable, while for C 1 > µ 1 k 1 there exists the semi-positive stationary state E 3 and then this state is locally stable.

Supercritical bifurcation and endemic state.
In this section the basic reproduction number R 0 will be treated as a threshold parameter and examined in details to investigate the stability conditions of the endemic stationary state near the threshold R 0 = 1. As it has been shown, the disease-free stationary state is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1. In general, for epidemic models there are two distinct bifurcations at R 0 = 1: forward (supercritical) and backward (subcritical). A forward bifurcation happens when there is no endemic stationary state near the disease-free stationary state which is locally asymptotically stable for R 0 < 1. Moreover a low level of endemicity occurs when R 0 is slightly above 1. For simple models it means that for R 0 < 1 there is no other stationary state except the disease-free one. On the other hand, a backward bifurcation occurs when the endemic stationary state exists for R 0 < 1 and disease-free stationary state may exist for R 0 > 1. Epidemiologically it means that reducing R 0 below the unity is no longer sufficient to guarantee disease elimination. This type of bifurcation allows multiple stable stationary states with fixed parameters, and large changes in stationary state behavior can be produced by small changes in parameters (cf. [12]). However, as we have shown, if there exists the endemic stationary state, then the disease-free stationary state loses stability. It means that in our case there is no backward bifurcation for R 0 < 1.
Applying the centre manifold theory we determine the existence and stability of the endemic stationary state near the threshold R 0 = 1. We follow the ideas from [11].
It is easy to verify that the Jacobian matrix for system (12) at z 0 has a simple zero eigenvalue for R 0 = 1. Hence, the result of Theorem 2 [11] can be applied. Thus, all other eigenvalues of the Jacobian matrix have negative real parts. Let v = (v 1 , v 2 , v 3 , v 4 ) T and u = (u 1 , u 2 , u 3 , u 4 ) T be left and right null vectors corresponding to the simple zero eigenvalue of the Jacobian at z 0 for R 0 = 1 such that vu = 1. Hence (J(z 0 ) − λI) u = 0 and J(z 0 ) T − λI v = 0.
If R 0 = 1, then there is a simple zero eigenvalue and both equations above can be written as a system of independent equations. Computations show that eigenvectors v and u can be chosen such that where v 1 , u 1 > 0 are such that vu = 1. Note that k 1 > 1, k 2 > η and C i < µ i k i , i = 1, 2, as R 0 = 1, thus v 1 , v 2 , u 1 , u 2 > 0 and u 3 , u 4 < 0. Let ξ be a bifurcation parameter such that R 0 < 1 for ξ < 0 and R 0 > 1 for ξ > 0. Let us consider the system z = f (z, ξ), where f = f (z, ξ) is the right hand side of system (12). We define The sign of a determines the nature of the endemic stationary state near the bifurcation point; cf. Lemma 3 and Theorem 4 in [11]. We can calculate that all second derivatives of f i appearing in formula (19) are zero at the disease-free stationary state except the following: Hence a = v 1 u 3 (u 1 + u 2 β 1 ) + v 2 u 4 (u 1 β 2 + u 2 ) < 0.
By Theorem 4 [11] we can conclude that there exists δ > 0 such that there is locally asymptotically stable endemic stationary state near z 0 for 0 < ξ < δ.
Proposition 4. For R 0 slightly greater than one there exists locally asymptotically stable endemic stationary state E 2 near the disease-free stationary state E 1 .
It means that a branch of super-threshold endemic stationary state exists and the direction of bifurcation is forward at R 0 = 1. The forward bifurcations can be also observed in a simple SIS model with and without recruitment and deaths, as well as in a simple multi-group SIS models with constant population size (cf. e.g. [17,14]). Backward direction of bifurcation can occur, in particular, due to vaccination or immunity. In such cases, reducing the value of R 0 below the unity is not a sufficient condition to eradicate the disease that is endemic due to hysteresis effect [5].
5. Numerical simulations. In this section we illustrate the dynamics of system (5) for the parameters fitted to the data from Warmian-Masurian province of Poland. We compared the actual epidemic data with the model to get the bestfitted parameters, which are summarized in Table 1. Demographic data used in Table 1. Parameters for the model described by system (5).

Name Definition
Value α 1 , α 2 Disease-related death rates 0.09 γ 1 , γ 2 Recovery coefficients 0.9 µ 1 , µ 2 Natural death rate 0.009 C 1 Constant inflow of humans into the subpopulation of the non-homeless 11000 C 2 Constant inflow of humans into the subpopulation of the homeless 60 β 11 Transmission coefficient 0.57566 · 10 −6 (estimated) β 12 Transmission coefficient 11.276 · 10 −6 (estimated) β 21 Transmission coefficient 3.7679 · 10 −6 (estimated) β 22 Transmission coefficient 98.249 · 10 −6 (estimated) the simulations were obtained from statistical yearbooks. The numbers of homeless were provided by the Regional Center for Social Policy, Office of the Marshall of the Warmian-Masurian province in the city of Olsztyn, Poland [19]. The epidemic data used for the study were anonymized by the Independent Public Tuberculosis and Lung Diseases Unit in Olsztyn, Poland, and only numerical results were used. The values of the parameters α 1 , α 2 , γ 1 , γ 2 , µ 1 , µ 2 , C 1 and C 2 are taken from [7]. A best-fit technique was used to estimate the values of the transmission coefficients β ij , i, j = 1, 2. All data are fully available without any restrictions. Comparison between the actual data and the criss-cross model is illustrated in Fig. 4. Graphs in Fig. 5 illustrate phase portraits for system (5) in the phase planes (S 1 , I 1 ) and (S 2 , I 2 ), for the specific parameter values summarized in Table 1. In this case the endemic stationary state E 2 exists and is locally stable.
6. Discussion. In this paper we have presented the analysis of the criss-cross model of tuberculosis dynamics in the heterogeneous population. The population is divided into two compartments with different risk of developing TB: non-homeless   Table 1.
and homeless people. The work is motivated by our previous results; cf. [6] as well as [23]. It should be pointed out that in the model presented and analyzed in this paper the Malthusian properties do not appear, as it was in [6,23]. It means that the situation when the subpopulation of non-homeless people goes to extinction while the subpopulation of homeless grows boundlessly, does not happen. We have investigated existence and local stability conditions for stationary states of the system. The basic reproduction number has been considered as a threshold parameter for the model. The disease-free stationary state is locally asymptotically stable if R 0 < 1 and loses stability for R 0 crossing 1. It means that if R 0 > 1, then the disease will spread and possibly persist within the host population, while if R 0 < 1, then the infection cannot maintain itself. Unfortunately, necessary conditions for stability of the endemic stationary state E 2 are difficult to determine. In the light of the above analysis it is concluded that: • if both C i < µ i k i , then either E 1 exists and is locally stable or both E 1 and E 2 exist, but E 1 is unstable, • in any other case E 2 exists and E 1 loses stability.
However, by the local analysis of the centre manifold theory we have obtained that a branch of super-threshold endemic stationary state exists near the disease-free stationary state. When R 0 is slightly larger than one, then there exists the endemic stationary state, which is locally asymptotically stable. It gives as a conclusion that reducing R 0 through one lowers the incidence of the disease until it is eliminated as R 0 passes below one. But, it is important to note that although the positive stationary state E 2 exists, it may be unstable, in particular, when one of recruitment rates is sufficiently large, whereas the other one is sufficiently small.
Note that for homogeneous population, the sufficient condition that guarantee disease elimination is C < kµ. This condition, however, does not guarantee that the disease will not spread in the heterogeneous population. If C i > µ i k i and C j < µ j k j , i = j, then the basic reproduction number is greater than one and the infection can spread. Moreover, if both C i < µ i k i , i = 1, 2, but (µ 1 k 1 − C 1 ) (µ 2 k 2 − C 2 ) < β 1 β 2 C 1 C 2 , the disease can still spread. Thus, our model suggests that disease transmission from one subpopulation to another may significantly affect the spread of the disease in a general population. Even if the basic reproduction number for one subpopulation is less than one and guarantee that the infection cannot maintain itself in the corresponding isolated subpopulation, the disease can be transmitted from one subpopulation to another and TB can invade the heterogeneous mixed population. It means that to control the spread of the disease in a heterogeneous population it is not enough to consider the disease dynamics in each subpopulation separately. In order to eliminate the infection it is required to consider a criss-cross dynamics of the disease spread.
Let us now return to the model described by system (5) and the original parameters. For the Warmian-Masurian Province of Poland, using the actual data and the best-fitted parameters (see Table 1), we get which is equivalent to conditions C i < µ i k i , i = 1, 2 for system (6). Moreover, we get that which is equivalent to (µ 1 k 1 − C 1 )(µ 2 k 2 − C 2 ) < β 1 β 2 C 1 C 2 . It means that the disease-free stationary state is unstable, whereas the endemic stationary state exists.
In fact, we get that in this case the endemic state is locally stable (cf. Fig. 5). The bifurcation diagram is shown in Fig. 6. Note that, if the value of C 2 is small enough to cross the black solid curve, then the endemic state bifurcates to the disease-free stationary state, which is locally stable. It means that regardless the transmission parameters the way to eliminate the disease spread in the general population is to limit the inflow into the subpopulation of homeless individuals. This confirms the hypothesis that homeless individuals may be a specific reservoir of the pathogen, and the disease may be transmitted from this subpopulation to the general population [25,18]. Figure 6. Bifurcation diagram for system (5). The solid line depicts the graph of µ 1 (γ 1 + α 1 + µ 1 ) − C 1 β 11 µ 2 (γ 2 + α 2 + µ 2 ) − C 2 β 22 − β 12 β 21 C 1 C 2 = 0. The dotted black lines depict C 1 = µ1(γ1+α1+µ1) β1 and C 2 = µ2(γ2+α2+µ2)

β2
. The red point depicts the point (C 1 , C 2 ) for the values of C 1 and C 2 taken from Table 1.
Note thatŷ i +µ i > 0,ŷ 1 (k 1 − 1) = y 1 (α 1 +µ 1 ) > 0 andŷ 2 (k 2 − η) = y 2 (α 2 +µ 2 ) > 0. Thus, the polynomial has four zeros with negative real parts. It means that regardless the model parameters, the stationary state E 4 is always locally stable. Thus, for β 1 , β 2 → 0, we can summarize the results as follows: • if both C i > µ i k i , then E 1 exists and is unstable, while E 4 exists and is locally stable, • if both C i < µ i k i , then E 1 exists and is locally stable, while E 4 does not exist, • in any other case E 1 exists and is unstable, while E 4 does not exist.
We therefore conclude that the (C 1 , C 2 ) parameter subspace is divided into four regions such that crossing from one region to another changes the stability of one stationary state. These results are summarized in Fig. 7. Note that if C i > µ i k i and C j < µ j k j , i = j, then for the first subpopulation there exists a locally stable endemic stationary state (x i ,ȳ i ),x i ,ȳ i > 0, whereas for the second subpopulation there exists a disease-free stationary state (x j , 0),x j > 0, which is locally stable. Figure 7. Bifurcation diagram for system (6) with β 1 , β 2 tending to zero. The solid black lines depict C 1 = µ 1 k 1 and C 2 = µ 2 k 2 .

7.
Conclusion. The increase of TB incidence by a high level of infectious population that have not been detected or treated causes that the challenge of TB control is becoming more and more important. Note that our model differs from the ones that have been widely studied in the previous years. The model considered by us is based on dividing the host population into two subgroups with different risk of developing TB and considering a criss-cross infection. Such class of epidemiological models can be further extended to many other classes of infective individuals and data for many other infectious diseases. What is more, we have fitted our model to the actual data and obtained best-fitted transmission parameters for both subpopulations. From the practical point of view, the model considered in this paper can be used to better understand the transmission behaviors of the disease and to forecast the disease trends, which can help to implement more effective preventive programs of TB (cf. e.g. [23,3,26]).