Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Using Machine Learning in Predicting the Impact of Meteorological Parameters on Traffic Incidents
Next Article in Special Issue
Ergodicity and Related Bounds for One Particular Class of Markovian Time—Varying Queues with Heterogeneous Servers and Customer’s Impatience
Previous Article in Journal
A New Fuzzy Reinforcement Learning Method for Effective Chemotherapy
Previous Article in Special Issue
Double-Sources Queuing-Inventory Systems with Finite Waiting Room and Destructible Stocks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Control over the Distribution of Ticks Based on the Extensions of the KISS Model

by
Vassili N. Kolokoltsov
1,2,3
1
Faculty of Computation Mathematics and Cybernetics, Moscow State University, 119991 Moscow, Russia
2
Faculty of Economics, National Research University Higher School of Economics, 101000 Moscow, Russia
3
Department of Statistics, University of Warwick, Coventry CV4 7AL, UK
Mathematics 2023, 11(2), 478; https://doi.org/10.3390/math11020478
Submission received: 8 November 2022 / Revised: 7 January 2023 / Accepted: 9 January 2023 / Published: 16 January 2023

Abstract

:
Ticks and tick-borne diseases present a well-known threat to the health of people in many parts of the globe. The scientific literature devoted both to field observations and to modeling the propagation of ticks continues to grow. To date, the majority of the mathematical studies have been devoted to models based on ordinary differential equations, where spatial variability was taken into account by a discrete parameter. Only a few papers use spatially nontrivial diffusion models, and they are devoted mostly to spatially homogeneous equilibria. Here we develop diffusion models for the propagation of ticks stressing spatial heterogeneity. This allows us to assess the sizes of control zones that can be created (using various available techniques) to produce a patchy territory, on which ticks will be eventually eradicated. Using averaged parameters taken from various field observations we apply our theoretical results to the concrete cases of the lone star ticks of North America and of the taiga ticks of Russia. From the mathematical point of view, we give criteria for global stability of the vanishing solution to certain spatially heterogeneous birth and death processes with diffusion.

1. Introduction

Ticks and tick-borne diseases present a well-known threat to the health of people in many parts of the globe. The scientific literature devoted both to field observations and to modeling the propagation of ticks continues to grow.
To date, the majority of the mathematical studies were devoted to models based on ordinary differential equations, where spatial variability was taken into account by a discrete parameter. One of the most basic models was developed in [1], which is an ODE in four parameters (number of ticks V, number of hosts N, number of infected ticks X, number of infected hosts Y), and where the standard quadratic term of logistic models (expressing the internal competition of species) was replaced by a term limiting the growth of ticks by the number of available hosts. Apart from the analysis of equilibria (and their stability) performed for this model, Gaff et al. [1] also suggested a variant with multiple patches (with some fixed rates of transmission between them), however, with only a numerical analysis for this advanced variant.
In White et al. [2], this model was further developed to include the dynamics of two pathogens. Very complicated extensions (with various variables and parameters) were suggested and analyzed in [3,4,5]. In [6], a fractional version of tick-host dynamics was developed.
A stochastic version (Markov-chain model) of the model [1] was developed in [7], and the extinction predictions for deterministic and stochastic models were compared. In Gaff et al. [8] the models from a previous work [1] were enhanced by introducing some control parameters and the related optimal control policies (for quadratic costs) were derived aimed at reducing the propagation of diseases.
In two other studies [9,10], models with delays were suggested and analyzed via various methods (Lyapunov functions, Hopf bifurcations). A cellular-automaton model for the spatial propagation of ticks was developed and analyzed numerically in [11].
In Mwambi et al. [12], a stage-structured population model was proposed and analyzed. In the present paper we extend partially this model by including spatial dynamics.
Tosato et al. [13] investigated how different types of host targeted tick-control strategies (combined acaricide-repellent control) affect tick population and disease transmission, leading to some remarkable counterintuitive conclusions.
A reaction-diffusion model for the propagation of lyme disease was suggested [14] extending a similar spatially trivial model [15]. In a paper by Caraco et al. [14], the dynamics involved 9 density-variables: susceptible mice M, infected mice m, questing larvae L, larvae infesting susceptible mice V, larvae infesting pathogen infected mice v, susceptible questing nymphs N, questing infectious nymphs n, uninfected adult ticks A, pathogen infected adult ticks a. Position-independent equilibria were obtained and their local stability was investigated. This model was further developed by Zhao et al. [16], where the global stability of spatially trivial equilibria was analyzed for the case of bounded domains and the travel wave solutions were constructed for the case of unbounded domains.
The exploitation of partial differential equations for the analysis of spatial distribution of various species is a standard tool in ecology [17,18]. The latter paper, in particular, discussed the advantages and disadvantages of using the standard diffusion models, as compared with the models based on the telegraph equation. Recent references include papers [19,20] that deal with an attraction-repulsion chemotaxis model.
One of the standard tools is the so-called KISS model (from Kierstead, Slobodkin and Skellam, see [21,22]). It was developed initially in the framework of the plankton propagation. It describes the minimal size of a patch, where the population can survive under killing (mathematically Dirichlet) conditions on the boundary. For instance, in a one-dimensional case we look at the equation y ˙ = a y + λ y , with a , λ > 0 (see (1) below for all notations), on an interval [ 0 , R ] with the conditions y ( 0 ) = y ( R ) = 0 . The minimal R, for which a solution will not necessarily tend to zero, as time goes to infinity, is referred to as the critical patch size or the KISS scale. More generally one considers the behavior of the population when the interval [ 0 , R ] is surrounded by an infinite non-beneficial territory, where the dynamics are described by the equation y ˙ = a y μ y , with μ > 0 . The critical size then decreases depending on the ratio λ / μ . A review of these results and their various extensions can be found in [23], see also [24].
Let us mention a very important development of the KISS model for the description of several interacting species (though less relevant to the present paper), in particular diffusive versions of the Lotka–Volterra equations, which dealt with patchiness (spatially nontrivial stationary solutions) and its stability, see, e.g., [25,26,27].
The KISS model is the starting point for our analysis. However, eventually we are interested in creating control zones (kind of barriers) with non-beneficial conditions that can stop the propagation of ticks. We first extend the KISS model to the case when the non-beneficial region is bounded, interpreting its size as an additional control parameter, and the diffusion coefficient is discontinuous at the interface of beneficial and non-beneficial regions together with the mortality coefficient. Then we consider the modifications of the model with various boundary conditions, most importantly with the periodic boundary conditions. Seemingly not very relevant from the practical point of view, they provide the key intermediate step to tackle our main questions: how to organize the sizes and properties of the regular patches of (intermittent) control zones that would prevent the growth of the population of ticks. Our main results are partial answers to this question. We give first rather precise results for a homogeneous population, and then extend them partially to the structured populations, with several stages, such as larvae, nymphs, and adult ticks.
We suppose that in the control zone one is able to achieve the required level of the death rates or to decrease the reproduction rates of ticks. In the literature, one can find various methods that can be used locally for this purpose. It can be achieved by the treatment of the soil (planned burns, cleaning and decreasing humidity levels, chemical treatment, see, e.g., [13,28] and references therein), by tick-targeted strategies such as TickBots (see, e.g., [13,29]), by specific actions on hosts, such as dipping of cattle with acaricide (see, e.g., [30,31]) or by biological methods such as infesting woods with ants (see [32] for the related experiments on ants in the USA and [33] for the experiments in Russia), or natural enemies of ticks, such as the insects Ixodiphagus hookeri (see, e.g., [28]).
From the mathematical point of view, what we are doing here is obtaining criteria for negativity of certain differential operators, or more precisely, looking for the precise estimates to their highest eigenvalue. This is one of the basic problems in the theory of pseudo-differential equations. Many general results are available, which are often given; however, for the operators with smooth (at least major) symbol and/or up to certain constants, see, e.g., [34,35]. Here we solve this problem for a class of second order operators with piecewise continuous symbols. The negativity of the diffusion operators that govern the propagation of ticks ensures that the population will eventually die out (in the region considered).
The content of the paper is as follows. In Section 2, our models are introduced, and the main mathematical results are formulated. The first introductory subsection deals with a simple scalar model of one-dimensional diffusion. We give the necessary and sufficient conditions for the negativity of the corresponding diffusion operators (that imply the eradication of ticks) with various boundary conditions under the additional complication of discontinuous diffusion coefficients. These results should be of no surprise to specialists (they might be even known, but the author did not find a proper reference), because they are related to the famous Kronig-Penney model (though there the diffusion coefficient is a constant) of quantum mechanics. These results are tailored in a way to be used for concrete calculations on ticks given below.
The main results are obtained further. They concern the case of vector-valued diffusions, which allows one to analyse the evolution of the densities of ticks on various stages of their development. We first suggest some general, though rather rough, sufficient conditions for negativity by the method of symmetrization. Then we give some more precise conditions, but only in the case of two-stage populations and under rather restrictive assumptions on the control parameters.
In Section 3, we apply our theoretical results to the two examples of tick populations, the lone star tick, Amblyoma Americanum, from the south of USA, and the taiga tick, Ixodes persulcatus, the major type in Russian Siberian forests. There exist quite a lot of field investigations devoted to these species and thus an abundance of the observations, see, e.g., [30,36] and references therein for the American ticks and [28,37] for the taiga tick. Much less information is available for the European tick, Ixodes ricinus. However, as was previously noted in the mathematical modeling literature (see, e.g., [1]), there is no consensus on many parameters required for modelling, and the available parameters are varying by the types of the territories (upper or bottom land, upland wooded, meadow, etc), humidity and temperature, seasonal oscillation of hosts, etc. Therefore, numeric estimates of the parameters are obtained by applying some reasonable averages from the parameters found in the literature. Section 4 presents the graphs of certain simulations illustrating the decay of ticks’ populations under numeric data of the previous section. Section 5 is devoted to the proofs of our results. In the last section, some conclusions are drawn and further perspectives are indicated.

2. Models and Theoretical Results

2.1. Warm-Up: One-Stage Modeling

In this paper, we address the following question: can one place some barriers, organized by some control zones that are not beneficially for the multiplication of ticks, which would practically stop their propagation beyond these control zones?
Since we are talking about spatial propagation, our model must be necessarily spatially nontrivial. Therefore, we shall apply a diffusion model. As was mentioned above, a diffusion model was developed by Caraco et al. [14]. Unlike the model of this paper, our model here will be simpler, in the sense that we will stick to linear modeling (unlike the reaction-diffusion setting of [14]), but, on the other hand, it will include additional complication of being spatially heterogeneous, and precisely the effective control of these heterogeneities will be the main objective.
Remark 1. 
By not including in the model the quadratic (logistic type) terms of the standard nonlinear model, we are not including the competition between species, thus allowing ticks to multiply and propagate under better conditions than for nonlinear models. Consequently, if we manage to achieve the decay of their propagation in our linear case (what we are doing here), then the same decay will hold as well for the corresponding models with additional nonlinear terms limiting their growth. Moreover, since we are effectively giving stability criteria for the vanishing solution of a linear equation in a strong asymptotic sense (all eigenvalues have negative real part), adding additional quadratic terms can neither spoil this stability, nor improve it. The latter means that if there exists a mode with a positive eigenvalue of the linear part, then the vanishing solution will not be stable even in the presence of quadratic terms. This justifies the use of the linear model for our purposes in general term. Finally, the competition in ticks is performed mostly via available spaces on the bodies of hosts, so that the quadratic term in the equations for the dynamics of ticks is proportional to V ( M N V ) (see [1]), with the number of ticks V, the number of hosts N, and the maximum number M of ticks per host. Therefore, when hosts are in abundance, as compared to the basic equilibrium concentration of ticks (for instance, this occurs in many areas of Siberian taiga, see [28,37]), quadratic terms can be considered as overall negligible as compared to linear terms.
In order to not overcomplicate the story, we will not distinguish between infected and noninfected ticks, aiming at their total eradication (at least locally) and we shall not include explicitly the dynamics of hosts (mice, deer, etc.), their overall influence being reflected just by specifying certain average speed of motion of ticks, which will be effectively modeled as the diffusion coefficient. Moreover, though the natural state space for the dynamics of ticks is two-dimensional, we shall stick here to a one-dimensional modeling having in mind the propagation in a certain direction, with the second coordinate averaged out.
The ticks are known to have several stages in their development (larvae, nymphs, adult). To better explain our ideas, we shall first consider the model, where all generations are averaged out, and next describe the modifications arising when several distinct generations are explicitly included in the model.
Thus, we start with the simplest diffusion model for the dynamic of ticks in some area of their habitat (we denote the derivatives with respect to time and space by a dot and a prime, respectively):
y ˙ = a y + λ y , y ( 0 , x ) = y 0 ( x ) .
Here a > 0 is a diffusion coefficient reflecting the average speed of their random wandering and the coefficient λ > 0 reflects the speed of their multiplication (usually represented as the difference of average birth rates and average mortality). The density y = y ( t , x ) on some interval x [ 0 , R ] evolves according to (1) subject to additional boundary conditions, which can be taken as the Dirichlet conditions ( y ( t , 0 ) = y ( t , R ) = 0 ), or as the Neumann conditions ( y ( t , 0 ) = y ( t , R ) = 0 ) (or more general mixing Robin conditions that we shall not touch here), or rather periodic conditions ( y ( t , 0 ) = y ( t , R ) , y ( t , 0 ) = y ( t , R ) = 0 ).
Remark 2. 
The Dirichlet conditions arise in the situation, when ticks cannot survive at the boundary, which can be given, for instance, as a water reservoir, or as a temperature barrier. The Neumann conditions arise in the situation with reflecting barrier, such as a natural or artificial hedge. In other words, this is the situation, when no flux of matter crosses the border of the system either from outside or from inside. Periodic conditions do not seem to have a natural interpretation, but they are often convenient for the analysis (and therefore are abundantly used in physics), and can be also used as an intermediate step. In particular, the Dirichlet conditions are clearly the most disastrous for ticks, and therefore, if one can prove the dying out of ticks for periodic conditions, then the same will hold for the Dirichlet conditions. It is also intuitively clear (and again widely used in physics) that for large intervals the difference in boundary conditions does not affect the solutions in any essential way.
The main point in our modeling is the possibility to attach to the main region [ 0 , R ] (the background or the beneficial zone) a control zone of (much smaller length) r, which is not beneficial for ticks, meaning that their reproduction coefficient becomes negative there. One may also be able to control the diffusion coefficient shifting it to a different value b. Thus, instead of a homogeneous model (1) we consider a model with two distinctive regions:
y ˙ b e n = a y b e n + λ y b e n , x ( 0 , R ) y ˙ n b = b y n b μ y n b , x ( R , R + r ) ,
where positive μ , b , r are our control parameters. Here y b e n stands for the background beneficial region for ticks, which is not under our control, and y n b stands for the small non-beneficial region, where we are supposed to have control over at least some parameters. To make the problem well posed we should add, as above, some boundary conditions for x = 0 and x = R + r and moreover, the standard gluing conditions for diffusions with discontinuous coefficients (see, e.g., [38]):
y b e n ( R ) = y n b ( R ) , a y b e n ( R ) = b y n b ( R ) .
Let the operator L be defined on the functions of the interval [ 0 , R + r ] by the formula
L y ( x ) = a y b e n ( x ) + λ y b e n ( x ) , x ( 0 , R ) b y n b ( x ) μ y n b ( x ) , x ( R , R + r ) .
It is a standard fact from the theory of diffusion operators that, restricted to the functions satisfying the gluing condition (3) and either Dirichlet, or Neumann or periodic conditions, the operator L becomes self-adjoint with a discrete spectrum in the Hilbert space of the square integrable functions L 2 ( [ 0 , R + r ] ) on [ 0 , R + r ] . Moreover, this operator is bounded from above, so that there exist a decreasing sequence of eigenvalues l 0 l 1 l 2 and the corresponding orthonormal basis { ξ j } in L 2 ( [ 0 , R + r ] ) (depending on the chosen boundary conditions) such that L ξ j = l j ξ j for all j. This implies that for any initial function y 0 from L 2 ( [ 0 , R + r ] ) the solution to (2) with this initial condition is given by the convergent series
y ( x ) = j e t l j ξ j ( x ) 0 R + r y 0 ( z ) ξ j ( z ) d z .
Consequently, if all l j are negative, i.e., the operator L is negative, then this solution tends to zero for any initial conditions. Thus, the condition of eventual eradication of ticks in our model is the condition of negativity of L, or equivalently, the condition of the absence of nonnegative eigenvalues. In the rare case of vanishing maximal eigenvalues, the solution converges to a finite limit as t .
Theorem 1. 
Assume the Dirichlet boundary conditions. (i) If
λ / a π 2 / R 2 ,
L has a positive eigenvalue (independently of μ , r !).
(ii) If
π 2 / ( 2 R ) 2 < λ / a < π 2 / R 2 ,
the operator L is negative if and only if
tanh ( r μ / b ) b μ > tan ( R λ / a ) a λ .
(iii) If
λ / a π 2 / ( 2 R ) 2 ,
L has no nonnegative eigenvalues (and is negative in the case of the strict inequality).
Statement (i) is the well-known (initial) result of the KISS model yielding the critical patch size
R c = π a / λ ,
see [21,22]. Statement (ii) was proved for the limit r in [22].
Evidently, for the Neumann and periodic boundary condition the situation is different. There is no critical size, as the population can survive on any small interval. However, as the following results show, by introducing appropriate control zones, one can eradicate the population in a similar way to the above.
Theorem 2. 
Assume the Neumann boundary conditions. (i) If
λ / a π 2 / ( 2 R ) 2 ,
L has a positive eigenvalue (independently of μ , r !).
(ii) Otherwise, L is negative if and only if
μ b tanh ( r μ / b ) > λ a tan ( R λ / a ) ,
so that one can achieve negativity by choosing appropriate μ and r.
Theorem 3. 
Assume the periodic boundary conditions. (i) If (5) holds, then L has a positive eigenvalue. (ii) Otherwise, L is negative if and only if
μ b tanh ( r 2 μ / b ) > λ a tan ( R 2 λ / a ) .
One sees that Theorem 3 provides the same condition as Theorem 2 though with R / 2 and r / 2 instead of R and r. On the other hand, the boundary KISS value of (5) is the same for the Dirichlet and the periodic conditions. Moreover, one can check explicitly from conditions (8) and (6) that, as expected (see Remark 2), if the operator L is negative under the periodic condition, then it is also negative under the Dirichlet conditions.
Theorems 1 and 3 provide exact information on how μ , b , r can be tuned in order to achieve negativity of L and hence the eventual eradication of ticks. We see, in particular, that if R does not exceed certain critical level, one can eradicate the ticks by introducing a control zone of arbitrary small length r, if a sufficiently high level of mortality μ can be imposed on this control zone.
On the other hand, the theorems show that if λ R 2 / a is large enough, then no control zone can efficiently influence the global growth of the tick population. This observation leads to a natural idea that, in the case of a large territory, one can fight with the growth of ticks by arranging several small control zones placed in a periodic (patched) fashion. As we are going to show, if this is organized in a way that each pair of adjacent zones satisfies the conditions of the previous theorem, the ticks will be eradicated in the whole patched territory.
Thus, assume that our territory of habitat is represented as the interval [ 0 , K ( R + r ) ] , with some natural K, and that it is decomposed into 2 K subintervals, odd I 1 = [ 0 , R ] , I 3 = [ R + r , 2 R + r ] , ⋯, I 2 K 1 = [ ( K 1 ) ( R + r ) , K R + ( K 1 ) r ] , so that I 2 k 1 = [ ( k 1 ) ( R + r ) , k R + ( k 1 ) r ] for k = 1 , , K , and even I 2 = [ R , R + r ] , , I 2 K = [ K R + ( K 1 ) r , K ( R + r ) ] , so that I 2 k = [ k R + ( k 1 ) r , k ( R + r ) ] for k = 1 , , K . Assume that on the long odd intervals we have some background parameters λ > 0 , a > 0 , and on the short even intervals we have some (controllable) parameters b > 0 , μ > 0 , so that the diffusion is given by the first and second equations of (4) on even and odd intervals, respectively; see Figure 1.
Evidently, the gluing conditions (3) are supposed to hold on the interface of all intervals. Theorem 3 (valid for the case K = 1 ) can now be extended to the case of arbitrary K as follows.
Theorem 4. 
Assume the periodic boundary conditions for such diffusion on [ 0 , K ( R + r ) ] . Then for any K the conditions of Theorem 3 provide also the conditions for the corresponding diffusion operator L on [ 0 , K ( R + r ) ] to be negative or not.
Remark 3. 
As we already mentioned, if all positive solutions for the periodic conditions tend to zero, as t , then the same holds for the Dirichlet conditions. It is intuitively clear, as the Dirichlet conditions are less beneficial for survival. Formally it follows from the representations of solutions in terms of the Feynmann–Kac formula.

2.2. Several-Stage Modeling: Trivial Case

Until now, we have considered the situation with all stages of ticks’ lives averaged out. In more precise modeling, one has to take into account the presence of several stages. For ticks these are eggs, larvae, nymphs, and adults. With some reasonable averaging one can reduce the consideration of the lifespan of ticks to the two basic periods, from eggs to nymphs, and from nymphs to hatching female adults. On the other hand, a more detail analysis, can include not only the stages, but their time developments. Namely, say, nymphs can develop to adults in the same season as their own molting takes place, or after a diapause (wintering), so one can distinguish not only the stage, but also whether it develops in one year or two years.
Let us consider the general case of n stages (could be also generations for other species). The basic Equation (1) is generalized to the vector-valued equation
y ˙ = a y + M y , y ( 0 , x ) = y 0 ( x ) ,
where y R n and M is the birth-and-death matrix showing the progression of ticks from their birth through their various stages of development. This matrix M = M n has the standard form (used by many authors, see, e.g., [12] or [16]) reflecting the sequential propagation through various stages. They have elements
M j j = m j , j = 1 , , n , M j + 1 , j = b j , j = 1 , , n 1 , M 1 n = b n ,
with all other elements vanishing, where all parameters m j and b j are positive.
For instance, for dimensions 2 and 3, these matrixes have the form
M 2 = m 1 b 2 b 1 m 2 , M 3 = m 1 0 b 3 b 1 m 2 0 0 b 2 m 3 .
It is easy to see that a matrix M of type (10) always has a real eigenvalue.
Let us assume now that we can set a control zone with the increased death rates for ticks. Namely, let us consider the extension of Equation (2) of the form
y ˙ b e n = a y b e n + M y b e n , x ( 0 , R ) y ˙ n b = b y n b + ( M μ 1 ) y n b , x ( R , R + r ) ,
where 1 is the unit matrix. Notice that μ here has a slightly different meaning as in (2), referring to the relative decrease in the death rates. The corresponding diffusion operator becomes matrix-valued:
L y ( x ) = [ a y ( x ) + M y ( x ) ] 1 [ 0 , R ] ( x ) + [ b y ( x ) + ( M μ 1 ) ] 1 [ R , R + r ] ( x ) .
It turns out that the results of Theorems 1–4 have a straightforward extension to this case with the role of λ played by the largest eigenvalue of M.
Theorem 5. 
Assume that the maximal real eigenvalue Λ 1 of the matrix M in (12) is positive (if all eigenvalues of M have negative real parts, then the population would die out even in the background territory, the case of no interest to us) and all other eigenvalues Λ j have negative real part and are different (the latter conditions are technical simplifications that are not essential). Assume also that μ > Λ 1 (otherwise the ticks could survive even in the control zone alone).
Assume the Dirichlet boundary conditions for L. (i) If
Λ 1 / a π 2 / R 2 ,
L has a positive eigenvalue. (ii) If
π 2 / ( 2 R ) 2 < Λ 1 / a < π 2 / R 2 ,
the operator L is negative if and only if
tanh ( r ( μ Λ 1 ) / b ) b ( μ Λ 1 ) > tan ( R Λ 1 / a ) a Λ 1 .
(iii) If
Λ 1 / a π 2 / ( 2 R ) 2 ,
L has no nonnegative eigenvalues (and is negative in the case of the strict inequality).
Assume now that the territory of habitat is represented by the interval [ 0 , K ( R + r ) ] (with some natural K), decomposed into 2 K subintervals in the same way as formulated before Theorem 4, and that on the odd and even intervals our diffusion follows the first and the second equation of (12), respectively, with the usual gluing condition on the interfaces.
Theorem 6. 
Let the conditions of Theorem 5 for M hold and K be arbitrary. Assume that the periodic boundary conditions are chosen. (i) If (13) holds, then L has a positive eigenvalue. (ii) Otherwise the operator L is negative if and only if
( μ Λ 1 ) b tanh ( r 2 ( μ Λ 1 ) / b ) > Λ 1 a tan ( R 2 Λ 1 / a ) .

2.3. Several-Stage Modeling: Advanced Case

This section contains our main theoretical results.
Until now, we have looked at the case when diffusion coefficients differ in beneficial and control zones, but are independent of the stage. Evidently, usual averaging allows one to apply this model for ticks. However, for many types of ticks, larvae and nymphs use small rodents as hosts, while adult ticks use large mammals, such as deer, or birds, so that the displacement, and hence the diffusion coefficient differs drastically for the adults and the earlier stages of ticks. Hence it is more natural to use the model with different diffusions on different stages. Moreover, the matrixes M specifying the birth and death rates, can be of course quite different for the background and control zones, thus differing not only by a multiple of the unit matrix, as in (12).
Let us start with a simple extension of the KISS model. Namely, let us consider the extension of Equation (9) with variable diffusion, i.e., the equation
y ˙ = A y + M y , y ( 0 , x ) = y 0 ( x ) ,
where A is a diagonal matrix with diagonal elements a j > 0 , j = 1 , , n .
Theorem 7. 
Let M be a matrix of type (10). Suppose that the maximal eigenvalue Λ 1 = Λ ( M , A ) of the matrix A 1 M is strictly positive. Then the critical patch size equals R c = π / Λ 1 . That is, there exists a nontrivial solution to the equation A y + M y = E y with the Dirichlet boundary conditions and with some positive E if and only if R > R c .
The story becomes more complicated when we put together the original and a control zones. It seems difficult to expect here explicit necessary and sufficient conditions for the negativity of the spectrum, such as in the cases analyzed above. However, reasonable sufficient conditions can be obtained, as we are going to show now.
Let us consider the system
y ˙ b e n = A b e n y b e n + M b e n y b e n , y ( R / 2 , R / 2 ) , y ˙ n b = A n b y n b + M n b y n b , y ( R / 2 , r + R / 2 ) ,
with symmetric positive matrixes A b e n , A n b and arbitrary matrixes M b e n , M n b , where matrixes M n b , A n b are supposed to depend on some control parameters. We assume the periodic boundary conditions and the usual gluing conditions ( y n b = y b e n , A n b y n b = A b e n y b e n ) on the interface.
Equations for the eigenvalues are
A b e n y b e n + M b e n y b e n = E y b e n A n b y n b + M n b y n b = E y n b .
Theorem 8. 
Assume exactly k out of n eigenvalues λ 1 λ n of the symmetric matrix N b e n = ( M b e n A b e n 1 + A b e n 1 M b e n T ) / 2 are positive, and all eigenvalues μ 1 μ n of the symmetric matrix N n b = ( M n b A n b 1 + A n b 1 M n b T ) / 2 are negative (by T we denote the transposition). If
R R c s y m = π / λ 1 ,
and
| μ 1 | sinh ( r | μ 1 | ) 1 + cosh ( r | μ 1 | ) > 2 R n k λ 1 1 + cos ( R λ 1 ) ) ,
then system (20) has no solutions with non-negative E (and periodic boundary conditions).
This result extends automatically to the case of the territory represented by the interval [ 0 , K ( R + r ) ] , with some natural K, which is decomposed into 2 K subintervals in the same way as formulated before Theorem 4 (see also Theorem 6), so that on the odd and even intervals the diffusion follows the first and the second equation of (19).
Notice that the l.h.s. of (22) increases as | μ 1 | as in our previous results, but the r.h.s. is of order λ 1 (for R far away from R c s y m ), unlike λ 1 in Theorem 3.
In Theorem 7, the condition of the existence of a positive eigenvalue is given in terms of the maximal eigenvalue Λ 1 of the matrix A 1 M and in Theorem 8 a sufficient condition for negativity is linked with the maximal eigenvalue of the corresponding symmetrized matrix. Hence from the combination of these two theorems nothing can be said for sizes R such that
π / λ 1 = R c s y m < R < R c = π / Λ 1 .
In order to work with these sizes, the effective method of symmetrization used in Theorem 8 cannot be applied, and the analysis in arbitrary dimension seems to be quite involved. Working in dimension n = 2 , where explicit formulas for eigenvectors allow for a rather detailed analysis, possible results seem to depend strongly on the structure of diffusion coefficients and the matrix M. We present one such result valid for a range of coefficients that can be applied to tick populations.
Let us look at the equations
y ˙ b e n = A y b e n + M b e n y b e n y ˙ n b = a A y n b + M n b y n b
and the corresponding eigenvalue problem
A y b e n + M b e n y b e n = E y b e n a A y n b + M n b y n b = E y n b ,
with the usual gluing condition: y b e n coincides with y n b and y b e n coincides with a y n b on the interface of two regions.
Let n = 2 . For a E 0 , denote by Λ 1 ( E ) Λ 2 ( E ) the eigenvalues of the matrix N b e n ( E ) = A 1 M b e n A 1 E , and by μ 1 ( E ) μ 2 ( E ) the eigenvalues of the matrix N n b ( E ) = ( A 1 M n b A 1 E ) / a . Let Λ 1 = Λ 1 ( 0 ) > 0 > Λ 2 ( 0 ) and E 0 > 0 be such that Λ 1 ( E 0 ) = 0 . Assume that Λ 1 ( E ) 0 > Λ 2 ( E ) and 0 > μ 1 ( E ) > μ 2 ( E ) for all E [ 0 , E 0 ] . Let { v j ( E ) } and { w j ( E ) } be the corresponding bases of eigenvectors of N b e n ( E ) and N n b ( E ) , and c ( E ) = ( c i j ( E ) ) be the matrix that takes the basis { v j ( E ) } to the basis { w j ( E ) } .
Theorem 9. 
Assume the bases of eigenvectors can be chosen in such a way that
c 12 ( E ) c 21 ( E ) 0 , c 11 ( E ) c 22 ( E ) 0
for E [ 0 , E 0 ] . Then, if
a | μ 1 ( 0 ) | tanh ( | μ 1 ( 0 ) | r 2 ) > Λ 1 tan ( Λ 1 R 2 ) ,
there are no solutions to Equation (25) with periodic boundary conditions and non-negative E.
Comparing (27) with (8) we see that under (26) our two-dimensional condition of negativity is the exact extension of the one-dimensional case.
Remark 4. 
The technically convenient assumptions (26) (and especially the first of these two) are not very natural and should not hold in a general situation. However, as we show below, they can be achieved by an appropriate proportional change of the elements of the birth-and-death matrix.
Let us now look at a concrete situation when (26) holds. Namely, consider the system (24) where a = 1 , the matrix A is diagonal with the diagonal elements a 1 , a 2   M b e n has the form M 2 from (11) and M n b has the same form with m ˜ j and b ˜ j instead of m j and b j . Thus
N b e n ( 0 ) = A 1 M b e n = a 1 1 m 1 a 1 1 b 2 a 2 1 b 1 a 2 1 m 2 , N n b ( 0 ) = A 1 M n b = a 1 1 m ˜ 1 a 1 1 b ˜ 2 a 2 1 b ˜ 1 a 2 1 m ˜ 2 .
Assume that we can control the non-beneficial zone by a proportional decrease of the reproduction coefficients, i.e., calibrating
b ˜ j = ω b j , j = 1 , 2 ,
by choosing an appropriate parameter ω ( 0 , 1 ) and choosing m ˜ j m j in such a way that
m ˜ 1 m ˜ 2 m 1 m 2 .
For real ticks, a 1 is usually much less than a 2 (see numeric examples below). Hence the assumptions of the next theorem are fully relevant.
Theorem 10. 
Assume that det N b e n ( 0 ) = m 1 m 2 b 1 b 2 < 0 (this enures that Λ 1 ( 0 ) > 0 > Λ 2 ( 0 ) ) and
a 1 1 a 2 1 , a 1 1 m 1 a 2 1 m 2 .
Assume that b ˜ j and m ˜ j satisfy (28) and (29). Then the conditions of Theorem 9 (namely, inequalities (26)) hold, so that (27) is sufficient for the absence of nonnegative eigenvalues.

3. Numerical Results with Real Life Data

3.1. One-Stage Modeling with the North American Ticks

In US, the predominant types of ticks are the lone star tick (Amblyoma Americanum) in the south, the blacklegged tick (Ixodes scapularis, formerly called the deer tick) and the Americal dog tick (Dermacenter Variabilis) in the north. One can find lots of experimental research on the various parameters needed for modeling, which depend on many factors. Here, to estimate the birth and death rates, we will employ the averaged parameters for the lone star tick from [1]. Namely, the time unit is a month. Taking 2000 as the average number of eggs (occurring once in two year), 70 % survival rate and sex ratio 1:1, we obtain the birth rate of female larvae per month to be 0.5 × 0.7 × 2000 / 24 . With the out of host survival average 0.85 and the probability to find the host 0.03 (recall that without finding a host for a full meal blood no further development of a tick is possible) the total survival rate (of females) per month is calculated in [1] as
0.5 × 0.7 × 2000 / 24 × 0.03 × 0.85 = 0.75 .
Death rate is estimated as 0.01 in woods and as 0.1 in grass. Sticking to the case of grass, we obtain the parameter λ in Equation (1) to be λ = 0.75 0.1 = 0.65 (measured in month−1).
To estimate the diffusion coefficient a is a more difficult task.
We shall use the standard method, used both in physics and ecology (see, e.g., [17]), where the mean squared displacement M 2 ( t ) during a time t is estimated as
M 2 ( t ) = 2 a t .
The lifetime of the majority of lone star ticks is known to be 2 years, but sometimes it is over in one year. Thus, we can assume approximately that it has two rides per year. The main hosts for the lone star ticks (and many other American Ticks) are the white tailed deer that could travel for many kilometers during the 3–4 days needed for a tick to obtain its blood meal. One can choose about 10 km as the reasonable estimate for an average distance per a ride (see [31]). Note that the displacement due to ticks’ own movement is negligibly small compared to their displacement by the hosts.
Then the total squared displacement is M 2 ( 12 ) = ( 2 × 10 ) 2 = 400 , and therefore a = 400:24 ≈ 16.67, measured in km2 per month. Thus, for the critical patch size we obtain
R c = π 16.67 / 0.65 15.9 ( km ) .
Theorems 1–4 can be used to define the exact relation between the length of a control zones and the death rate in it, which is needed for the eradication of ticks.
For instance, let us apply Theorem 4.
Let us choose R = 14 (km) (which is reasonably close to the critical size R c ) and the length of the control zone r = 1 (km), and let b = a for simplicity. Then condition (8) for the eradication of ticks obtains the following numeric form
16.67 μ tanh ( 0.5 μ / 16.67 ) > 17.03 ,
so that the minimal required death rate μ is rather high, about 1960. These numbers are summarized in Figure 2.

3.2. One-Stage Modeling with the Taiga Ticks

Let us exploit here another time unit, choosing it to be one year. The life cycle of the taiga tick varies from 3 to 6, so that four years can be taken as an approximate average. Recall that in order to have a molting and to turn from one stage to another (larvae to a nymph, nymph to an adult, and finally to hatch) a tick must have a ride on a host with a full meal. Hence with 3–4 year cycle a tick makes on average a single ride per year.
Taking, as above, 10 km as an average distance per a meal (and thus per year) we obtain from (31) that a = 100 / 2 = 50 , measured in km2 per year.
As the mean survival rate we choose λ = 2 , measured in year−1, which is the approximate value used in [11], where it was shown to produce a reasonable fit to the experimental data available. Thus, for the critical patch size we obtain
R c = π 50 / 2 15.7 ( km ) ,
giving approximately the same result (!), as in the calculations above, based on the American experimental data.

3.3. Two-Stage Modeling with the Taiga Ticks

One of the ways to estimate the average travel distance of ticks on hosts can be obtained using chemical treatments of a controlled territory and looking for how far ticks can penetrate into the treated territory from the uncontrolled zone. These studies indicate (see Section VIII.3 of [28]) that, for a taiga tick, one can choose about 10 km as the average displacement for adult ticks during their meals (supporting the number used above and taken from the literature on American ticks) and about 1.5 km for nymphs and larvae. Therefore the diffusion coefficient a = 100 / 2 = 50 , used above for the whole population of ticks, in a more detailed analysis becomes the diffusion coefficient a 2 = 50 for the adult ticks only, while for all lower stages it can be estimated as only a 1 = 1 . 5 2 / 2 1.1 .
Let us work with the two-stage model, where the first stage represents the species grown up to the well-fed nymph, and the second stage represents adults. We shall use here the data for the taiga tick from the Sayans mountains of Siberia, Krasnoyarsk region, as presented in the very detailed observations of [37]. These observations show that only about 2.4 % of female tick obtain their meal and that about 1 / 2 of well-fed females produce eggs in the amount of about 5 thousand each. Taking also the standard sex ratio as 1:1, it follows that the potential number for the next generation is
p = 5000 × 0.024 × 0 . 5 2 = 30
per an adult tick. From this amount about 43 % survive producing hungry nymphs in the amount of 0.43 p = 12.9 . Only 8.2 % of the potential p survive to the stage of an adult well-fed nymphs (because of the high death rate on this stage), i.e., in total remain 0.082 p = 2.46 per an original adult tick.
Next, the death rate from a well-fed nymph to an adult tick is about 14.6 % , and the death rate during a winter is about 39 % . Thus, from a well-fed nymph one can expect about
0.854 × 0.61 0.52
hungry adults to appear next spring. Thus, the birth-and-death matrix M and the diffusion matrix A in (18) take the concrete form
M = 1 2.46 0.52 1 , A = 1.1 0 0 50 .
Therefore
A 1 M = 0.91 2.46 × 0.91 0.52 × 0.02 0.02 = 0.91 2.24 0.01 0.02 .
For the positive eigenvalue Λ 1 of this matrix we obtain approximately Λ 1 = 0.067 and thus the critical patch size of Theorem 7 equals
R c = π / 0.067 46.9 .
This value is essentially larger than the values (32) and (34), calculated for one-stage models and data (and taken from different sources). This can be expected. In fact, more detailed observations and calculations based on the regions of approximate equilibrium should reflect only slight possible average growth rates and hence weaker requirements for control zones. Theorem 10 can be used to assess the required decrease in the reproduction rates of control zones that would ensure the eradication.
Choosing, say, R = 40 , r = 1 , the condition of negativity (27) obtains the following numeric form:
μ 1 ( 0 ) tanh ( μ 1 ( 0 ) / 2 ) > 0.067 tan ( 1.34 ) = 0.29 .
Therefore | μ 1 | must be larger than about 0.64 . This condition yields the corresponding estimates for the required m ˜ j , b ˜ j . This setting is given schematically on Figure 3.
On the other hand,
1 2 ( A 1 M + M T A 1 ) = 0.91 1.125 1.125 0.02 .
Consequently, one obtains for the highest eigenvalue λ 1 of this symmetrized matrix that λ 1 0.86 . Therefore
R c s y m = π / 0.863 3.64 .
This value is much less than values (32) and (34). This corroborates the idea that the results of Theorem 8, which are very natural, theoretically, must be used cautiously in situations with highly nonsymmetric birth-and-death matrixes, as the symmetrization leads to strong distortions for such data.
Remark 5. 
One of the ways to enhance the practical application of Theorem 8 might be the introduction of additional life stages (increasing dimension) that would make the birth-and-death matrix more symmetric.

4. Simulations

Many insightful simulations (numeric solutions of Equations (2) and (19)) can be performed to support our findings and to observe the details of ticks’ population decay under various assumptions, boundary conditions and initial data. As illustrations we give three graphs related to real numeric data of the previous section. Everywhere the Dirichlet boundary conditions are assumed. The space and time coordinates are denoted by x and t, respectively.
The first graph, Figure 4 below, shows the behavior of the sum of the coordinates of the solution u 1 + u 2 for the situation given visually on Figure 3, in the case of just one pair of background and control zones. The initial function was taken as having both coordinates sin ( x / 45 ) .
The second graph, Figure 5 below, shows the behavior of the sum of the coordinates of the solution u 1 + u 2 for the situation given visually on Figure 3, in the case of three pairs of background and control zones. The process of dying out is seen to be slower on this much larger territory. Nevertheless the decay is seen to start rapidly after a short period of initial perturbations.
Finally we give a picture, Figure 6 below, illustrating the decay of ticks’ population in the simple model with all generations averaged out, which is visually presented on Figure 2. Again three pairs of background and control zones are used. The decay is seen to occur almost monotone and rather fast in this situation of rather poor complexity.

5. Proofs

Proof Theorem 1. 
Negativity of the operator L means that L has no non-negative eigenvalues. If such an eigenvalue E exists, then it solves the stationary problem L y = E y . Since the operator y y with Dirichlet boundary conditions is negative, it follows that E [ 0 , λ ) . Then the general solutions in the two domains are
y b e n = C sin λ E a x + D cos λ E a x , x ( 0 , R ) ,
y n b = A exp { μ + E b x } + B exp { μ + E b x } , x ( R , R + r ) .
The Dirichlet boundary condition y b e n ( 0 ) = 0 implies that D = 0 . Since eigenfunctions are defined up to a multiplicative constant, we can set C = 1 , so that
y b e n = sin λ E a x .
The Dirichlet boundary condition y n b ( R + r ) = 0 yields
B = A exp { 2 μ + E b ( R + r ) } .
By the gluing conditions y b e n ( R ) = y n b ( R ) and a y b e n ( R ) = b y n b ( R ) on the interface,
sin λ E a R = A exp { μ + E b R } + B exp { μ + E b R } ,
a ( λ E ) cos λ E a x = A b ( μ + E ) exp { μ + E b R } B b ( μ + E ) exp { μ + E b R } .
Thus
sin λ E a R = A exp { μ + E b R } exp { μ + E b ( R + 2 r ) } ,
a ( λ E ) cos λ E a x
= A b ( μ + E ) exp { μ + E b R } + exp { μ + E b ( R + 2 r ) } ,
or
tanh r μ + E b b ( μ + E ) = tan λ E a R a ( λ E ) .
Changing to x = ( λ E ) / a ( 0 , λ / a ] yields
tanh ( r ( λ + μ a x ) / b ) b ( λ + μ a x ) = tan ( R x ) a x .
The l.h.s. is negative decreasing on [ 0 , λ / a ] . Thus, if λ / a π 2 / ( 4 R 2 ) , there are no solutions x ( 0 , λ / a ] , because the r.h.s. is positive on this interval.
Moreover, the l.h.s. is decreasing on [ 0 , λ / a ] from
tanh ( r ( λ + μ ) / b ) b ( λ + μ ) to tanh ( r μ / b ) b μ .
Thus, if λ / a π 2 / R 2 , there is a solution x ( 0 , λ / a ] , because the l.h.s. necessarily intersects with the part of the r.h.s. increasing from to 0. Finally, if
π 2 4 R 2 < λ / a < π 2 R 2 ,
a solution exists if and only if the value of the l.h.s of (41) at x = λ / a lies below the graph of the r.h.s. yielding condition (6). □
Proof of Theorem 2. 
Similar calculations to Theorem 1 show that the condition for E [ 0 , λ ) to be an eigenvalue writes down as the equation
a ( λ E ) tan λ E a R = b ( μ + E ) tanh μ + E b r ,
or, in terms of x = ( λ E ) / a , as
( μ + λ a x ) b tanh r μ + λ a x b = a x tan ( R x ) .
The r.h.s. is positive increasing from 0 to on [ 0 , π 2 / 4 R 2 ) . The l.h.s. is positive decreasing on [ 0 , λ ] between two positive values (or between a positive value and 0 if μ = 0 ). Hence if
λ / a π 2 / ( 2 R ) 2 ,
there is a positive eigenvalue for any μ , r . Otherwise, a positive eigenvalue does not exist if and only (7) holds. □
Proof of Theorem 3. 
We again look at the system of equations
y ˙ b e n = a y b e n + λ y b e n , y ˙ n b = b y n b μ y n b .
It is convenient to choose coordinates so that the first equation holds for x ( R / 2 , R / 2 ) and the second for x ( R / 2 , R / 2 + r ) . Here the key point for the analysis is the symmetry. Namely, our system writes down as the equation y ˙ = L y on the torus y [ ( R + r ) / 2 , ( R + r ) / 2 ] with the identified right and left points. Alternatively one can think of y as periodic functions (with period R + r ) on the whole R . Here L is the linear operator such that L R = R L , where R is the reflection operator: R f ( x ) = f ( x ) . Hence if y is an eigenfunction of L, then so is also the function R y . However, L is an operator of the Schrödinger type, and it is well known from quantum mechanics (and actually very easy to prove, see, e.g., [39]) that the eigenvalues of one-dimensional Scrödinger operators are non-degenerate, i.e., there may exist only one eigenfunction to each eigenvalue, up to a multiplier. Hence R y = ± y . However, the minus sign would contradict the continuity and thus R y = y . Thus, an eigenfunction of L, which is the solution to the equation L Y = E y must be an even function. Similarly it has to be symmetric under the reflection with respect to the point ( R + r ) / 2 .
The equation L y = E y writes down more explicitly as the system
a y b e n + λ y b e n = E y b e n , b y n b μ y n b = E y n b ,
or equivalently
y b e n + λ E a y b e n = 0 , y n b μ + E b y n b = 0 .
As usual we are looking for the solutions with E [ 0 , λ ] . Then the functions y 1 and y 2 are of type (39) and (40) (now on intervals ( R / 2 , R / 2 ) , ( R / 2 , r + R / 2 ) ).
However, it is easy to show that if
A sin ( λ ( ξ + x ) ) + B cos ( λ ( ξ + x ) )
is an even function of x, then it is of the form C cos ( λ x ) (and the same holds for the linear combinations of sinh and cosh). Consequently, by the symmetry mentioned above, we can conclude that
y b e n = A cos ( λ E a x ) , y n b = B cosh μ + E b ( x ( R + r ) / 2 ) .
Gluing the functions and their derivatives on the interface x = R / 2 (recall (3)) yields
A cos ( λ E a R 2 ) = B cosh μ + E b r 2 ,
A a λ E a sin ( λ E a R 2 ) = B b μ + E b sinh μ + E b r 2 ,
and thus
a ( λ E ) tan λ E a R 2 = b ( μ + E ) tanh μ + E b r 2 ,
which is the same as (42) but with R / 2 and r / 2 instead of R ad r. Hence the proof is completed as in Theorem 2. □
Proof of Theorem 4. 
Now a solution y to the eigenvalue problem L y = E y is obtained by gluing solutions y o d d and y e v e n defined by formulas such as (39) and (40) on odd and even subintervals. The key point is, that under the periodic condition the operator L commutes with the shift operator T f ( x ) = f ( x + R + r ) , where the addition is understood modulo R + r . This means that for any eigenfunction y of L, the function T y is also an eigenfunction with the same eigenvalue. Using the nondegeneracy of eigenvalues of L (as in the proof of the previous theorem) we can conclude that T y = ω y for some ω . However, T k is the identity operator, so that ω k = 1 . Since our eigenfunctions are real it follows that ω = ± 1 . Again the multiplication by 1 would contradict the continuity and thus only ω = 1 is allowed. Hence we are directly in the setting of Theorem 3 and the proof is complete. □
Proof of Theorems 5 and 6. 
Under the conditions of the Theorems, there exists an invertible matrix C such that C M C 1 = D is diagonal with the elements Λ j on the diagonal. Then in the variables z = C y all equations and boundary conditions remain the same, but with M substituted by D. Hence our system decomposes into n independent equations for the coordinates z j of the vector-valued function z. For all z j with j > 1 there can be no positive eigenvalues, because all operators involved are strictly negative. Furthermore, the equations for z 1 (and all boundary conditions) become identical to the equations for y of the one-dimensional case, with Λ 1 instead of λ and the new μ being equal to λ + μ of the one-dimensional case. □
Proof of Theorem 7. 
The eigenvalue problem to Equation (18) with the Dirichlet boundary condition writes down as A y + M y = E y subject to y ( 0 ) = y ( R ) = 0 . Equivalently this equation rewrites as
y + ( A 1 M A 1 E ) y = 0 .
It is easy to see that the highest real eigenvalue of the matrix A 1 M A 1 E decreases with the increase of E. Hence, according to Theorem 5, Equation (44) is solvable for some positive E if and only if the highest real eigenvalue Λ 1 of the matrix A 1 M satisfies the condition Λ 1 > π 2 / R 2 implying the claim of the theorem. □
Proof of Theorem 8. 
Change the variable in (20) to z so that to z b e n = A b e n y , z n b = A n b y . Then we obtain
z b e n + ( M b e n E ) A b e n 1 z b e n = 0 z n b + ( M n b E ) A n b 1 z n b = 0 ,
and the gluing condition for z is the continuity of z and z on the interfaces.
Let L ( E ) denote the operator on the r.h.s. of (45), so that system (45) can be concisely written as L ( E ) z = 0 . We are interested in a criterion that ensures that this equation has no solution with E 0 . This would be the case, if one could show that ( z , L ( E ) z ) < 0 for all z. However, ( z , L ( E ) z ) = ( z , L S ( E ) z ) , where L S ( E ) = ( L ( E ) + L T ( E ) ) / 2 is the symmetrization. Thus, it is sufficient to show that L S ( E ) is a negative symmetric operator for any positive E. However, L S ( E ) decreases with the increase of E. Thus it is sufficient to show that L = L ( 0 ) is negative, i.e., L z = P z has no solutions with non-negative P.
Furthermore, this means that there are no solutions to the system
z b e n + ( M b e n S P ) z b e n = 0 z n b + ( M n b S P ) z n b = 0 ,
with P 0 and M b e n S = ( M b e n A b e n 1 + A b e n 1 M b e n T ) / 2 , M n b S = ( M n b A n b 1 + A n b 1 M n b T ) / 2 .
It is clear that the required sufficient condition is sufficient to show for the case when all eigenvalues of M b e n S and M n b S are different, as the general case is obtained by a straightforward limiting procedure.
For simplicity, let us assume that M b e n S has one positive eigenvalue (i.e., k = 1 ) and all other are negative: λ 1 > 0 λ 2 > > λ n .
Clearly P λ 1 are impossible in (46), so that only P [ 0 , λ 1 ] have to be analyzed.
Due to the periodic conditions, solutions to (46) are given by the formulas (see beginning of the proof of Theorem 3)
z b e n = A 1 cos ( λ 1 P x ) v 1 + j = 2 n A j cosh ( P λ j x ) v j ,
z n b = j = 1 n B j cosh ( P μ j ( x ( R + r ) / 2 ) ) w j ,
where v j and w j are orthonormal spectral bases for M b e n S and M n b S .
Let w j = k c j k v k with an orthogonal matrix C = ( c i j ) . Then
z n b = B k cosh ( P μ k ( x ( R + r ) / 2 ) ) c k j v j .
Thus, gluing z b e n and z n b at x = R / 2 yields
A 1 cos ( λ 1 P R / 2 ) = j B j cosh ( P μ j r / 2 ) c j 1 ,
and
A i cosh ( P λ j R / 2 ) = j B j cosh ( P μ j r / 2 ) c j i , i > 1 .
Next, since
R / 2 r + R / 2 ( z ( x ) , z ( x ) ) d x = R / 2 r + R / 2 ( z ( x ) , z ( x ) ) d x ,
which follows from the integration by parts and the gluing (and periodic) conditions, we have that, for a solution to (46),
0 b e n ( z b e n , ( M b e n S P ) z b e n ) d x + n b ( z n b , ( M n b S P ) z n b ) d x
= A 1 2 ( λ 1 P ) b e n cos 2 ( λ 1 P x ) d x + j > 1 A j 2 ( λ j P ) b e n cosh 2 ( P λ j x ) d x
+ j 1 B j 2 ( μ j P ) n b cosh 2 ( P μ j ( x ( R + r ) / 2 ) ) d x .
Since cos ( 2 x ) = 2 cos 2 x 1 and cosh ( 2 x ) = 2 cosh 2 x 1 ,
ξ ξ cos 2 ( α x ) d x = ξ + 1 2 ξ ξ cos ( 2 α x ) d x = ξ + sin ( 2 α ξ ) 2 α ,
and thus
A 1 2 ( λ 1 P ) R 2 + sin ( R λ 1 P ) 2 λ 1 P j 1 B j 2 ( P μ j ) r 2 + sinh ( r P μ j ) 2 P μ j .
In particular,
A 1 2 B j 2 P μ j λ 1 P r 2 + sinh ( r P μ j ) 2 P μ j R 2 + sin ( R λ 1 P ) 2 λ 1 P
for any j.
One could work with this estimate, but it makes the final result more transparent (though a bit rougher), if one simplifies this estimate to the following one:
A 1 2 B j 2 P μ j sinh ( r P μ j ) 2 R ( λ 1 P ) .
On the other hand,
| A 1 | cos ( λ 1 P R 2 ) j | B j | cosh ( P μ j r 2 ) c j 1 ,
and thus
A 1 2 cos 2 ( λ 1 P R 2 ) n j B j 2 cosh 2 ( P μ j r 2 ) c j 1 2 n max j B j 2 cosh 2 ( P μ j r 2 ) ,
or
A 1 2 ( 1 + cos ( λ 1 P R ) ) n max j B j 2 ( 1 + cosh ( P μ j r ) ) .
Thus, for a j realizing the maximum we have
n 1 + cosh ( P μ j r ) 1 + cos ( λ 1 P R ) P μ j sinh ( r P μ j ) 2 R ( λ 1 P ) ,
or
R ( λ 1 P ) ( 1 + cos ( R λ 1 P ) ) P μ j sinh ( r P μ j ) 2 n ( 1 + cosh ( r P μ j ) ) .
This can be rewritten as
1 R g 1 ( R λ 1 P ) 1 2 n r g 2 ( r P μ j )
with
g 1 ( x ) = x 2 1 + cos x , g 2 ( x ) = x sinh x 1 + cosh x .
In order to make (49) impossible (and thus to ensure the negativity of our operator) it is sufficient to assume that
1 R g 1 ( R x ) < 1 2 n r g 2 ( r λ 1 μ 1 x )
for x [ 0 , λ 1 ] , where x = λ 1 P . For λ 1 π 2 / R , this cannot be true. Hence assume λ 1 < π 2 / R .
Since it is straightforward to check that the function
g 1 ( R x ) g 2 ( r λ 1 μ 1 x )
is increasing for x [ 0 , λ 1 ] , condition (51) is equivalent to this condition on the endpoint:
1 R g 1 ( R λ 1 ) < 1 2 n r g 2 ( r | μ 1 | ) ,
or explicitly
2 R n λ 1 1 + cos ( R λ 1 ) ) < | μ 1 | sinh ( r | μ 1 | ) 1 + cosh ( r | μ 1 | ) ,
which coincides with (22) for k = 1 .
Modifications required for arbitrary k are straightforward, and we omit them. □
Proof of Theorem 9. 
The eigenvalue Equation (25) can be written as
y b e n + N b e n ( E ) y b e n = 0 y n b + N n b ( E ) y n b = 0 .
As above we can write
y b e n = A cos ( Λ 1 x ) v 1 + B cosh ( | Λ 2 | x ) v 2 ,
y n b = C cosh ( | μ 1 | ( x ( R + r ) / 2 ) ) w 1 + D cosh ( | μ 2 | ( x ( R + r ) / 2 ) ) w 2 .
The last equation is rewritten as
y n b = C cosh ( | μ 1 | ( x ( R + r ) / 2 ) ) ( c 11 v 1 + c 12 v 2 ) + D cosh ( | μ 2 | ( x ( R + r ) / 2 ) ) ( c 21 v 1 + c 22 v 2 ) .
Evidently, all terms here depend on E.
Furthermore, in order to shorten the formulas, we shall write sin, cos, tan if the argument is Λ 1 R / 2 and sinh, cosh, tanh if the argument is | Λ 2 | R / 2 . Gluing the solutions at the interface x = R yields
A cos = C cosh ( | μ 1 | r / 2 ) c 11 + D cosh ( | μ 2 | r / 2 ) c 21 ,
B cosh = C cosh ( | μ 1 | r / 2 ) c 12 + D cosh ( | μ 2 | r / 2 ) c 22 .
and
A Λ 1 sin = a ( C | μ 1 | sinh ( | μ 1 | r / 2 ) c 11 + D | μ 2 | sinh ( | μ 2 | r / 2 ) c 21 ) ,
B Λ 1 sinh = a ( C | μ 1 | sinh ( | μ 1 | r / 2 ) c 12 + D | μ 2 | sinh ( | μ 2 | r / 2 ) c 22 ) .
Equalising A from the 1st equations yields
Λ 1 tan [ C cosh ( | μ 1 | r / 2 ) c 11 + D cosh ( | μ 2 | r / 2 ) c 21 ]
= a [ C | μ 1 | sinh ( | μ 1 | r / 2 ) c 11 + D | μ 2 | sinh ( | μ 2 | r / 2 ) c 21 ] .
Equalising B from the 2nd equations yields
| Λ 2 | tanh [ C cosh ( | μ 1 | r / 2 ) c 12 + D cosh ( | μ 2 | r / 2 ) c 22 ]
= a [ C | μ 1 | sinh ( | μ 1 | r / 2 ) c 12 + D | μ 2 | sinh ( | μ 2 | r / 2 ) c 22 ] .
The condition for the existence of a solution to this system of two equations with two unknown C , D writes down as the equation
c 11 c 22 [ Λ 1 tan a | μ 1 | tanh ( | μ 1 | r / 2 ) ] [ | Λ 2 | tanh + a | μ 2 | tanh ( | μ 2 | r / 2 ) ]
c 12 c 21 [ Λ 1 tan a | μ 2 | tanh ( | μ 2 | r / 2 ) ] [ | Λ 2 | tanh + a | μ 1 | tanh ( | μ 1 | r / 2 ) ] = 0 ,
or as Λ 1 tan = a B / A with
A = c 11 c 22 [ | Λ 2 | tanh + a | μ 2 | tanh ( | μ 2 | r / 2 ) ]
c 12 c 21 [ | Λ 2 | tanh + a | μ 1 | tanh ( | μ 1 | r / 2 ) ] ,
B = c 11 c 22 | μ 1 | tanh ( | μ 1 | r / 2 ) [ | Λ 2 | tanh + a | μ 2 | tanh ( | μ 2 | r / 2 ) ]
c 12 c 21 | μ 2 | tanh ( | μ 2 | r / 2 ) ] [ | Λ 2 | tanh + a | μ 1 | tanh ( | μ 1 | r / 2 ) ] .
Since c 12 c 21 0 and c 11 c 22 0 (and | μ 2 | | μ 1 | ) , it follows that
B det c [ a | μ 2 | tanh ( | μ 2 | r / 2 ) + Λ 2 tanh ] | μ 1 | tanh ( | μ 1 | r / 2 ) ,
0 < A < det c [ | Λ 2 | tanh + a | μ 2 | tanh ( | μ 2 | r / 2 ) ] ,
so that
a B / A a | μ 1 | tanh ( | μ 1 | r / 2 ) .
Hence
Λ 1 tan ( λ 1 R / 2 ) a | μ 1 | tanh ( | μ 1 | r / 2 ) .
Consequently, if (27) holds, this is impossible. It remains to notice that, by monotonicity, this condition holds for all E if and only if it holds for E = 0 . □
Proof of Theorem 10. 
As is easily seen, if inequalities (30) hold for N b e n ( 0 ) and N n b ( 0 ) , then they hold for N b e n ( E ) and N n b ( E ) with all positive E. Hence it is sufficient to check condition (26) for E = 0 .
For a matrix
A = m 1 b 1 b 2 m 2 ,
the eigenvalues are
Λ 1 , 2 = 1 2 ( m 1 + m 2 ) ± ( m 1 m 2 ) 2 + 4 b 1 b 2 .
Thus, Λ 2 < 0 and Λ 1 > 0 .
Looking for eigenvectors v 1 , 2 , let us choose v 1 1 = 1 , v 2 2 = 1 . Then
v 1 2 = Λ 1 + m 1 b 1 = b 2 Λ 1 + m 2 ,
v 2 1 = b 1 Λ 2 + m 1 = Λ 2 + m 2 b 2 .
Thus
V = 1 v 2 1 v 1 2 1 , W = 1 w 2 1 w 1 2 1 .
Consequently,
C = W V 1 = 1 det V 1 v 2 1 w 1 2 w 2 1 v 2 1 v 1 2 + w 1 2 1 v 1 2 w 2 1 .
Assuming m 1 > m 2 , we see that
v 2 1 = ( m 1 m 2 ) + ( m 1 m 2 ) 2 + 4 b 1 b 2 2 b 2 < 0 ,
v 1 2 = ( m 1 m 2 ) + ( m 1 m 2 ) 2 + 4 b 1 b 2 2 b 1 < 0 .
Similarly, w 2 1 < 0 and w 1 2 > 0 . They have the same expression, but with m ˜ j and b ˜ j instead of m j , b j . Hence c 11 , c 22 1 > 0 . Moreover,
c 12 c 21 = ( w 1 2 v 1 2 ) ( | w 2 1 | | v 2 1 | ) .
Under (28) and (29), it follows that, for any ω ( 0 , 1 ) ,
w 1 2 v 1 2 , | w 2 1 | | v 2 1 | ,
and thus c 12 c 21 < 0 , which completes the proof. □

6. Conclusions

In this paper, we developed the diffusion model for the propagation of ticks with variable discontinuous diffusion coefficients. We extended the famous KISS model (originated in plankton research) to describe the dynamics on patchy territories with intermittent background (beneficial) and control (non-beneficial) zones and gave conditions for the control parameters that ensure the eradication of ticks’ population. Mathematically, these conditions imply the negativity of the corresponding diffusion operator. We started with rather trivial (though complete) results for the scalar case (one-stage populations) and then presented some advanced theorems for the vector-valued non-symmetric diffusions stressing serious theoretical difficulties that arise in this extension. Finally we used various published resources reporting on the concrete observations of ticks in order to extract some concrete numeric values for the parameters of background zones and thus to provide approximate numeric values to the required parameters of the control zones. In particular, we found that the roughest scalar model yields closer values of about 16 km for the KISS critical patch sizes related to both the lone star and the taiga ticks. On the other hand, using more detailed information combined with a more detailed multi-stage model can provide quite different output parameters. Furthermore, concrete numeric data showed that the estimates based on the symmetrization of the birth-and-death matrix (vary natural from the theoretical point of view) can underestimate strongly the lower bounds for the KISS scales.
Let us note some further perspectives. We used in this paper a one-dimensional diffusion model aiming at the approximate description of the propagation in a certain direction. Evidently, it would be natural to extend the results to a more realistic (for ticks and many other species) two-dimensional diffusion model. On the other hand, we discussed the solutions in a global closed territory. Of interest would be a more local analysis of the solutions in a (possibly intermittent) small territory placed in contact with a large reservoir of an expanding background population. Such solutions could be possibly analyzed based on the Fokas method, see [38,40], or via the method of estimating the heat content, see, e.g., [41]. It would be also of interest to see whether the results of Theorem 9 have any reasonable extensions beyond the restrictive assumptions (26) and whether one can obtain some multi-dimensional analogs. Furthermore, Remark 5 can be used for practical calculations with multi-dimensional models. A realistic vector-valued diffusion model can include not only ticks in various stages, but also the hosts.
Practically, the main goal of this paper was to give exact recommendations on the properties of control zones needed for the eradication of ticks. However, the numbers obtained in our analysis are based on very rough estimates, and our calculations show that they are quite sensitive to the input data. A much more detailed program of concrete observations in various places is needed to make the estimates more precise and more specific for concrete regions.
In this paper, we only used the diffusion models. Further investigations are needed to see whether more appropriate approximations can be achieved using other models, such as the telegraph equations, or presently very popular sub-diffusions or super-diffusions.
Please note finally that our models with periodic intermittent beneficial/non-beneficial zones present much simplified versions of the parabolic Anderson models with random potentials, to which an immense number of studies were devoted, see, e.g., [42] and references therein.

Funding

This work was supported by the Ministry of Science and Higher Education of the Russian Federation Grant ID: 075-15-2020-928.

Data Availability Statement

No data were used apart from those taken from published resources.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Gaff, H.D.; Gross, L.J. Modeling Tick-Borne Disease: A Metapopulation Model. Bull. Math. Biol. 2007, 69, 265–288. [Google Scholar] [CrossRef]
  2. White, A.; Schaefer, E.; Thompson, C.W.; Kribs, C.M.; Gaff, H. Dynamics of two pathogens in a single tick population. Lett. Biomath. 2019, 6, 50–66. [Google Scholar] [CrossRef]
  3. Nah, K.; Wu, J. Long-term transmission dynamics of tick-borne diseases involving seasonal variation and co-feeding transmission. J. Biol. Dyn. 2021, 15, 269–286. [Google Scholar] [CrossRef]
  4. Switkes, J.; Nannyonga, B.; Mugisha, J.Y.T.; Nakakawa, J. A mathematical model for Crimean-Congo haemorrhagic fever: Tick-borne dynamics with conferred host immunity. J. Biol. Dyn. 2016, 10, 59–70. [Google Scholar] [CrossRef] [Green Version]
  5. Zhang, X.; Sun, B.; Lou, Y. Dynamics of a periodic tick-borne disease model with co-feeding and multiple patches. J. Math. Biol. 2021, 82, 27. [Google Scholar] [CrossRef]
  6. dos Santos, J.P.C.; Cardoso, L.C.; Monteiro, E.; Lemes, N.H.T. A Fractional-Order Epidemic Model for Bovine Babesiosis Disease and Tick Populations. In Abstract and Applied Analysis; Hindawi: London, UK, 2015. [Google Scholar]
  7. Maliyoni, M.; Chirove, F.; Gaff, H.D.; Govinder, K.S. A Stochastic Tick-Borne Disease Model: Exploring the Probability of Pathogen Persistence. Bull. Math. Biol. 2017, 79, 1999–2021. [Google Scholar] [CrossRef]
  8. Gaff, H.D.; Schaefer, E.; Lenhart, S. Use of optimal control models to predict treatment time for managing tick-borne disease. J. Biol. Dyn. 2011, 5, 517–530. [Google Scholar] [CrossRef] [Green Version]
  9. Kashkynbayev, A.; Koptleuova, D. Global dynamics of tick-borne diseases. Math. Biosci. Eng. (MBE) 2020, 17, 4064–4079. [Google Scholar] [CrossRef]
  10. Shu, H.; Xu, W.; Wang, X.-S.; Wu, J. Complex dynamics in a delay differential equation with two delays in tick growth with diapause. J. Differ. Equ. 2020, 269, 10937–10963. [Google Scholar] [CrossRef]
  11. Vshivkova, O.A.; Komarov, A.S.; Frolov, P.V.; Khlebopros, R.G. The role of the heterogeneity of the habitat under the control of the ixode ticks number: Cellular-automaton model. Control Sci. 2013, 4, 57–63. (In Russian) [Google Scholar]
  12. Mwambi, H.G.; Baumgärtner, J.; Hadeler, K.P. Development of a stage-structured analytical population model for strategic decision making: The case of ticks and tick-borne diseases. Riv. Mat. Univ. Parma 2000, 3, 157–169. [Google Scholar]
  13. Tosato, M.; Nah, K.; Wu, J. Are host control strategies effective to eradicate tick-borne diseases (TBD)? J. Theor. Biol. 2021, 508, 110483. [Google Scholar] [CrossRef]
  14. Caraco, T.; Glavanakov, S.; Chen, G.; Flaherty, J.; Ohsumi, T.K.; Szymanski, B.K. Stage-structured infection transmission and a spatial epidemic: A model for Lyme disease. Am. Nat. 2002, 160, 348–359. [Google Scholar] [CrossRef]
  15. Caraco, T.; Gardener, G.; Maniatty, W.; Deelman, E.; Szymanski, B.K. Lyme Disease: Self regulation and Pathogen Invasion. J. Theor. Biol. 1998, 193, 561–575. [Google Scholar] [CrossRef]
  16. Zhao, X. Global dynamics of a reaction and diffusion model for Lyme disease. J. Math. Biol. 2012, 65, 787–808. [Google Scholar] [CrossRef]
  17. Holmes, E.E. Are diffusion models too simple? A comparison with telegraph models of invasion. Am. Nat. 1993, 142, 779–795. [Google Scholar] [CrossRef]
  18. Holmes, E.E.; Lewis, M.A.; Banks, J.E.; Veit, R.R. Partial differential equations in ecology: Spatial interactions and population dynamics. Ecology 1994, 75, 17–29. [Google Scholar] [CrossRef]
  19. Frassu, S.; Galván, R.R.; Viglialoro, G. Uniform in time L-estimates for an attraction-repulsion chemotaxis model with double saturation. Discret. Contin. Dyn. Syst. Ser. B 2022, 28, 1886–1904. [Google Scholar] [CrossRef]
  20. Frassu, S.; Li, T.; Viglialoro, G. Improvements and generalizations of results concerning attraction-repulsion chemotaxis models. Math. Methods Appl. Sci. 2022, 45, 11067–11078. [Google Scholar] [CrossRef]
  21. Kierstead, H.; Slobodkin, L.B. The Size of Water Masses Containing Plankton Bloom. J. Mar. Res. 1953, 12, 141–147. [Google Scholar]
  22. Skellam, J.G. Random Dispersal in theoretical biology. Bull. Math. Biol. 1991, 53, 135–165, Reprinted from Biometrika 1951, 38, 196–218. [Google Scholar] [CrossRef]
  23. Okubo, A. Critical patch size for plankton and patchiness. In Mathematical Ecology; Levin, S.A., Hallam, T.G., Eds.; Lecture Notes in Biomathematics 54; Springer: Berlin/Heidelberg, Germany, 1986; pp. 456–477. [Google Scholar]
  24. Murray, J.D. Mathematical Biology II: Spatial Models and Biomedical Applications, 3rd ed.; Springer: Berlin/Heidelberg, Germany, 2003. [Google Scholar]
  25. Baron, J.W.; Galla, T. Dispersal-induced instability in complex ecosystems. Nat. Commun. 2020, 11, 6032. [Google Scholar] [CrossRef]
  26. Gourley, S.A.; Britton, N.F.; Chaplain, M.A.J.; Byrne, H.M. Mechanisms for stabilisation and destabilisation of systems of reaction-diffusion equations. J. Math. Biol. 1996, 34, 857–877. [Google Scholar] [CrossRef]
  27. Mimura, M.; Murray, J.D. On a Diffusive Prey-Predator Model which Exhibits Patchiness. J. Theor. Biol. 1978, 75, 249–262. [Google Scholar] [CrossRef] [PubMed]
  28. Filippova, N.A. (Ed.) Taiga Tick Ixodes Persulcatus Schulze (Acarina, Ixodidae): Morphology, Systematics, Ecology, Medical Significance; Nauka: Moscow, Russia, 1985. (In Russian) [Google Scholar]
  29. Gaff, H.D.; White, A.; Leas, K.; Kelman, P.; Squire, J.C.; Livingston, D.L.; Sullivan, G.A.; Baker, E.W.; Sonenshine, D.E. TickBot: A novel robotic device for controlling tick populations in the natural environment. Ticks Tick-Borne Dis. 2015, 6, 146–151. [Google Scholar] [CrossRef] [PubMed]
  30. Mount, G.A.; Halle, D.G.; Barnard, D.R.; Daniels, E. New Version of LSTSIM for Computer Simulation of Amblyomma americanum (Acari: Ixodidae) Population Dynamics. J. Med. Entomol. 1993, 30, 843–857. [Google Scholar] [CrossRef] [PubMed]
  31. Ostfeld, R.S. The Ecology of Lyme-Disease Risk: Complex interactions between seemingly unconnected phenomena determine risks of exposure to this expanding disease. Am. Sci. 1997, 85, 338–346. [Google Scholar]
  32. Harris, G.W.; Burns, E.C. Predation on the lone star tick by the imported fire ant. Environ. Entomol. 1972, 1, 362–365. [Google Scholar] [CrossRef]
  33. Vshivkova, O.A. Mathematical simulation of ant influence on the ixode ticks number in Euroasian ecosystems. Izvestia Samarskogo Nauchnogo Cent. Russ. Acad. Sci. 2009, 11, 1631–1633. (In Russian) [Google Scholar]
  34. Fefferman, C.L. The uncertainty principle. Bulletin of the Americal Math. Society 1983, 9, 129–206. [Google Scholar]
  35. Sweezy, C. Relating different conditions for the positivity of the Schrödinger operator. Rocky Mt. J. Maths. 1993, 23, 353–366. [Google Scholar] [CrossRef]
  36. Mount, G.A.; Halle, D.G.; Daniels, E. Simulation of Blacklegged Tick (Acari: Ixodidae) Population Dynamics and Transmission of Borrelia burgdorferi. J. Med. Entomol. 1997, 34, 461–484. [Google Scholar] [CrossRef] [PubMed]
  37. Korotkov, U.S. Life cycle of the taiga tick Ixodes Persulcatus in the conifer forests of the bottom land of the East Sayan mountain range. Parasotologya 2014, 48, 20–36. (In Russuan) [Google Scholar]
  38. Sheils, N.E.; Deconinck, B. Initial-to-Interface Maps for the Heat Equation on Composite Domains. Stud. Appl. Math. 2016, 137, 140–154. [Google Scholar] [CrossRef] [Green Version]
  39. Landau, L.D.; Lifshitz, E.M. Quantum Mechanics, 2nd ed.; Pergamon Press: Oxford, UK, 1965. [Google Scholar]
  40. Fokas, A.S. A Unified Approach to Boundary Value Problems; SIAM: Philadelphia, PA, USA, 2008. [Google Scholar]
  41. Berg, M.v.; Gittins, K. Uniform bounds for the heat content of open sets in Euclidean space. Differ. Geom. Appl. 2015, 40, 67–85. [Google Scholar] [CrossRef]
  42. Makarova, Y.; Han, D.; Molchanov, S.; Yarovaya, E. Branching random walks with immigration. Lyapunov stability. Markov Process. Relat. Fields 2019, 25, 683–708. [Google Scholar]
Figure 1. General scheme of periodic control zones.
Figure 1. General scheme of periodic control zones.
Mathematics 11 00478 g001
Figure 2. Parameters of control zones for lone star ticks under a simple model.
Figure 2. Parameters of control zones for lone star ticks under a simple model.
Mathematics 11 00478 g002
Figure 3. Parameters of control zones for taiga ticks under a two-stage model.
Figure 3. Parameters of control zones for taiga ticks under a two-stage model.
Mathematics 11 00478 g003
Figure 4. Decay of ticks in the case of one pair of background/control zones.
Figure 4. Decay of ticks in the case of one pair of background/control zones.
Mathematics 11 00478 g004
Figure 5. Decay of ticks in the case of three pairs of background/control zones.
Figure 5. Decay of ticks in the case of three pairs of background/control zones.
Mathematics 11 00478 g005
Figure 6. Decay in the simplest one-stage model.
Figure 6. Decay in the simplest one-stage model.
Mathematics 11 00478 g006
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Kolokoltsov, V.N. On the Control over the Distribution of Ticks Based on the Extensions of the KISS Model. Mathematics 2023, 11, 478. https://doi.org/10.3390/math11020478

AMA Style

Kolokoltsov VN. On the Control over the Distribution of Ticks Based on the Extensions of the KISS Model. Mathematics. 2023; 11(2):478. https://doi.org/10.3390/math11020478

Chicago/Turabian Style

Kolokoltsov, Vassili N. 2023. "On the Control over the Distribution of Ticks Based on the Extensions of the KISS Model" Mathematics 11, no. 2: 478. https://doi.org/10.3390/math11020478

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop