Study of a singularly perturbed tuberculosis model E. Tropkina1, E. Shchepakina1 1 Samara National Research University, 34 Moskovskoe Shosse, 443086, Samara, Russia Abstract A detailed analysis of the dynamic model of the tuberculosis epidemic was conducted. It is shown that the dynamic model contains several time scales and can be represented in a singularly perturbed form. With help of the integral manifolds theory and the reduction principle, the reduction of the modeling system was justified and carried out. This approach allow us to replace the original system by another system of a lower order on an integral manifold whose dimension is equal to that of the slow subsystem and which, at the same time, preserves the essential properties of the original system. The conditions for the stabilization of the epidemiology based on the selection of necessary treatment and preventive measures are determined. Keywords: singular perturbations; integral manifold; order reduction; stability; epidemiology; tuberculosis 1. Introduction It is well known that tuberculosis is a deadly disease, the fight against which is still relevant today. In 1882, Robert Koch, along with the discovery of the tubercle bacillus, found that this disease is transmitted by aerogenic means [1]. Consequently, people who frequently contact individuals with an active form of tuberculosis (the infectious stage of the disease) have a much higher risk of infection. Majority of infected people remain latent carriers throughout their life. The average duration of the latent period (the period of latent infection) ranges from several months to dozens of years. However, the risk of progression to the active form of tuberculosis increases dramatically in the presence of concomitant infections that weaken the immune system. In the absence of treatment for tuberculosis of respiratory organs, mortality is about 50%. Currently, about 3 million people die from tuberculosis every year worldwide [2]. But in most cases, tuberculosis is curable. The current methods of treating this disease require long-term courses of treatment (from six months to several years), the violation of which often leads to the return of the disease and the development of drug resistance. Treatment of drug-resistant tuberculosis is much more dangerous for the patient, less successful and more costly than treatment of the disease caused by common strains of the pathogen. All this shows the necessity of developing a new approach of identifying and treating patients with tuberculosis, as well as timely and effective treatment. One of the most effective methods for solving such problems is the construction and study of a mathematical model describing the processes of infection spread in the population, the development of the disease and the impact of anti-tuberculosis measures. Based on this study we can determine effective measures to combat this dangerous phenomenon. In the present work, a detailed analysis of the dynamic model of the tuberculosis epidemic based on the cluster approach has been carried out. The presence of several time scales made it possible to apply the geometric theory of singular perturbations for its qualitative study. This approach made it possible to determine the conditions for the stabilization of epidemiological conditions based on the selection of the necessary treatment and preventive measures. 2. Cluster model of the tuberculosis epidemic Mycobacterium drops of tuberculosis get into the air when coughing or sneezing infected people. Tuberculosis bacillus, spreading by such drops, lives in the air for a short period of time (about two hours) and, therefore, it is believed that casual contact with persons with active tuberculosis (infected individuals) rarely lead to the spread of the disease, and that most relapses are the result of prolonged and close contacts with primary carriers of infection. Latently infected individuals become infectious after some, usually long, time period. This period of transition to an active form of infection is called latent. Latent periods vary from several months to dozens of years. Most infected individuals never go to the active form of tuberculosis. On the other hand, the average length of the infection period is relatively short (several months). This indicator is decreasing in countries with affordable treatment. A common scheme for analyzing the spread of the tuberculosis epidemic is based on the division of the population into certain classes. Consider the basic mathematical model of the spread of tuberculosis [2]:  dS I  dt      S  1  S  N ,   dE    S  I  (   k  r )  E    E  I    R  I ,  dt 1 N 1 2 N 3 N  (1) dI I   k  E    E   (  d  r )  I ,  dt 2 N 2  dR I   r2  I  r1  E    R  3  R  .  dt N 3rd International conference “Information Technology and Nanotechnology 2017” 117 Mathematical Modeling / E. Tropkina, E. Shchepakina Model (1) was derived under assumption that the total population is divided into four epidemiological classes: susceptible (S), i.e., uninfected, but susceptible to infection individuals, those in whose body the pathogens of tuberculosis have not yet penetrated; carriers of latent infection (E), i.e., individuals in whose body tuberculosis pathogens are present, in equilibrium with the immune system, such individuals are characterized by the absence of any external manifestations of the disease; infected individuals (I) with clinical manifestations of tuberculosis caused by sufficiently extensive tissue damage as a result of the activity of mycobacteria in their organisms; recovery individuals (R) who have completed treatment and recovered from the disease. The parameter  reflects the influx of young people into the model population The parameters 1 , 2 , 3 are the transmission rates of tuberculosis infection for the respective classes;  is natural mortality rate; k is the tuberculosis progression rate; d is the tuberculosis induced mortality rate; r1 , r2 denote the treatment rates for latent class and infectious class, respectively; N  S  E  I  R is the total population size. The mixing of classes in this model is homogeneous, that is, there are no assumed differences between individuals while tuberculosis transmission depends on the rate of infection. The basic reproduction number, one of the most important characteristics in mathematical biology, defined as the average number of secondary cases produced by typical infected individuals, mainly in a susceptible population, is described by 1 k k 0HM    Q0  , (   d  r2 ) ( k    r1 ) (k    r1 ) 1 Q0  ,     d  r2 ,  where Q0 is the number of secondary latent infected individuals produced by a typical infected individual during the mean infectious period 1  , f  k / (k    r1 ) shows the probability of survival during the transition from the latent to the active infectious stage. It is assumed that only individuals who have frequent and prolonged interactions with infected people have a high risk of tuberculosis infection. Newly infected individuals activate clusters (groups of individuals who come into regular and close contact with people in an active form of tuberculosis), increasing the risk of tuberculosis infection for susceptible individuals of each cluster [2]. Let us make a number of variables changes. We set a constant n be the average size of the cluster; the risk of infection with tuberculosis in the cluster will be determined by the parameter  . Further, the population of uninfected people in the cluster will be given by N1 (t )  n  I (t ), where N1 includes two subpopulations, namely, susceptible S1 and latent infected E1 , i.e., N1 (t )  n  I (t )  S1 (t )  E1 (t ). The population of persons who do not belong to the cluster at the time t is denoted as N2 . This population consists only of . the sensitive  S2  and the latent infected  E 2  individuals. The subpopulations of the sick are not included in this model for the sake of simplicity. We assume that n  k  E2 (S2 / N2 ) individuals go to the class S1 per unit of time, while the individuals n  k  E2 ( E2 / N2 ) go to the class E1 per unit time. In addition, since infected individuals are cured or die (with a speed   I ), that   I is the rate at which clusters become inactive (or die). We suppose that n    I  S1 / N individuals returns to the class S2 and n    I  E1 / N returns to the class E2 for the unit of time, respectively. Assuming a low level of distribution of individuals with an active form of tuberculosis, we consider the case N1 N 2 , hence we can neglect the birth rate in the population. The above assumptions lead (1) to the following basic cluster model [2]:  dS1 S2  dt  (    )  S1  N  n  k  E2 ,  2  dE1 E2     S1    E1   n  k  E2 ,  dt N2  dI   k  E2    I , (2)  dt  dS 2 S2  dt      S 2    S1  N  n  k  E2 ,  2  dE2 E     E1  (   k )  E2  2  n  k  E2 .  dt N 2 The basic reproduction number for system (2) is:  n k c0    Q0  f . (   ) (  k ) Hence, Q0    n / (   ). is the expected number of infections produced by one infected individual in his cluster. Only a part f  k / (  k ) among the infected individuals will survive during the latent period. 3. Dimensionless model Disease and dynamics of populations have characteristic time scales. Dynamics of influenza in humans is "super fast" at the individual and community levels in comparison with the dynamics of the carriers of infection (people). This is because the 3rd International conference “Information Technology and Nanotechnology 2017” 118 Mathematical Modeling / E. Tropkina, E. Shchepakina average life expectancy of an infected person is approximately 4000-8000 times the average duration of an influenza infection and, therefore, the largest outbreaks of influenza occur in local communities before any significant demographic change can occur (several months). Consequently, when studying the dynamics of the flu epidemic, two typical time scales are often used: the time scale of the disease and the life span of the carrier of the infection [3]. Tuberculosis is usually described as a slow disease due to its long and varied distribution of the latent period and due to the short and relatively narrow distribution of its infectious period [1]. Most latently infected with tuberculosis do not become actively infected, i.e., there is no transition of latent form to active. Some become actively infected during a five-year period, while others become active only after a longer period of time (perhaps decades). On the other hand, infected individuals remain so for relatively short periods of time, in part because of the use of antibiotics (an average of six months). Since secondary infections are formed from infected individuals, individuals with an active form of tuberculosis have a relatively short period for possible infection of other people. Consequently, the infection of tuberculosis-susceptible individuals occurs on the same time scale as the recovery of persons with active tuberculosis. The disease progresses, moving from the latent stage to the active one. This occurs on a time scale that is of the same order as the average lifespan of the carrier of the infection (human). Tuberculosis can be acquired at random (individual level), that is, as a result of accidental contacts, or through members of an epidemiologically active cluster that includes at least one actively infected person. The choice of these levels of distribution is not accidental, it is associated with the observed statistical data of the spread of tuberculosis. Since in the original system there is no clear separation into fast and slow variables, but the process speeds have different orders, it is necessary to bring the system (2) to a dimensionless form. To do this, we introduce the following dimensionless variables and parameters: d S E    S1   k  t , dt  , x1  2 , x2  2 , y1   , k   k     E1   I k   y2   , y3   ,  , m , B , k  k      k where  is a small positive parameter. With new variables system (2) has a form:  dx1 x1  x2  d  B  (1  x1 )  (1  m)  y1  n  x  x ,  1 2  dx2 x 2   (1  m)  y2  (1  B)  x2  n  2 ,  d x1  x2  dy x  x     y1  n  1 2 , 1 (3)  d  x1  x2  dy x2   2  m  y1  (1  m)  y2  n  2 ,  d x1  x2    dy3  x  (1  m)  y ,  d 2 3 where y1 , y2 and y3 are the fast variables, while x1 and x2 are the slow variables; x1 corresponds to a population of sensitive individuals not belonging to the cluster; x2 is a population of latently infected individuals not belonging to the cluster; y1 reflects a population of susceptible people in the cluster; y2 is a latently infected population; y3 is the population of infected individuals belonging to the cluster. The generated system [4, 5] for (3) is:  dx1 x1  x2  d  B  (1  x1 )  (1  m)  y1  n  x  x ,  1 2  dx2 x22   (1  m)  y2  (1  B)  x2  n  ,  d x1  x2  x1  x2 0   y1  n  , (4)  x1  x 2  x22 0  m  y1  (1  m)  y2  n  ,  x1  x2  0  x2  (1  m)  y3 .  The last three equations in (4) determine the zeroth order approximation of the slow manifold (the slow surface) of system (3). From these equations we have the expressions for y1 , y2 and y3 : 3rd International conference “Information Technology and Nanotechnology 2017” 119 Mathematical Modeling / E. Tropkina, E. Shchepakina x(t )1  x2 (t ) y1 (t )  n  , x1 (t )  x2 (t ) x2 (t ) Q  x (t )  n  x2 (t ) (5) y2 (t )  n   0 1 , x1 (t )  x2 (t ) 1 m x2 (t ) y3 (t )  . (1  m) Since the Jacobi matrix of the fast subsystem of system (3):  g1 g1 g1     y1 y2 y3   1 0 0   g 2 g2 g 2    D    m (1  m) 0   y1 y2 y3    g 0 0 (1  m)  g3 g3   3   y   1 y2 y3  has the negative eigenvalues, then the slow surface (5) is stable [4, 5], hence, we can replace the original system (3) by the reduced one. The reduced system is the projection of the original system on the slow surface (5) with preservation of the essential qualitative features of the dynamics of the complete system (see, for example, [6-17]). The reduced system has the form  dx1 x1  x2  d  B  (1  x1 )  Q0  x  x ,  1 2  (6)  dx x  x 2  Q0  1 2  (1  B )  x2 ,  d x1  x2 where Q0  n  m    n (   ) is a number of secondary infections produced by one infected individual in a population each member of which is susceptible. System (6) is a homogeneous mixed model in which the infection spread parameter is a function of the parameters: Q0    n (   ). Recall that the basic reproduction number is defined as c0  Q0  k (  k ). It is easy to see that  c0 is the threshold parameter for the dynamic model (6). 4. Analysis System (6) has two equilibriums P1 and P2 :  B B  B 2  BQ0  P1 (1, 0), P2  , .  Q0  1 (1  B)(Q0  1)  Our goal is to find such conditions on the parameters of the system, reflecting the methods of treatment and preventive measures, in which the infected and latent infected populations are stabilized at the minimum value. From a mathematical point of view, this means that the equilibrium of the system is globally asymptotically stable, and their coordinates x1 , x2 should be as small as possible. It should be noted that the coordinates of the point P1 correspond to the following situation in real life: all persons are susceptible to tuberculosis, but not infected (either in active or latent form). In other words, this is the most favorable situation from the point of view of epidemiology. Therefore, if this point is globally asymptotically stable, then this will be the most favorable outcome. The Jacobi matrix of the system (6)  x22 x12    B  Q0  Q 0   ( x1  x2 ) 2 ( x1  x2 ) 2 J     (7) x22 x12  Q0  (1  B )  Q0    ( x1  x2 ) 2 ( x1  x2 ) 2   at the point P1 (1,0) is  B Q0  J P1   .  0 (1  B )  Q0  with the eigenvalues 1  В, 2  1  В  Q0 . Thus, according to the reduction principle, the condition Q0  1  В is the criteria for the globally asymptotic stability of the point P1 . This condition can be written in terms of the main reproductive number as c0  1. The Jacobi matrix (7) at the equilibrium point P2 has the form 3rd International conference “Information Technology and Nanotechnology 2017” 120 Mathematical Modeling / E. Tropkina, E. Shchepakina  В 2  В (Q0  2)  (Q0  1) 2  (1  В)2    Q0 Q0 J P2   . (1  В  Q0 ) 2 (1  B)(1  В  Q0 )     Q0 Q0    For the asymptotic stability of the equilibrium P2 , it is necessary and sufficient that the trace of this matrix be negative and the determinant positive, i.e.,  В 2  В(Q0  2)  (Q0  1) 2  (1  B)(Q0  1  В)  0,  2  В  В(Q0  2)  (Q0  1)  (1  B)(Q0  1  В)  (1  В) (Q0  1  В)  0. 2 2 2 The last system can be written as  Q0  1 Q0  0,  (1  В)(Q0  1  В)  Q0  1 Q0  0. Taking into account the physical meaning of the parameters, this condition is equivalent to the inequality Q0  1  B. Thus, we have the following statement. Theorem. If c0  1 a disease-free equilibrium P1 (1, 0) (i.e., the point determining the absence of infection in the population) is globally asymptotically stable. If c0  1, that point P1 (1, 0) is unstable and the equilibrium  B B  B 2  BQ0  P2  ,   Q0  1 (1  B )(Q0  1)  is globally asymptotically stable. It should be noted that in [1] this statement was obtained on the basis of the Hoppenshtadt theorem. Although the singular point P2 corresponds to the situation when the infected individuals in the population are present, but for typical values of the parameters, the ratio of the quantity of latently infected individuals to the amount of susceptible individuals is insignificant. Hence, the situation when this point is asymptotically stable is not so bad from the point of view of the real situation. In other words, the latently infected individuals will exist but their quantity will not be so high, that is, the epidemic threshold will not be reached. 5. Correctness of model reduction The solutions of the original and reduced systems were plotted with the help of Wolfram Mathematica 10.3 and NetBeans 8.0.1. Figures 1-4 show the graphs for systems (3) and (6) with the following parameters values:   1 / 60,   2,     10 5 , n  20,   1. Fig. 1. The solutions of system (3) with ε=0.00053. Using a program written in the Java language in the NetBeans 8.0.1 environment, the errors in the deviation of the plots of the solutions to the complete and reduced systems were determined. For the three cases considered, the errors between the graphs of the solutions of the original and reduced systems are 0.00010, 0.00012 and 0.00015, respectively. In Figure 5 one can see the solutions of the original and reduced systems. The figure clearly demonstrates the almost complete coincidence of the solutions plots, which confirms the correctness of the reduction performed. Thus, the conclusions drawn about the qualitative behavior of the solutions to the reduced system can be transferred to the original model (3). 3rd International conference “Information Technology and Nanotechnology 2017” 121 Mathematical Modeling / E. Tropkina, E. Shchepakina Fig. 2. The solutions of system (3) with ε=0.00204. Fig. 3. The solutions of system (3) with ε=0.00476. Fig. 4. The solutions of system (6). 6. Conclusion In the present work, the tuberculosis model has been investigated via methods of qualitative analysis. It has been shown that the model has several time scales, so it can be represented as a singularly perturbed system of ODEs. The reduction of the 3rd International conference “Information Technology and Nanotechnology 2017” 122 Mathematical Modeling / E. Tropkina, E. Shchepakina system has been carried out, as a result of which, the original system of five differential equations was replaced by its projection on the slow integral manifold. It should be noted that, due to the stability of the slow integral manifold, this reduction is correct, and the reduced system of two differential equations preserves the essential properties of the original model. The conditions under which the system has the globally asymptotically stable equilibrium have been found. This result means that under the appropriate selection of treatment and preventive measures, the spread of the tuberculosis epidemic can be completely suppressed. Fig. 5. The solutions x1 (t ), x2 (t ) of systems (3) and (6) with ε=0.00053. Acknowledgements This study was supported by the Russian Foundation for Basic Research and Samara region (grant 16-41-630529-p) and the Ministry of Education and Science of the Russian Federation as part of a program of increasing the competitiveness of SSAU in the period 2013–2020. References [1] Song B, Castillo-Chavez C, Aparicio JP. Tuberculosis models with fast and slow dynamics: the role of close and casual contacts. Mathematical Biosciences 2002; 180: 187–205. [2] Castillo-Chavez C, Song B. Dynamical models of tuberculosis and their applications. Мathematical biosciences and engineering 2004; 361–404. [3] Vynnycky E, Fine PE. The long-term dynamics of tuberculosis and other diseases with long serial intervals: implications of and for changing reproduction numbers . Epidemiol Infect 1998; 121(2): 309–324. [4] Sobolev VA, Shchepakina EA. Model Reduction and Critical Phenomena in Macrokinetics. Moscow: Energoatomizdat Publisher, 2010; 320 p. (in Russian) [5] Shchepakina E, Sobolev V, Mortell MP. Singular Perturbations: Introduction to System Order Reduction Methods with Applications. Springer Lecture Notes in Mathematics 2014; 2114: 212 + XIII p. [6] Strygin VV, Sobolev VA. Effect of geometric and kinetic parameters and energy dissipation on orientation stability of dual-spin satellites. Cosmic Research 1976; 14(3): 331–335. [7] Gavin C, Pokrovskii A, Prentice M, Sobolev V. Dynamics of a Lotka-Volterra type model with applications to marine phage population dynamics. J. Phys.: Conf. Series 2006; 55(1): 80–93. [8] Pokrovskii A, Shchepakina E, Sobolev V. Canard doublet in a Lotka-Volterra type model. J. Phys.: Conf. Series 2008; 138: 012019. [9] Sazhin SS, Shchepakina E, Sobolev V. Order reduction of a non-Lipschitzian model of monodisperse spray ignition. Mathematical and Computer Modelling 2010; 52(3-4): 529–537. [10] Pokrovskii A, Rachinskii D, Sobolev V, Zhezherun A. Topological degree in analysis of canard-type trajectories in 3-D systems. Applicable Analysis 2011; 90(7): 1123–1139. [11] Sobolev VA, Tropkina EA. Asymptotic expansions of slow invariant manifolds and reduction of chemical kinetics models. Computational Mathematics and Mathematical Physics 2012; 52(1): 75–89. [12] Korobeinikov A, Archibasov A, Sobolev V. Order reduction for an RNA virus evolution model. Journal Mathematical Biosciences and Engineering 2015; 12(5): 1007–1016. [13] Archibasov AA, Sobolev VA, Korobeinikov A. Asymptotic expantions of solutions in a singularly perturbed model of virus evolution. Computational Mathematics and Mathematical Physics 2015; 55(2): 240–250. [14] Korobeinikov A, Archibasov A, Sobolev V. Multi-scale problem in the model of RNA virus evolution. J. Phys.: Conf. Series 2016; 727(1): 012007. [15] Sobolev VA. Canards and the effect of apparent disappearance. CEUR Workshop Proceedings 2015; 1490: 190–197. [16] Lapshova MA, Shchepakina EA. Study of the dynamical model of HIV. CEUR Workshop Proceedings 2016; 1638: 600–609. [17] Korobeinikov A, Shchepakina E, Sobolev V. Paradox of enrichment and system order reduction: bacteriophages dynamics as case study. Math Med Biol. 2016; 33(3): 359–369. 3rd International conference “Information Technology and Nanotechnology 2017” 123