Modeling of Infectious Disease Dynamics under the Conditions of Spatial Perturbations and Taking into account Impulse Effects

Abstract

The infectious disease mathematical model by Marchuk for conditions of diffusion perturbations and taking into account impulse influences is generalized. The corresponding singularly perturbed model problem with delays is reduced to a sequence of problems without delay, for which the corresponding asymptotic developments of solutions are obtained. The results of numerical experiments characterizing the impulse effects of infectious disease factors on the immune response development in the conditions of spatially distributed diffusion perturbations are presented. The model decrease of the antigens maximum level in the infection epicenter due to their diffusion "erosion" in the viral disease process development is illustrated. Infectious disease model, dynamic systems, asymptotic methods, singularly perturbed ORCID: 0000-0001-5528-4192 (A. Bomba ); 0000-0002-8056-2980 (S. Baranovsky); ORCID: 0000-0001-7491-0371 (M. Pasychnyk);

1. Introduction

The simplest (basic) infectious disease model which describes the most general laws of organism humoral immune response to the viral antigens found is resulted in [ 1,2,3 ]. The infectious disease process development in the model is determined by the nonlinear differential equations system with delay, describing the rate of change in the viral antigens number, plasma cells, antibodies and the extent of damage to the target organ. The relatively small number of active factors of the infectious disease simplest model allows to establish strictly justified properties of its solutions, in particular, the stability of inpatient solutions. In [ 1 ] it was shown that the stationary solution, which describes the healthy organism state under certain conditions is asymptotically stable and retains this kind of resistance when infecting a healthy organism with a dose of antigen V 0 , that does not exceed a certain level V * of immunological barrier. The infectious disease basic model and its modifications in identifying their parameters according to clinical observations allow to predict the nature of the course and outcome of infectious disease, to investigate the general patterns of external influence on the process dynamics, analyze and evaluate various treatment procedures. The generalization of an infectious disease simplest model is antiviral and antibacterial immune response mathematical models[ 2,3 ]. In contrast to the simplest model, in addition to the humoral immune response with the antibodies production, the cellular type of immunity with the cytotoxic T-lymphocyte effectors accumulation is taken into account. As mentioned in [ 2,3 ], antibodies are able to neutralize viral antigens that circulate freely in the blood or lymph, but can’t penetrate into infected cells and neutralize viruses that multiply in them. Detection and destruction of infected cells is carried out by cytotoxic T-killer lymphocytes. The antiviral immune response model, as well as the basic model, is represented by a nonlinear differential equations system

2020 Copyright for this paper by its authors. with delay, describing the rates of change of free circulating in the body viruses Vf , antibodies F, infected with viruses of target cells CV, T-killer lymphocytes Е, T-helper lymphocytes of cellular immunity НЕ, T-lymphocytes helpers of humoral immunity HB, B-lymphocytes B, plasma P cells, stimulated MV macrophages and non-functioning part of the virus-affected organ m.

The infectious disease basic model, the antiviral immune response model and other immunology models [ 7,8,9 ] are built under the assumption that the "organism" is a homogeneous environment in which all the process components are instantly mixed and, as a result, evenly distributed. On the other hand, the antigens detection and the launch of the immune system of appropriate mechanisms of response to them does not occur immediately after the antigens penetration into the body. That is, some of the antigens that were not immediately neutralized by the immune system will penetrate into the cells, multiply in them and spread further in the body. As a result, infection foci with higher antigens concentrations forms around the affected cells. The antigens generated in the body will eventually be redistributed from the initial infection foci to the surrounding uninfected areas, increasing the affected area and decreasing the antigen concentrations values in the respective infection epicenters. An approach for taking into account small spatially distributed diffusion effects on population dynamics is presented in [ 4,5,6 ]. In particular, in [ 6 ] when modifying the simplest infectious disease model to take into account the impact on the dynamics of the disease of certain drugs introduced into the body, articles describing diffusion perturbations of the process active factors were added. The model decrease of the maximum antigens concentration in the infection epicenter due to their diffusion "erosion" in the infectious disease process development is illustrated. It is emphasized that even if the initial antigens concentration in some infection area exceeds a certain critical value (immunological barrier), diffusion "redistribution" for a certain period of time will reduce above critical values of antigen concentrations to below critical level, and their subsequent disposal can be provided with the level of organism immune protection available before infection.

But as in traditional infectious disease models of antiviral (antibacterial) immune response, and in their modifications [ 5,6 ], which take into account spatially distributed diffusion perturbations, impulse effects are not considered. The purpose of this work is to "fill" this gap.

2. Problem Statement

To take into account spatially distributed diffusion perturbations of infectious disease development, it is proposed to modify the basic models by introducing components that describe small diffusion spatially distributed effects ("redistributions"). Let us generalize this kind of modification of basic models by G. Marchuk by introducing additional terms describing the influences that are close to the impulse character. The spatio-temporal dynamics of the infectious disease process model components, taking into account close to impulse influences in the domain GZ  x,t :    x  ; 0  t   will be described by such a singularly perturbed dimensionless nonlinear differential equations system: in conditions vtx,t  h  h2 f x,tv x,t  h9vxx x,t  uV x,t,

1 stx,t  mh3v x,t   f x,t    h5 s x,t 1  2h10sxx x,t, ftx,t  h4 s x,t  f x,t  h8v x,t f x,t  h11 fxx x,t  uF x,t, mtx,t  h6vx,t  h7mx,t  2h12mxx x,t , sx,t0   s0 x , mx,t0   m0 x, vx,t   v0 x,t , f x,t   f 0 x,t , t0   t  t0 , (1) (2) where vx,t V x,t Vm , s x,t  С x,t C* , f x,t  F x,t F  ; V x,t , С x,t , F x,t , mx,t – respectively, the concentration of antigens, plasma cells, antibodies and the value of the damage degree to the target organ at point x at time t; Vm is some scale factor for the antigens concentration, for example, the biologically acceptable antigens concentration in the body; C* , F  are the plasma cells and antibodies concentration in a healthy organism; h1   , h2  F* , h3 Vm F* C* , h4   f , h5  C , h6  Vm , h7  m , h8  Vm ;  is the antigens reproduction rate;  is the coefficient that takes into account the result of the antigens interaction with antibodies;  is the period of time (delay) required to form a plasma cells cascade; С is value inverse to the plasma cells lifespan;  is immune system stimulation factor;  is the rate of antibodies production by one plasma cell;  f is the value inverse of the antibodies duration;  is the cost of antibodies to neutralize one antigen;  is the rate of the target organ cells damage;  m is the rate of the target organ recovery,  h9 ,  h11 ,  2h10 ,  2h12 are diffusion redistribution coefficients of antigens, antibodies, plasma and affected cells, respectively,  is a small parameter that characterizes the respective components small impact compared to the dominant components of the process. The function (m) takes into account the effect of reducing the antibody production productivity in significant damage to the target organ. If is the maximum value of the degree of the target organ damage, at which the immune system normal functioning is still possible, then on the segment 0  m  m* value  m is equal to one, regardless of the lesion, the immunological organs function fully. If m*  m 1 , the body efficiency is rapidly declining. Functions uV x,t , uF x,t , describing a close to pulse change, respectively, of the antigens and antibodies concentrations with maximum values at points xV j , xF j at times tV j , tF j will be represented as uV x,t  nV AV jeV jxxV j 2 eV j ttV j 2 , uF x,t  nF AF jeF jxxF j2 eF jttF j2

j1 j1 where AV j , AF j , V j ,  F j , V j ,  F j are parameters that determine the pulse intensity and "duration".

Note that in the future as initial, we will take, in particular, the conditions that characterize the stationary solution, which corresponds to the healthy organism state, namely

sx,t0  1, mx,t0   0, vx,t   0, f x,t  1, t0   t  t0 .

At the beginning we consider the case when the level of damage to the target organ by antigens remains such that it does not lead to a decrease in the productivity of antibody production,  m 1 . Then the solution of problem (1) - (2) with delay is reduced to a sequence of problems without delay [ 10 ]:

v0 x,t0   v0 x,t0  , f0 x,t0   f 0 x,t0  , t0  t  t0  ; ………………………………………………………………………………… vkt x,t  h  h2 fk x,tvk x,t  h9vkxx x,t  uV x,t ,

1  skt x,t  h3vk1 x,t   fk1 x,t    h5 sk x,t 1  2h10skxx x,t , s0 x,t0   s0 x , m0 x,t0   m0 x , v0t x,t  h  h2 f0 x,tv0 x,t  h9v0xx x,t  uV x,t,

1  s0t x,t  h3v0 x,t   f 0 x,t    h5 s0 x,t 1  2h10s0xx x,t,  f0t x,t  h4 s0 x,t  f0 x,t  h8v0 x,t f0 x,t  h11 f0xx x,t  uF x,t,  m0t x,t  h6v0 x,t  h7m0 x,t  2h12m0xx x,t,  fkt x,t  h4 sk x,t  fk x,t  h8vk x,t fk x,t  h11 fkxx x,t  uF x,t,  mkt x,t  h6vk x,t  h7mk x,t  2h12mkxx x,t , sk x,t0  k   sk1 x,t0  k  , mk x,t0  k   mk1 x,t0  k  ,  vk x,t0  k   vk1 x,t0  k  , Fk x,t0  k   Fk1 x,t0  k  , t0  k  t  t0  k 1 , k 1,2,... . ………………………………………………………………………………… (3) (4) (5) (6)

To ensure sufficient smoothness of the corresponding solutions at t  t0  , t  t0  2 , …, t  t0  n , in addition to the traditional smoothness conditions with respect to the initial conditions functions in the infectious disease model, it is necessary to impose conditions of their consistency at t  t0  , t  t0 , … [ 11 ]. In particular, the condition must be met

s0t x,t0  h3v0 x,t0   f 0 x,t0    h5 s0 x,t0 1  2h10s0xx x,t0.

Given that the diffusion redistributions of active factors are small compared to other components of the infectious disease process, we use the asymptotic method to solve the corresponding singularly perturbed model problems (4) - (5) [ 11,12 ]. In particular, the solutions of problems (4) - (5) are formally n n represented as asymptotic series v j (x,t)   ivij (x,t)  Rnvj x,t ,  , s j (x,t)   i sij (x,t)  Rnsj x,t ,  , io io f j (x,t)  n  i fij (x,t)  Rnfj x,t ,  , mj (x,t)  n  imij (x,t)  Rnmj x,t ,  as perturbations of the io io corresponding degenerate problems solutions [ 4 ], where j  0,1,...,k ,... , vij (x,t), sij (x,t) , fij (x,t) , mij (x,t) are the required functions (members of the asymptotics) Rnvj x,t ,  , Rnsj x,t ,  , Rnfj x,t ,  , Rnmj x,t ,  are relevant residual members. After substituting the asymptotic series and performing the standard procedure of "equalization", we obtain the following problems to determine the functions vij (x,t) , sij (x,t) , fij (x,t) , mij (x,t) i  0,1,...,n , j  0,1,...,k ,... : v0,0 t x,t  h  h2 f0,0 x,tv0,0 x,t  uV x,t ,

1 s0,0 t x,t  h3v0 x,t   f 0 x,t    h5 s0,0 x,t 1,  f0,0 t x,t  h4 s0,0 x,t  f0,0 x,t  h8v0,0 x,t f0,0 x,t  uF x,t,  m0,0 t x,t  h6v0,0 x,t  h7m0,0 x,t , s0,0 x,t0   s0 x , m0,0 x,t0   m0 x , v0,0 x,t0   v0 x,t0  , f0,0 x,t0   f 0 x,t0  , t0  t  t0  , v1,0 x,t0   0, f1,0 x,t0   0, t0  t  t0  , v1,0 t x,t  h1v1,0 x,t  h2 a0,0 x,t f1,0 x,t  b0,0 x,tv1,0 x,t  v1,0 x,t , s  1,0 t x,t  h3 a0,0 x,t   f1,0 x,t    b0,0 x,t   f1,0 x,t    h5s1,0 x,t ,  f1,0 t x,t  h4 s1,0 x,t  f1,0 x,t  h8 a0,0 x,t f1,0 x,t  b0,0 x,tv1,0 x,t   f 1,0 x,t ,  m1,0 t x,t  h6v1,0 x,t  h7m1,0 x,t , s  1,0 x,t0   0, m1,0 x,t0   0 ,

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi,0 x,t0   0, fi,0 x,t0   0, t0  t  t0  , vi,0 t x,t  h1vi,0 x,t  h2 a0,0 x,t fi,0 x,t  b0,0 x,tvi,0 x,t  V i,0 x,t, s  i,0 t x,t  h3 a0,0 x,t   fi,0 x,t    b0,0 x,t  vi,0 x,t    h5si,0 x,t  si,0 x,t,  fi,0 t x,t  h4 si,0 x,t  fi,0 x,t  h8 a0,0 x,t fi,0 x,t  b0,0 x,tvi,0 x,t   f i,0 x,t,  mi,0 t x,t  h6vi,0 x,t  h7mi,0 x,t  m i,0 x,t , s  i,0 x,t0   0, mi,0 x,t0   0 , (7) (8) (9) v0,k t x,t  h  h2 f0,k x,tv0,k x,t ,

1 s0,k t x,t  h3v0,k1 x,t   f0,k1 x,t    h5 s0,k x,t 1,  f0,k t x,t  h4 s0,k x,t  f0,k x,t  h8v0,k x,t f0,k x,t,  m0,k t x,t  h6v0,k x,t  h7m0,k x,t,  s0,k x,t0  k   s0,k1 x,t0  k  , m0,k x,t0  k   m0,k1 x,t0  k  , v0,k x,t0  k   v0,k1 x,t0  k  , f0,k x,t0  k   f0,k1 x,t0  k  ,  t0  k  t  t0  k 1 , v1,k t x,t  h1v1,k x,t  h2 a0,k x,t f1,k x,t  b0,k x,tv1,k x,t  v1,k x,t, s  1,k t x,t  h3 a0,k x,t   f1,k1 x,t    b0,k x,t  v1,k1 x,t    h5s1,k x,t,  f1,k t x,t  h4 s1,k x,t  f1,k x,t  h8 a0,k x,t f1,k x,t  b0,k x,tv1,k x,t   f 1,k x,t,  m1,k t x,t  h6V1,k x,t  h7m1,k x,t , s1,k x,t0  k   s1,k1 x,t0  k  , m1,k x,t0  k   m1,k1 x,t0  k , v1,k x,t0  k   v1,k1 x,t0  k  , f1,k x,t0  k   f1,k1 x,t0  k  ,  t0  k  t  t0  k 1 ,

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi,k t x,t  h1vi,k x,t  h2 a0,k x,t fi,k x,t  b0,k x,tvi,k x,t  vi,k x,t, s  i,k t x,t  h3 a0,k x,t   fi,k1 x,t    b0,k x,t  vi,k1 x,t    h5si,k x,t  si,k x,t,  fi,k t x,t  h4 si,k x,t  fi,k x,t  h8 a0,k x,t fi,k x,t  b0,k x,tvi,k x,t   f i,k x ,t ,  mi,k t x,t  h6vi,k x,t  h7mi,k x,t  m i,k x,t, si,k x,t0  k   si,n1 x,t0  k  , mi,k x,t0  k   mi,k1 x,t0  k , vi,k x,t0  k   vi,k1 x,t0  k  , fi,k x,t0  k   fi,k1 x,t0  k ,  t0  k  t  t0  k 1 , (10) (11) (12) where a0,j x,t  v0,j x,t, b0,j x,t  f0,j x,t ; 2v0,j x,t

x2 v1,j x,t  h9

,  f 1,j x,t  h11 i1 vi,j x,t  h2 vk ,j x,t fik ,j x,t  h9 l1 2 f0,j x,t

x2 2vi1,j x,t x2 ; , i1 s i,j x,t   h3vk ,j x,t   fik ,j x,t    h10 l1 2si2,j x,t x2

, m i,j x,t  h12 2mi2,j x,t x2

, i  2,3,...,n , j  0,1,...,k ,... .

Note that the proposed approach is easy to transfer to other, in particular, finite domains GZ . In this case, of course, instead of those described above, more complex schedules should be used (see, for example, [ 9,10,11 ]). Estimation of residual terms Rnvj x,t ,  , Rnsj x,t ,  , Rnfj x,t , , Rnmj x,t ,  and establishment of spatio-temporal intervals of convergence for forecasting of concrete processes are carried out on the basis of the principle of type of a maximum similarly to [ 4,11,12 ].

3. Numerical Experiments Results

Numerical experiments within this model investigated the features of the body's humoral immune response to viral antigens and the corresponding spatiotemporal dynamics of infectious disease for different situational conditions under diffusion disturbances and taking into account close to pulse changes in antigen and antibody concentrations in certain areas of the body.

Fig. 1 presents the spatial and temporal dynamics of antigen concentrations with the development of infectious disease in the chronic form according to model (1) - (2) in cases without taking into account (Fig. 1, a)) and taking into account (Fig. 1, b)) small spatially distributed diffusion influences under conditions that the initial distribution of antigens concentration is uneven in space vx,t0  v0 x e (x )2 (there is a separate center of infection of an organism with a maximum antigens concentration in a point x0  ). These results show that in the case without diffusion "redistribution", the development of the process according to the "scenario" of chronic disease is in some way "localized" in some area, which corresponds to the area with higher relative to some immunological barrier values of antigen concentration at the initial time. The influence of diffusion "redistribution" of the antigens initial concentration smooths out such "localization" of the model process. The corresponding model change with time of antigen concentration in the conditions of the chronic form of the disease at different intensities of diffusion "redistribution" in the infection epicenter is shown in Fig. 2. Under conditions without diffusion "redistribution" (ε = 0) the antigens v

v x a) t

t x b) concentration in the body increases to some maximum level, then decreases and is established over time at some stationary level. As the intensity of diffusion "redistribution" increases, the maximum value of antigen concentration in the epicenter decreases, and starting from some value of intensity the maximum antigens concentration will not increase, ie the level of immune protection adopted in the model in the presence of diffusion redistribution moment of time to prevent a model increase in the maximum antigens concentration in the infection epicenter and over time without "exacerbations" to reduce their concentration to some stationary level.

Spatial and temporal dynamics of antigen concentration with the development of infectious disease in a situation where at the initial time the values of the active factors of the process correspond to the values of the stationary solution, which characterizes the state of a healthy organism is shown in Fig. 3. The change in the concentration of antigens in the body is close to the pulse nature with the maximum value at some point xV at time tV . As already mentioned, the humoral type of immune response provides antibody neutralization of viral antigens that circulate freely in the blood or lymph. Depending on the immune system state, individual antigens can enter the cells of the target organ, where they can multiply and cause its destruction. As a result, a cell with a high antigens concentration appears at the site of the destroyed cell, which causes a close to impulse effect. The obtained results, as in the case of the initial condition with uneven distribution of antigen concentration, illustrate a certain "localization" in some area of the "scenario" of the disease in the chronic form (Fig. 3, a)) in the case without diffusion "redistribution". As in the previous case, the diffusion "redistribution" of the antigens concentration in the region of the impulse effect smooths the "localization" of the model process. Note that the presented generalization of the mathematical model of infectious disease taking into account the impulse effects under diffusion perturbations allows us to investigate the effects caused by several close to pulse sources of antigens with maximum values at different points xV j and in different time tV j .

Quite an effective procedure for the treatment of infectious diseases is the use of immunotherapy [ 6 ]. Donor antibodies can be administered by injection, which in this model will be presented as close to pulsed sources of donor antibodies with maximum values at points xF j in time tF j . Figure 4 presents the spatiotemporal dynamics of antigen concentrations with the development of infectious disease in the chronic form in the presence of close to pulsed sources of antigens and donor antibodies with maximum values at one point ( xV  xF ), but different time ( tV  tF ) under conditions excluding diffusion "redistribution" (Fig. 4, a)) and taking into account the diffusion "redistribution" (Fig. 4, b). The results illustrate, as expected, the decrease in antigen concentration due to the introduction of donor antibodies in the appropriate area. In conditions without diffusion "redistribution", the effect of donor antibodies in the cell of their introduction is longer. In terms of diffusion "redistribution" introduced donor antibodies over time "blur" to a larger "territory" of the body, resulting in faster consumption of antigen neutralization and their impact on the disease process is less long.

Fig. 5 presents the spatial and temporal dynamics of viral antigens concentrations with the development of infectious disease in the chronic form in the presence of one close to the pulse source of antigens and two close to pulse sources of donor antibodies with maximum values at point ( xV  xF ) and at different time ( tV  tF1  tF2 ), without taking into account (Fig. 5, a)) and taking into account the diffusion "redistribution" (Fig. 5, b)). As can be seen from the presented results, if the diffusion "redistribution" is taken into account, the intensity and duration of action of donor antibodies introduced into the body are smaller than in the model situation without such "redistribution" caused by diffusion "erosion" of donor antibodies from their injection site.

Figure 6 illustrates the change in the model dynamics of active infectious disease factors in the chronic form in the infection epicenter in a situation with one pulsed source of viral antigens and two pulsed sources of donor antibodies depending on the intensity of diffusion "redistribution" (parameter ε). The presented results show a decrease in the maximum value of the number of antigens, antibodies, plasma and affected cells in the infection epicenter with increasing intensity of diffusion "redistribution", which leads to a decrease in model "severity" of infectious disease. The introduction of close to the pulse of several sources of donor antibodies allows in this model to further reduce the antigens concentration in the infection epicenter. In particular, in a situation without taking into account the diffusion "redistribution", the effect of donor antibodies causes a decrease in the concentration of antigens to values close to zero. Given the diffusion effect, the action of donor antibodies due to their "redistribution" is less effective and leads to a smaller decrease in the antigens concentration. Therefore, to achieve the desired therapeutic effect, it is necessary to change the treatment procedure, for example, to increase the frequency of administration of antibodies, or their number in one injection.

4. Conclusions

Based on the modification of the simplest infectious disease model, an approach is presented to take into account close to impulse influences on the development of an infectious disease in the conditions of small spatially distributed diffusion perturbations. The corresponding model problem with delay is reduced to a sequence of problems without delay, for which representations of the required functions in the form of asymptotic series as perturbation of solutions of the corresponding degenerate problems are constructed.

The numerical experiments results illustrate the decrease in the maximum value of the antigens concentration in the infection epicenter due to their diffusion "redistribution", including for different situational conditions. It has been shown that even when the initial antigen concentration V 0 or the intensity of the pulsed antigen source in some area of the infection zone exceeds a certain critical value V * the diffusion "redistribution" over a period of time can reduce the critical antigen concentration to a level below the critical the reduction can be provided by the available level of antibodies, as well as a more economical mode of administration of donor antibodies by injection. That is, under this model, the "severity" of the infectious disease in such cases will decrease, so to speak, at low cost.

The developed computational procedure can be an element of designing specialized expert systems for making a wide range of decisions such as: can we in this case according to the values of relevant input data, in particular, on the intensity of diffusion "redistribution", on the level of immune protection available in the body, or, otherwise, to carry out external therapeutic effects. And in a situation where it is decided that such an effect should be exercised, in particular by means of injections of donor antibodies, to establish the most rational frequency of their introduction and an acceptable concentration of antibodies for each injection.

It is also promising to take into account such impulse effects in spatially distributed diffusion perturbations in the study of viral and bacterial diseases on the basis of more general models, in particular, models of antiviral and antibacterial immune response by Marchuk and Petrov [ 5 ].

5. References

References

[1] Marchuk G. I. Mathematical models in immunology , Moskow: Nauka , 1980 .

[2] Marchuk G. I. Mathematical models in immunology. Computational methods and experiments , Moskow: Nauka , 1991 .

[3] Marchuk G.L. Mathematical models of immune response in infectious diseases , Dordrecht: Kluwer Press, 1997 .

[4] Bomba A. Ya. Baranovsky S . V. Singular spatially distributed diffusion perturbations of one class of dynamic processes, Visnyk Natsionalnoho universytetu Vodnoho hospodarstva ta prirodokorystuvannia 3 (87) ( 2019 ): 54 - 65 .

[5] Bomba A. Ya. , Baranovsky S. V. , Pasichnyk O. V. , Pryschepa O. V. Modeling small -scale spatial distributed influences on the development of infectious disease process , Mathematical modeling and computing 2 ( 2020 ): 310 - 321 .

[6] Bomba A. Ya. , Baranovsky S. V. , Modeling of small spatially distributed influences on the dynamics of infectious disease in conditions such as pharmacotherapy . Journal of Computational and Applied Mathematics 1 ( 133 ) ( 2020 ): 5 - 17 .

[7] Nowak M.A. , May R. M. Virus dynamics . Mathematical principles of immunology and virology . Oxford University Press, 2000 .

[8] Wodarz D. Killer Cell Dynamics Mathematical and Computational Approaches to Immunology. Springer Science + Business Media, LLC, 2007 .

[9] Matt J. Keeling , Pejman Rohani . Modeling Infectious Diseases in Humans and Animals . Princeton University Press, 2008 .

[10] Elsholts L. E. , Norkin S . B. Introduction to the theory of differential equations with a deviating argument . Moskow: Nauka , 1971 .

[11] Bomba A. Ya. , Baranovsky S. V. , Prisiazniuk I. M. Nonlinear singularly perturbed problems of the “convection-diffusion” type . Rivne: NUWEE , 2008 .

[12] Bomba A. Ya . On the asymptotic method for solving one mass transfer problem during filtration in a porous medium , Ukrainian Mathematical Journal , Vol. 4 . 4. ( 1982 ): 493 - 496 .