# On the control over the distribution of ticks based on the extensions of the
KISS model

Vassili N. Kolokoltsov Department of Statistics, University of Warwick,
Coventry CV4 7AL UK, and Higher School of Economics Moscow, ORCID
0000-0001-6574-4615, Email<EMAIL_ADDRESS>
###### 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.
So far 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. Only 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.

Key words: propagation of ticks, ticks’ growth control, critical patch size,
KISS model, control zones, vector-valued diffusion models, discontinuous
diffusion coefficients.

Mathematics Subject Classification: 35K10, 35K40, 35P15, 35Q92, 92D25, 92D40.

##  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.

So far 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
[8], which is an ODE in 4 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, paper [8] also suggests a variant
with multiple patches (with some fixed rates of transmission between them),
however with only numerical analysis for this advanced variant.

In [38] this model was further developed to include the dynamics of two
pathogens. Very complicated extensions (with various variables and parameters)
were suggested and analysed in [33], [39] and [26]. In [4] a fractional
version of tick-host dynamics was developed.

A stochastic version (Markov-chain model) of the model from [8] was developed
in [20], and the extinction predictions for deterministic and stochastic
models were compared. In [9] the models from [8] 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 [15] and [30] models with delays were suggested and analysed via various
methods (Lyapunov functions, Hopf bifurcations). A cellular-automaton model
for the spatial propagation of ticks was developed and analysed numerically in
[37].

In [25] a stage structured population model was proposed and analysed. In the
present paper we extend partially this model by including spatial dynamics.

Paper [34] investigates 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
in [3] extending a similar spatially trivial model of [2]. In [3] 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 investigated. This model was further developed in [40], where the
global stability of spatially trivial equilibria was analysed 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, see e.g. [14]
and [13]. The latter paper, in particular, discusses the advantages and
disadvantages of using the standard diffusion models, as compared with the
models based on the telegraph equation. One of the standard tool is the so-
called KISS model (from Kierstead, Slobodkin and Skellam, see [16] and [31]).
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
the killing (mathematically Dirichlet) conditions on the boundary. For
instance, in one dimensional case we look at the equation
$\dot{y}=ay^{\prime\prime}+\lambda y$, with $a,\lambda>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 is described by the equation $\dot{y}=ay^{\prime\prime}-\mu y$, with
$\mu>0$. The critical size then decreases depending on the ratio
$\lambda/\mu$. A review of these results and their various extensions can be
found in [27], 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-Volterrra
equations, which dealt with patchiness (spatially nontrivial stationary
solutions) and its stability, see e.g. [1], [11], [21],

The KISS model is the starting point for our analysis. But 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. Being seemingly not much relevant from the practical
point of view, they provide the key intermediate step to tackle our main
questions: how to organise 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, like larvae,
nymphs and adults in 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
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. [34],
[6] and references therein), by tick-targeted strategies such as TickBots (see
e.g. [10] and [34]), by specific actions on hosts, like dipping of cattle with
acaricide (see e.g. [28] and [22]) or by biological methods like infesting
woods with ants (see [12] for the related experiments on ants in USA and [36]
for the experiments in Russia) or the natural enemies of ticks like the
insects Ixodiphagus hookeri (see e.g. [6]).

From the mathematical point of view, what we are doing here is getting
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. [5] and [32]. 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 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 ticks
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. [22] and [23] and
references therein for the American ticks and [6], [17] for the taiga tick.
Much less information is available for the European tick, Ixodes ricinus.
However, as was many times noted in the mathematical modelling literature (see
e.g. [8]), 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 is devoted to the proofs of our results. In the last section some
conclusions are drawn and further perspectives 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,
organised 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 in [3]. 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 [3]), 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 control their propagation in
our linear case (what we are doing here), then it will produce at least the
same effect for the corresponding models with additional nonlinear terms
limiting their growth.

In order to not overcomplicate the story, we will not distinguish infected and
noninfected ticks, aiming at their total eradication (at least locally) and we
shall not include explicitly the dynamics of hosts (mice, deers, 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 of their developments (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):

$\dot{y}=ay^{\prime\prime}+\lambda y,\quad y(0,x)=y_{0}(x).$ (1)

Here $a>0$ is a diffusion coefficient reflecting the average speed of their
random wandering and the coefficient $\lambda>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\in[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^{\prime}(t,0)=y^{\prime}(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^{\prime}(t,0)=y^{\prime}(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, like 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:

$\displaystyle\dot{y}_{ben}=ay^{\prime\prime}_{1}+\lambda y_{ben},\quad
x\in(0,R)$ (2) $\displaystyle\dot{y}_{nb}=by^{\prime\prime}_{nb}-\mu
y_{nb},\quad x\in(R,R+r),$

where positive $\mu,b,r$ are our control 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. [29]):

$y_{ben}(R)=y_{nb}(R),\quad ay^{\prime}_{ben}(R)=by^{\prime}_{nb}(R).$ (3)

Let the operator $L$ be define on the functions of the interval $[0,R+r]$ by
the formula

$Ly(x)=\left\\{\begin{aligned} &ay^{\prime\prime}_{ben}(x)+\lambda
y_{ben}(x),\quad x\in(0,R)\\\ &by^{\prime\prime}_{nb}(x)-\mu y_{nb}(x),\quad
x\in(R,R+r).\end{aligned}\right.$ (4)

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 exists a decreasing sequence of eigenvalues $l_{0}\geq l_{1}\geq
l_{2}\geq\cdots$ and the corresponding orthonormal basis $\\{\xi_{j}\\}$ in
$L_{2}([0,R+r])$ (depending on the chosen boundary conditions) such that
$L\xi_{j}=l_{j}\xi_{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)=\sum_{j}e^{tl_{j}}\xi_{j}(x)\int_{0}^{R+r}y_{0}(z)\xi_{j}(z)dz.$

Consequently, if all $\lambda_{j}$ are negative, that is 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 eigenvalue, the
solution converges to a finite limit as $t\to\infty$.

######  Theorem 2.1.

Assume the Dirichlet boundary conditions. (i) If

$\lambda/a\geq\pi^{2}/R^{2},$ (5)

$L$ has a positive eigenvalue (independently of $\mu,r$!).

(ii) If

$\pi^{2}/(2R)^{2}<\lambda/a<\pi^{2}/R^{2},$ (6)

the operator $L$ is negative if and only if

$-\frac{\tanh(r\sqrt{\mu/b})}{\sqrt{b\mu}}>\frac{\tan(R\sqrt{\lambda/a})}{\sqrt{a\lambda}}.$
(7)

(iii) If

$\lambda/a\leq\pi^{2}/(2R)^{2},$ (8)

$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}=\pi\sqrt{a/\lambda},$ (9)

see [16] and [31]. Statement (ii) was proved for the limit $r\to\infty$ in
[31].

Of course 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.2.

Assume the Neumann boundary conditions. (i) If

$\lambda/a\geq\pi^{2}/(2R)^{2},$ (10)

$L$ has a positive eigenvalue (independently of $\mu,r$!).

(ii) Otherwise, $L$ is negative if and only if

$\sqrt{\mu b}\tanh(r\sqrt{\mu/b})>\sqrt{\lambda a}\tan(R\sqrt{\lambda/a}),$
(11)

so that one can achieve negativity by choosing appropriate $\mu$ and $r$.

######  Theorem 2.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

$\sqrt{\mu b}\tanh(\frac{r}{2}\sqrt{\mu/b})>\sqrt{\lambda
a}\tan(\frac{R}{2}\sqrt{\lambda/a}).$ (12)

One sees that Theorem 2.3 provides the same condition as Theorem 2.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 (12) and (7) 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 2.1 and 2.3 provide exact information on how $\mu,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 $\mu$ can be
imposed on this control zone.

On the other hand, the theorems show that if $\lambda R^{2}/a$ is large
enough, then no control zone can efficiently influence the global growth of
ticks population. This observation leads to a natural idea that, in 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 organised 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 $2K$
subintervals, odd $I_{1}=[0,R]$, $I_{3}=[R+r,2R+r]$, $\cdots$,
$I_{2K-1}=[(K-1)(R+r),KR+(K-1)r]$, so that $I_{2k-1}=[(k-1)(R+r),kR+(k-1)r]$
for $k=1,\cdots,K$, and even $I_{2}=[R,R+r],\cdots,I_{2K}=[KR+(K-1)r,K(R+r)]$,
so that $I_{2k}=[kR+(k-1)r,k(R+r)]$ for $k=1,\cdots,K$. Assume that on the
long odd intervals we have some background parameters $\lambda>0,a>0$, and on
the short even intervals we have some (controllable) parameters $b>0,\mu>0$,
so that the diffusion is given by the first and second equations of (4) on
even and odd intervals, respectively. Of course the gluing conditions (3) are
supposed to hold on the interface of all intervals. Theorem 2.3 (valid for the
case $K=1$) can now be extended to the case of arbitrary $K$ as follows.

######  Theorem 2.4.

Assume the periodic boundary conditions for such diffusion on $[0,K(R+r)]$.
Then for any $K$ the conditions of Theorem 2.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\to\infty$, 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

So far 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 generalised to the vector-valued
equation

$\dot{y}=ay^{\prime\prime}+My,\quad y(0,x)=y_{0}(x),$ (13)

where $y\in\mathbf{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. [25] or [40]) reflecting the sequential propagation through
various stages. They have elements

$M_{jj}=-m_{j},\,j=1,\cdots,n,\quad M_{j+1,j}=b_{j},\,j=1,\cdots,n-1,\quad
M_{1n}=b_{n},$ (14)

with all other elements vanishing, where all parameters $\mu_{j}$ and $b_{j}$
are positive.

For instance, for dimensions $2$ and $3$, these matrices have the form

$M_{2}=\left(\begin{aligned} &-m_{1}\quad b_{2}\\\ &\,\,\,b_{1}\quad-
m_{2}\end{aligned}\right),\quad M_{3}=\left(\begin{aligned}
&-m_{1}\quad\,\,\,0\quad\quad b_{3}\\\ &\quad b_{1}\quad-m_{2}\quad\,\,0\\\
&\quad 0\quad\quad\,\,b_{2}\quad-m_{3}\end{aligned}\right).$ (15)

It is easy to see that a matrix $M$ of type (14) 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

$\displaystyle\dot{y}_{ben}=ay^{\prime\prime}_{ben}+My_{ben},\quad x\in(0,R)$
(16)
$\displaystyle\dot{y}_{nb}=by^{\prime\prime}_{nb}+(M-\mu\mathbf{1})y_{nb},\quad
x\in(R,R+r),$

where $\mathbf{1}$ is the unit matrix. Notice that $\mu$ 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:

$Ly(x)=[ay^{\prime\prime}(x)+My(x)]\mathbf{1}_{[0,R]}(x)+[by^{\prime\prime}(x)+(M-\mu\mathbf{1})]\mathbf{1}_{[R,R+r]}(x).$

It turns out that the results of Theorems 2.1-2.4 have a straightforward
extension to this case with the role of $\lambda$ played by the largest
eigenvalue of $M$.

######  Theorem 2.5.

Assume that the maximal real eigenvalue $\Lambda_{1}$ of the matrix $M$ in
(16) 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 $\Lambda_{j}$ have negative real
part and are different (the latter conditions are technical simplifications
that are not essential). Assume also that $\mu>\Lambda_{1}$ (otherwise the
ticks could survive even in the control zone alone).

Assume the Dirichlet boundary conditions for $L$. (i) If

$\Lambda_{1}/a\geq\pi^{2}/R^{2},$ (17)

$L$ has a positive eigenvalue. (ii) If

$\pi^{2}/(2R)^{2}<\Lambda_{1}/a<\pi^{2}/R^{2},$ (18)

the operator $L$ is negative if and only if

$-\frac{\tanh(r\sqrt{(\mu-\Lambda_{1})/b})}{\sqrt{b(\mu-\Lambda_{1})}}>\frac{\tan(R\sqrt{\Lambda_{1}/a})}{\sqrt{a\Lambda_{1}}}.$
(19)

(iii) If

$\Lambda_{1}/a\leq\pi^{2}/(2R)^{2},$ (20)

$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 $2K$ subintervals in the
same way as formulated before Theorem 2.4, and that on the odd and even
intervals our diffusion follows the first and the second equation of (16),
respectively, with the usual gluing condition on the interfaces.

######  Theorem 2.6.

Let the conditions of Theorem 2.5 for $M$ hold and $K$ be arbitrary. Assume
that the periodic boundary conditions are chosen. (i) If (17) holds, then $L$
has a positive eigenvalue. (ii) Otherwise the operator $L$ is negative if and
only if

$\sqrt{(\mu-\Lambda_{1})b}\tanh(\frac{r}{2}\sqrt{(\mu-\Lambda_{1})/b})>\sqrt{\Lambda_{1}a}\tan(\frac{R}{2}\sqrt{\Lambda_{1}/a}).$
(21)

###  2.3 Several-stage modeling: advanced case

This section contains our main theoretical results.

So far we have looked at the case when diffusion coefficients differ in
beneficial and control zones, but are independent of the stage. Of course,
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, like deers, or birds, so that the displacement, and
hence the diffusion coefficient differ 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 matrices $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 (16).

Let us start with a simple extension of the KISS model. Namely, let us
consider the extension of equation (13) with variable diffusion, that is, the
equation

$\dot{y}=Ay^{\prime\prime}+My,\quad y(0,x)=y_{0}(x),$ (22)

where $A$ is a diagonal matrix with diagonal elements $a_{j}>0$,
$j=1,\cdots,n$.

######  Theorem 2.7.

Let $M$ be a matrix of type (14). Suppose that the maximal eigenvalue
$\Lambda_{1}=\Lambda(M,A)$ of the matrix $A^{-1}M$ is strictly positive. Then
the critical patch size equals $R_{c}=\pi/\sqrt{\Lambda_{1}}$. That is, there
exists a nontrivial solution to the equation $Ay^{\prime\prime}+My=Ey$ 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, like in the cases,
analyzed above. However, reasonable sufficient conditions can be obtained, as
we are going to show now.

Let us consider the system

$\displaystyle\dot{y}_{ben}=A_{ben}y^{\prime\prime}_{ben}+M_{ben}y_{ben},\quad
y\in(-R/2,R/2),$ (23)
$\displaystyle\dot{y}_{nb}=A_{nb}y^{\prime\prime}_{nb}+M_{nb}y_{nb},\quad
y\in(R/2,r+R/2),$

with symmetric positive matrices $A_{ben}$, $A_{nb}$ and arbitrary matrices
$M_{ben}$, $M_{nb}$, where matrices $M_{nb}$, $A_{nb}$ are supposed to depend
on some control parameters. We assume the periodic boundary conditions and the
usual gluing conditions ($y_{nb}=y_{ben}$,
$A_{nb}y^{\prime}_{nb}=A_{ben}y^{\prime}_{ben}$) on the interface.

Equations for the eigenvalues are

$\displaystyle A_{ben}y^{\prime\prime}_{ben}+M_{ben}y_{ben}=Ey_{ben}$ (24)
$\displaystyle A_{nb}y^{\prime\prime}_{nb}+M_{nb}y_{nb}=Ey_{nb}.$

######  Theorem 2.8.

Let exactly $k$ out of $n$ eigenvalues $\lambda_{1}\geq\cdots\geq\lambda_{n}$
of the symmetric matrix
$N_{ben}=(M_{ben}A_{ben}^{-1}+A_{ben}^{-1}M_{ben}^{T})/2$ are positive, and
all eigenvalues $\mu_{1}\geq\cdots\geq\mu_{n}$ of the symmetric matrix
$N_{nb}=(M_{nb}A_{nb}^{-1}+A_{nb}^{-1}M_{nb}^{T})/2$ are negative (by $T$ we
denote the transposition). If

$R\leq R_{c}^{sym}=\pi/\sqrt{\lambda_{1}},$ (25)

and

$\frac{\sqrt{|\mu_{1}|}\sinh(r\sqrt{|\mu_{1}|})}{1+\cosh(r\sqrt{|\mu_{1}|})}>\frac{2Rnk\lambda_{1}}{1+\cos(R\sqrt{\lambda_{1}}))},$
(26)

then system (24) has no solutions with non-negative $E$.

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
$2K$ subintervals in the same way as formulated before Theorem 2.4 (see also
Theorem 2.6), so that on the odd and even intervals the diffusion follows the
first and the second equation of (23).

Notice that the r.h.s. of (26) increases as $\sqrt{|\mu_{1}|}$ like in our
previous results, but the l.h.s. is of order $\lambda_{1}$ (for $R$ far away
from $R_{c}^{sym}$), unlike $\sqrt{\lambda_{1}}$ in Theorem 2.3.

In Theorem 2.7 the condition of the existence of a positive eigenvalue is
given in terms of the maximal eigenvalue $\Lambda_{1}$ of the matrix $A^{-1}M$
and in Theorem 2.8 a sufficient condition for negativity is linked with the
maximal eigenvalue of the corresponding symmetrised matrix. Hence from the
combination of these two theorems nothing can be said for sizes $R$ such that

$\pi/\sqrt{\lambda_{1}}=R_{c}^{sym}<R<R_{c}=\pi/\sqrt{\Lambda_{1}}.$ (27)

In order to work with these sizes the effective method of symmetrization used
in Theorem 2.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

$\displaystyle\dot{y}_{ben}=Ay^{\prime\prime}_{ben}+M_{ben}y_{ben}$ (28)
$\displaystyle\dot{y}_{nb}=aAy^{\prime\prime}_{nb}+M_{nb}y_{nb}$

and the corresponding eigenvalue problem

$\displaystyle Ay^{\prime\prime}_{ben}+M_{ben}y_{ben}=Ey_{ben}$ (29)
$\displaystyle aAy^{\prime\prime}_{nb}+M_{nb}y_{nb}=Ey_{nb},$

with the usual gluing condition: $y_{ben}$ coincides with $y_{nb}$ and
$y^{\prime}_{ben}$ coincides with $ay^{\prime}_{nb}$ on the interface of two
regions.

Let $n=2$. For a $E\geq 0$, denote by $\Lambda_{1}(E)\geq\Lambda_{2}(E)$ the
eigenvalues of the matrix $N_{ben}(E)=A^{-1}M_{ben}-A^{-1}E$, and by
$\mu_{1}(E)\geq\mu_{2}(E)$ the eigenvalues of the matrix
$N_{nb}(E)=(A^{-1}M_{nb}-A^{-1}E)/a$. Let
$\Lambda_{1}=\Lambda_{1}(0)>0>\Lambda_{2}(0)$ and $E_{0}>0$ be such that
$\Lambda_{1}(E_{0})=0$. Assume that $\Lambda_{1}(E)\geq 0>\Lambda_{2}(E)$ and
$0>\mu_{1}(E)>\mu_{2}(E)$ for all $E\in[0,E_{0}]$. Let $\\{v_{j}(E)\\}$ and
$\\{w_{j}(E)\\}$ be the corresponding bases of eigenvectors of $N_{ben}(E)$
and $N_{nb}(E)$, and $c(E)=(c_{ij}(E))$ be the matrix that takes the basis
$\\{v_{j}(E)\\}$ to the basis $\\{w_{j}(E)\\}$.

######  Theorem 2.9.

Assume the bases of eigenvectors can be chosen in such a way that

$c_{12}(E)c_{21}(E)\leq 0,\quad c_{11}(E)c_{22}(E)\geq 0$ (30)

for $E\in[0,E_{0}]$. Then, if

$a\sqrt{\mu_{1}(0)}\tanh(\sqrt{\mu_{1}(0)}\frac{r}{2})>\sqrt{\Lambda}_{1}\tan(\sqrt{\Lambda_{1}}\frac{R}{2}),$
(31)

there are no positive solutions to equation (29).

Comparing (31) with (12) we see that under (30) our two-dimensional condition
of negativity is the exact extension of the one-dimensional case.

######  Remark 4.

The technically convenient assumptions (30) (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 (30) holds. Namely, consider the
system (28) where $a=1$, the matrix $A$ is diagonal with the diagonal elements
$a_{1},a_{2}$ $M_{ben}$ has the form $M_{2}$ from (15) and $M_{nb}$ has the
same form with $\tilde{m}_{j}$ and $\tilde{b}_{j}$ instead of $m_{j}$ and
$b_{j}$. Thus

$N_{ben}(0)=A^{-1}M_{ben}=\left(\begin{aligned} &-a_{1}^{-1}m_{1}\quad
a_{1}^{-1}b_{2}\\\ &\,\,\,a_{2}^{-1}b_{1}\quad-
a_{2}^{-1}m_{2}\end{aligned}\right),\quad
N_{nb}(0)=A^{-1}M_{nb}=\left(\begin{aligned} &-a_{1}^{-1}\tilde{m}_{1}\quad
a_{1}^{-1}\tilde{b}_{2}\\\ &\,\,\,a_{2}^{-1}\tilde{b}_{1}\quad-
a_{2}^{-1}\tilde{m}_{2}\end{aligned}\right).$

Assume that we can control the non-beneficial zone by a proportional decrease
of the reproduction coefficients, that is calibrating

$\tilde{b}_{j}=\omega b_{j},\quad j=1,2,$ (32)

by choosing an appropriate parameter $\omega\in(0,1)$ and choosing
$\tilde{m}_{j}\geq m_{j}$ in such a way that

$\tilde{m}_{1}-\tilde{m}_{2}\geq m_{1}-m_{2}.$ (33)

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 2.10.

Assume that $\det N_{ben}(0)=m_{1}m_{2}-b_{1}b_{2}<0$ (this enures that
$\Lambda_{1}(0)>0>\Lambda_{2}(0)$) and

$a_{1}^{-1}\geq a_{2}^{-1},\quad a_{1}^{-1}m_{1}\geq a_{2}^{-1}m_{2}.$ (34)

Assume that $\tilde{b}_{j}$ and $\tilde{m}_{j}$ satisfy (32) and (33). Then
the conditions of Theorem 2.9 (namely, equation (30)) hold, so that (31) is
sufficient for the negativity of all eigenvalues (and hence for the
eradication of ticks).

##  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 [8]. 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 get the birth rate of female larvae per month to be
$0.5\times 0.7\times 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 developments of ticks is possible) the total
survival rate (of females) per month is calculated in [8] as

$0.5\times 0.7\times 2000/24\times 0.03\times 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 get the parameter $\lambda$ in equation (1) to be
$\lambda=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.
[13]), where the mean squared displacement $M^{2}(t)$ during a time $t$ is
estimated as

$M^{2}(t)=2at.$ (35)

The life time of the majority of lone star ticks is known to be 2 years, but
sometimes it is completed 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 deers that could travel for many
kilometers during 3 - 4 days needed for a tick to get its blood meal. One can
choose about 10 km as the reasonable estimate for an average distance per a
ride (see [28]). 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\times 10)^{2}=400$, and
therefore $a=400:24\approx 16.67$, measured in km2 per month. Thus for the
critical patch size we get

$R_{c}=\pi\sqrt{16.67/0.65}\approx 15.9\,\text{(km)}.$ (36)

Theorems 2.1 \- 2.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 2.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 (12) for the eradication of ticks gets the
following numeric form

$\sqrt{16.67\mu}\tanh(0.5\sqrt{\mu/16.67})>17.03,$ (37)

so that the minimal required death rate $\mu$ is rather high, about 1958.

###  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 4 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 get from (35) that $a=100/2=50$, measured in km2 per year.

As the mean survival rate we choose $\lambda=2$, measured in year-1, which is
the approximate value used in [37], where it was shown to produce a reasonable
fit to the experimental data available. Thus for the critical patch size we
get

$R_{c}=\pi\sqrt{50/2}\approx 15.7\,\text{(km)},$ (38)

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 by using chemical treatments of a controlled territory and looking
for how far ticks can penetrate into the treated territory from the
uncontrolled zone. Theses studies indicate (see Section VIII.3 of [6]) 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\approx 1.1$.

Let us work with the two-stage model, where the first stage represent the
species grown up to the well-fed nymph, and the second stage represent 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
[17]. These observations show that only about $2.4\%$ of female tick get 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\times 0.024\times 0.5^{2}=30$

per an adult tick. From this amount about $43\%$ survive producing hungry
nymphs in the amount of $0.43p=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), that is in total remain $0.082p=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 deathrate during a winter is about $39\%$. Thus from a well-fed nymph
one can expect about

$0.854\times 0.61\approx 0.52$

hungry adults to appear next spring. Thus the birth-and-death matrix $M$ and
the diffusion matrix $A$ in (22) take the concrete form

$M=\left(\begin{aligned} &-1\quad 2.46\\\
&0.52\quad-1\end{aligned}\right),\quad A=\left(\begin{aligned} &1.1\quad 0\\\
&\,\,0\quad 50\end{aligned}\right).$ (39)

Therefore

$A^{-1}M=\left(\begin{aligned} &-0.91\quad 2.46\times 0.91\\\ &0.52\times
0.02\quad-0.02\end{aligned}\right)=\left(\begin{aligned} &-0.91\quad 2.24\\\
&\,\,\,0.01\quad-0.02\end{aligned}\right)$

For the positive eigenvalue $\Lambda_{1}$ of this matrix we get approximately
that $\sqrt{\Lambda_{1}}=0.067$ and thus the critical patch size of Theorem
2.7 equals

$R_{c}=\pi/0.067\approx 46.9.$ (40)

This value is essentially larger than the values (36) and (38), 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 2.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 (31) gets the following
numeric form:

$\sqrt{\mu_{1}(0)}\tanh(\sqrt{\mu_{1}(0)}/2)>0.067\tan(1.34)=0.29.$ (41)

Therefore $|\mu_{1}|$ must be larger than about $0.64$. This condition yields
the corresponding estimates for the required $\tilde{m}_{j},\tilde{b}_{j}$.

On the other hand,

$\frac{1}{2}(A^{-1}M+M^{T}A^{-1})=\left(\begin{aligned} &-0.91\quad 1.125\\\
&\,\,\,1.125\quad-0.02\end{aligned}\right).$

Consequently, one gets for the highest eigenvalue $\lambda_{1}$ of this
symmetrised matrix that $\sqrt{\lambda}_{1}\approx 0.86$. Therefore

$R_{c}^{sym}=\pi/0.863\approx 3.64.$ (42)

This value is much less than values (36) and (38). This corroborates the idea
that the results of Theorem 2.8, which are very natural theoretically, must be
used cautiously in situations with highly nonsymmetric birth-and-death
matrices, as the symmetrization leads to strong distortions for such data.

######  Remark 5.

One of the ways to enhance the practical application of Theorem 2.8 might be
the introduction of additional life stages (increasing dimension) that would
make the birth-and-death matrix more symmetric.

##  4 Proofs

Proof of Theorem 2.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
$Ly=Ey$. Since the operator $y\to y^{\prime\prime}$ with Dirichlet boundary
conditions is negative, it follows that $E\in[0,\lambda)$. Then the general
solutions in the two domains are

$y_{ben}=C\sin\left(\sqrt{\frac{\lambda-E}{a}}x\right)+D\cos\left(\sqrt{\frac{\lambda-E}{a}}x\right),\quad
x\in(0,R),$ (43)
$y_{nb}=A\exp\\{\sqrt{\frac{\mu+E}{b}}x\\}+B\exp\\{-\sqrt{\frac{\mu+E}{b}}x\\},\quad
x\in(R,R+r).$ (44)

The Dirichlet boundary condition $y_{ben}(0)=0$ implies that $D=0$. Since
eigenfunctions are defined up to a multiplicative constant, we can set $C=1$,
so that

$y_{ben}=\sin\left(\sqrt{\frac{\lambda-E}{a}}x\right).$

The Dirichlet boundary condition $y_{nb}(R+r)=0$ yields

$B=-A\exp\\{2\sqrt{\frac{\mu+E}{b}}(R+r)\\}.$

By the gluing conditions $y_{ben}(R)=y_{nb}(R)$ and
$ay^{\prime}_{ben}(R)=by^{\prime}_{nb}(R)$ on the interface,

$\sin\left(\sqrt{\frac{\lambda-E}{a}}R\right)=A\exp\\{\sqrt{\frac{\mu+E}{b}}R\\}+B\exp\\{-\sqrt{\frac{\mu+E}{b}}R\\},$
$\sqrt{a(\lambda-E)}\cos\left(\sqrt{\frac{\lambda-E}{a}}x\right)=A\sqrt{b(\mu+E)}\exp\\{\sqrt{\mu+E}R\\}-B\sqrt{b(\mu+E)}\exp\\{-\sqrt{\mu+E}R\\}.$

Thus

$\sin\left(\sqrt{\frac{\lambda-E}{a}}R\right)=A\left[\exp\\{\sqrt{\frac{\mu+E}{b}}R\\}-\exp\\{\sqrt{\frac{\mu+E}{b}}(R+2r)\\}\right],$
$\sqrt{a(\lambda-E)}\cos\left(\sqrt{\frac{\lambda-E}{a}}x\right)$
$=A\sqrt{b(\mu+E)}\left[\exp\\{\sqrt{\frac{\mu+E}{b}}R\\}+\exp\\{\sqrt{\frac{\mu+E}{b}}(R+2r)\\}\right],$

or

$-\frac{\tanh\left(r\sqrt{\frac{\mu+E}{b}}\right)}{\sqrt{b(\mu+E)}}=\frac{\tan\left(\sqrt{\frac{\lambda-E}{a}}R\right)}{\sqrt{a(\lambda-E)}}.$

Changing to $x=(\lambda-E)/a\in(0,\lambda/a]$ yields

$-\frac{\tanh(r\sqrt{(\lambda+\mu-ax)/b})}{\sqrt{b(\lambda+\mu-
ax)}}=\frac{\tan(R\sqrt{x})}{a\sqrt{x}}.$ (45)

The l.h.s. is negative decreasing on $[0,\lambda/a]$. Thus, if
$\lambda/a\leq\pi^{2}/(4R^{2})$, there is no solutions $x\in(0,\lambda/a]$,
because the r.h.s. is positive on this interval.

Moreover, the l.h.s. is decreasing on $[0,\lambda/a]$ from

$-\frac{\tanh(r\sqrt{(\lambda+\mu)/b})}{\sqrt{b(\lambda+\mu)}}\,\,\text{to}\,\,-\frac{\tanh(r\sqrt{\mu/b})}{\sqrt{b\mu}}.$

Thus, if $\lambda/a\geq\pi^{2}/R^{2}$, there is a solution
$x\in(0,\lambda/a]$, because the l.h.s. necessarily intersects with the part
of the r.h.s. increasing from $-\infty$ to $0$. Finally, if

$\frac{\pi^{2}}{4R^{2}}<\lambda/a<\frac{\pi^{2}}{R^{2}},$

a solution exists if and only if the value of the l.h.s of (45) at
$x=\lambda/a$ lies below the graph of the r.h.s. yielding condition (7).

Proof of Theorem 2.2. Similar calculations to Theorem 2.1 show that the
condition for $E\in[0,\lambda)$ to be an eigenvalue writes down as the
equation

$\sqrt{a(\lambda-E)}\tan\left(\sqrt{\frac{\lambda-E}{a}}R\right)=\sqrt{b(\mu+E)}\tanh\left(\sqrt{\frac{\mu+E}{b}}r\right),$
(46)

or, in terms of $x=(\lambda-E)/a$, as

$\sqrt{(\mu+\lambda-ax)b}\tanh\left(r\sqrt{\frac{\mu+\lambda-
ax}{b}}\right)=a\sqrt{x}\tan(R\sqrt{x}).$ (47)

The r.h.s. is positive increasing from $0$ to $\infty$ on
$[0,\pi^{2}/4R^{2})$. The l.h.s. is positive decreasing on $[0,\lambda]$
between two positive values (or between a positive value and $0$ if $\mu=0$).
Hence if

$\lambda/a\geq\pi^{2}/(2R)^{2},$

there is a positive eigenvalue for any $\mu,r$. Otherwise, a positive
eigenvalue does not exist if and only (11) holds.

Proof of Theorem 2.3. We again look at the system of equations

$\dot{y}_{ben}=ay^{\prime\prime}_{ben}+\lambda
y_{ben},\quad\dot{y}_{nb}=by^{\prime\prime}_{nb}-\mu y_{nb}.$

It is convenient to choose coordinates so that the first equation holds for
$x\in(-R/2,R/2)$ and the second for $x\in(R/2,R/2+r)$. Here the key point for
the analysis is the symmetry. Namely, our system writes down as the equation
$\dot{y}=Ly$ on the torus $y\in[-(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 $\mathbf{R}$. Here $L$ is the linear operator
such that $L\mathcal{R}=\mathcal{R}L$, where $\mathcal{R}$ is the reflection
operator: $\mathcal{R}f(x)=f(-x)$. Hence if $y$ is an eigenfunction of $L$,
then so is also the function $\mathcal{R}y$. But $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. [18]) that the eigenvalues of one-dimensional
Scrödinger operators are non-degenerate, that is, there may exist only one
eigenfunction to each eigenvalue, up to a multiplier. Hence $\mathcal{R}y=\pm
y$. But the minus sign would contradict the continuity and thus
$\mathcal{R}y=y$. Thus an eigenfunction of $L$, which is the solution to the
equation $LY=Ey$ must be an even function. Similarly it has to be symmetric
under the reflection with respect to the point $(R+r)/2$.

The equation $Ly=Ey$ writes down more explicitly as the system

$ay^{\prime\prime}_{ben}+\lambda y_{ben}=Ey_{ben},\quad
by^{\prime\prime}_{nb}-\mu y_{nb}=Ey_{nb},$

or equivalently

$y^{\prime\prime}_{ben}+\frac{\lambda-E}{a}y_{ben}=0,\quad
y^{\prime\prime}_{nb}-\frac{\mu+E}{b}y_{nb}=0.$

As usual we are looking for the solutions with $E\in[0,\lambda]$. Then the
functions $y_{1}$ and $y_{2}$ are of type (43) and (44) (now on intervals
$(-R/2,R/2)$, $(-R/2,r+R/2)$).

But it is easy to show that if

$A\sin(\sqrt{\lambda}(\xi+x))+B\cos(\sqrt{\lambda}(\xi+x))$

is an even function of $x$, then it is of the form $C\cos(\sqrt{\lambda}x)$
(and the same holds for the linear combinations of $\sinh$ and $\cosh$).
Consequently, by the symmetry mentioned above, we can conclude that

$y_{ben}=A\cos(\sqrt{\frac{\lambda-E}{a}}x),\quad
y_{nb}=B\cosh\left(\sqrt{\frac{\mu+E}{b}}(x-(R+r)/2)\right).$

Gluing the functions and their derivatives on the interface $x=R/2$ (recall
(3)) yields

$A\cos(\sqrt{\frac{\lambda-E}{a}}\,\frac{R}{2})=B\cosh\left(\sqrt{\frac{\mu+E}{b}}\,\frac{r}{2}\right),$
$Aa\sqrt{\frac{\lambda-E}{a}}\sin(\sqrt{\frac{\lambda-E}{a}}\,\frac{R}{2})=Bb\sqrt{\frac{\mu+E}{b}}\sinh\left(\sqrt{\frac{\mu+E}{b}}\,\frac{r}{2}\right),$

and thus

$\sqrt{a(\lambda-E)}\tan\left(\sqrt{\frac{\lambda-E}{a}}\,\frac{R}{2}\right)=\sqrt{b(\mu+E)}\tanh\left(\sqrt{\frac{\mu+E}{b}}\,\frac{r}{2}\right),$

which is the same as (46) but with $R/2$ and $r/2$ instead of $R$ ad $r$.
Hence the proof is completed as in Theorem 2.2.

Proof of Theorem 2.4.

Now a solution $y$ to the eigenvalue problem $Ly=Ey$ is obtained by gluing
solutions $y_{odd}$ and $y_{even}$ defined by formulas like (43) and (44) on
odd and even subintervals. The key point is, that under the periodic condition
the operator $L$ commutes with the shift operator $Tf(x)=f(x+R+r)$, where the
addition is understood modulo $R+r$. This means that for any eigenfunction $y$
of $L$, the function $Ty$ is also an eigenfunction with the same eigenvalue.
Using the nondegeneracy of eigenvalues of $L$ (like in the proof of the
previous theorem) we can conclude that $Ty=\omega y$ for some $\omega$. But
$T^{k}$ is the identity operator, so that $\omega^{k}=1$. Since our
eigenfunctions are real it follows that $\omega=\pm 1$. Again the
multiplication by $-1$ would contradict the continuity and thus only
$\omega=1$ is allowed. Hence we are directly in the setting of Theorem 2.3 and
the proof is complete.

Proof of Theorem 2.5 and 2.6 . Under the condition of the Theorem, there
exists an invertible matrix $C$ such that $CMC^{-1}=D$ is diagonal with the
elements $\Lambda_{j}$ on the diagonal. Then in the variables $z=Cy$ 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. And the equations for $z^{1}$ (and all boundary conditions)
become identical to the equations for $y$ of the one-dimensional case, with
$\Lambda_{1}$ instead of $\lambda$ and the new $\mu$ being equal to
$\lambda+\mu$ of the one-dimensional case.

Proof of Theorem 2.7 . The eigenvalue problem to equation (22) with the
Dirichlet boundary condition writes down as $Ay^{\prime\prime}+My=Ey$ subject
to $y(0)=y(R)=0$. Equivalently this equation rewrites as

$y^{\prime\prime}+(A^{-1}M-A^{-1}E)y=0.$ (48)

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 2.5, equation (48) is solvable for some positive $E$ if and only if
the highest real eigenvalue $\Lambda_{1}$ of the matrix $A^{-1}M$ satisfies
the condition $\Lambda_{1}>\pi^{2}/R^{2}$ implying the claim of the theorem.

Proof of Theorem 2.8 .

Change the variable in (24) to $z$ so that to $z_{ben}=A_{ben}y$,
$z_{nb}=A_{nb}y$. Then we get

$\displaystyle z^{\prime\prime}_{ben}+(M_{ben}-E)A_{ben}^{-1}z_{ben}=0$ (49)
$\displaystyle z^{\prime\prime}_{nb}+(M_{nb}-E)A_{nb}^{-1}z_{nb}=0,$

and the gluing condition for $z$ is the continuity of $z$ and $z^{\prime}$ on
the interfaces.

Let $L(E)$ denote the operator on the r.h.s. of (49), so that system (49) 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\geq 0$. This would be the
case, if one could show that $(z,L(E)z)<0$ for all $z$. But
$(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$. But $L^{S}(E)$ decreases with the
increase of $E$. Thus it is sufficient to show that $L=L(0)$ is negative, that
is $Lz=Pz$ has no solutions with non-negative $P$.

And this means that there is no solutions to the system

$\displaystyle z^{\prime\prime}_{ben}+(M^{S}_{ben}-P)z_{ben}=0$ (50)
$\displaystyle z^{\prime\prime}_{nb}+(M^{S}_{nb}-P)z_{nb}=0,$

with $P\geq 0$ and
$M^{S}_{ben}=(M_{ben}A_{ben}^{-1}+A_{ben}^{-1}M_{ben}^{T})/2$,
$M^{S}_{nb}=(M_{nb}A_{nb}^{-1}+A_{nb}^{-1}M_{nb}^{T})/2$.

It is clear that the required sufficient condition is sufficient to show for
the case when all eigenvalues of $M^{S}_{ben}$ and $M^{S}_{nb}$ are different,
as the general case is obtained by a straightforward limiting procedure.

For simplicity, let us assume that $M^{S}_{ben}$ has one positive eigenvalue
(that is $k=1$) and all other are negative:
$\lambda_{1}>0\geq\lambda_{2}>\cdots>\lambda_{n}$.

Clearly $P\geq\lambda_{1}$ are impossible in (50), so that only
$P\in[0,\lambda_{1}]$ have to be analyzed.

Solutions to (50) are given by the formulas

$z_{ben}=A_{1}\cos(\sqrt{\lambda_{1}-P}x)v_{1}+\sum_{j=2}^{n}A_{j}\cosh(\sqrt{P-\lambda_{j}}x)v_{j},$
$z_{nb}=\sum_{j=1}^{n}B_{j}\cosh(\sqrt{P-\mu_{j}}(x-(R+r)/2))w_{j},$

where $v_{j}$ and $w_{j}$ are orthonormal spectral bases for $M^{S}_{ben}$ and
$M^{S}_{nb}$.

Let $w_{j}=\sum_{k}c_{jk}v_{k}$ with an orthogonal matrix $C=(c_{ij})$. Then

$z_{nb}=\sum B_{k}\cosh(\sqrt{P-\mu_{k}}(x-(R+r)/2))c_{kj}v_{j}.$

Thus gluing $z_{ben}$ and $z_{nb}$ at $x=R/2$ yields

$A_{1}\cos(\sqrt{\lambda_{1}-P}R/2)=\sum_{j}B_{j}\cosh(\sqrt{P-\mu_{j}}r/2)c_{j1},$

and

$A_{i}\cosh(\sqrt{P-\lambda_{j}}R/2)=\sum_{j}B_{j}\cosh(\sqrt{P-\mu_{j}}r/2)c_{ji},\quad
i>1.$

Next, since

$\int_{-R/2}^{r+R/2}(z(x),z^{\prime\prime}(x))\,dx=-\int_{-R/2}^{r+R/2}(z^{\prime}(x),z^{\prime}(x))\,dx,$

which follows from the integration by parts and the gluing (and periodic)
conditions, we have that, for a solution to (50),

$0\leq\int_{ben}(z_{ben},(M^{S}_{ben}-P)z_{ben})\,dx+\int_{nb}(z_{nb},(M^{S}_{nb}-P)z_{nb})\,dx$
$=A_{1}^{2}(\lambda_{1}-P)\int_{ben}\cos^{2}(\sqrt{\lambda_{1}-P}x)\,dx+\sum_{j>1}A_{j}^{2}(\lambda_{j}-P)\int_{ben}\cosh^{2}(\sqrt{P-\lambda_{j}}x)\,dx$
$+\sum_{j\geq
1}B_{j}^{2}(\mu_{j}-P)\int_{nb}\cosh^{2}(\sqrt{P-\mu_{j}}(x-(R+r)/2))\,dx.$

Since $\cos(2x)=2\cos^{2}x-1$ and $\cosh(2x)=2\cosh^{2}x-1$,

$\int_{\xi}^{\xi}\cos^{2}(\alpha
x)\,dx=\xi+\frac{1}{2}\int_{\xi}^{\xi}\cos(2\alpha
x)\,dx=\xi+\frac{\sin(2\alpha\xi)}{2\alpha},$

and thus

$A_{1}^{2}(\lambda_{1}-P)\left(\frac{R}{2}+\frac{\sin(R\sqrt{\lambda_{1}-P})}{2\sqrt{\lambda_{1}-P}}\right)\geq\sum_{j\geq
1}B_{j}^{2}(P-\mu_{j})\left(\frac{r}{2}+\frac{\sinh(r\sqrt{P-\mu_{j}})}{2\sqrt{P-\mu_{j}}}\right).$
(51)

In particular,

$A_{1}^{2}\geq
B_{j}^{2}\frac{P-\mu_{j}}{\lambda_{1}-P}\frac{\left(\frac{r}{2}+\frac{\sinh(r\sqrt{P-\mu_{j}})}{2\sqrt{P-\mu_{j}}}\right)}{\left(\frac{R}{2}+\frac{\sin(R\sqrt{\lambda_{1}-P})}{2\sqrt{\lambda_{1}-P}}\right)}$

for any $j$.

One could work with this estimate, but it makes the final result more
transparent (though a bit more rough), if one simplifies this estimate to the
following one:

$A_{1}^{2}\geq
B_{j}^{2}\frac{\sqrt{P-\mu_{j}}\sinh(r\sqrt{P-\mu_{j}})}{2R(\lambda_{1}-P)}.$

On the other hand,

$|A_{1}|\cos(\sqrt{\lambda_{1}-P}\frac{R}{2})\leq\sum_{j}|B_{j}|\cosh(\sqrt{P-\mu_{j}}\frac{r}{2})c_{j1},$

and thus

$A_{1}^{2}\cos^{2}(\sqrt{\lambda_{1}-P}\frac{R}{2})\leq
n\sum_{j}B_{j}^{2}\cosh^{2}(\sqrt{P-\mu_{j}}\frac{r}{2})c^{2}_{j1}\leq
n\max_{j}B_{j}^{2}\cosh^{2}(\sqrt{P-\mu_{j}}\frac{r}{2}),$

or

$A_{1}^{2}(1+\cos(\sqrt{\lambda_{1}-P}R))\leq
n\max_{j}B_{j}^{2}(1+\cosh(\sqrt{P-\mu_{j}}r)).$

Thus for a $j$ realising the maximum we have

$n\frac{1+\cosh(\sqrt{P-\mu_{j}}r)}{1+\cos(\sqrt{\lambda_{1}-P}R)}\geq\frac{\sqrt{P-\mu_{j}}\sinh(r\sqrt{P-\mu_{j}})}{2R(\lambda_{1}-P)},$

or

$\frac{R(\lambda_{1}-P)}{(1+\cos(R\sqrt{\lambda_{1}-P}))}\geq\frac{\sqrt{P-\mu_{j}}\sinh(r\sqrt{P-\mu_{j}})}{2n(1+\cosh(r\sqrt{P-\mu_{j}}))}.$
(52)

This can be rewritten as

$\frac{1}{R}g_{1}(R\sqrt{\lambda_{1}-P})\geq\frac{1}{2nr}g_{2}(r\sqrt{P-\mu_{j}})$
(53)

with

$g_{1}(x)=\frac{x^{2}}{1+\cos x},\quad g_{2}(x)=\frac{x\sinh x}{1+\cosh x}.$
(54)

In order to make (53) impossible (and thus to ensure the negativity of our
operator) it is sufficient to assume that

$\frac{1}{R}g_{1}(R\sqrt{x})<\frac{1}{2nr}g_{2}(r\sqrt{\lambda_{1}-\mu_{1}-x})$
(55)

for $x\in[0,\lambda_{1}]$, where $x=\lambda_{1}-P$. For
$\lambda_{1}\geq\pi^{2}/R$, this cannot be true. Hence assume
$\lambda_{1}<\pi^{2}/R$.

Since it is straightforward to check that the function

$\frac{g_{1}(R\sqrt{x})}{g_{2}(r\sqrt{\lambda_{1}-\mu_{1}-x})}$

is increasing for $x\in[0,\lambda_{1}]$, condition (55) is equivalent to this
condition on the endpoint:

$\frac{1}{R}g_{1}(R\sqrt{\lambda_{1}})<\frac{1}{2nr}g_{2}(r\sqrt{|\mu_{1}|}),$
(56)

or explicitly

$\frac{2Rn\lambda_{1}}{1+\cos(R\sqrt{\lambda_{1}}))}<\frac{\sqrt{|\mu_{1}|}\sinh(r\sqrt{|\mu_{1}|})}{1+\cosh(r\sqrt{|\mu_{1}|})},$
(57)

which coincides with (26) for $k=1$.

Modifications required for arbitrary $k$ are straightforward, and we omit
them.

Proof of Theorem 2.9 .

The eigenvalue equations (29) can be written as

$\displaystyle y^{\prime\prime}_{ben}+N_{ben}(E)y_{ben}=0$ (58) $\displaystyle
y^{\prime\prime}_{nb}+N_{nb}(E)y_{nb}=0.$

As above we can write

$y_{ben}=A\cos(\sqrt{\Lambda_{1}}x)v_{1}+B\cosh(\sqrt{|\Lambda_{2}|}x)v_{2},$
$y_{nb}=C\cosh(\sqrt{|\mu_{1}|}(x-(R+r)/2))w_{1}+D\cosh(\sqrt{|\mu_{2}|}(x-(R+r)/2))w_{2}.$

The last equation rewrites as

$y_{nb}=C\cosh(\sqrt{|\mu_{1}|}(x-(R+r)/2))(c_{11}v_{1}+c_{12}v_{2})+D\cosh(\sqrt{|\mu_{2}|}(x-(R+r)/2))(c_{21}v_{1}+c_{22}v_{2}).$

Of course, all terms here depend on $E$.

Further on, in order to shorten the formulas, we shall write $\sin$, $\cos$,
$\tan$ if the argument is $\sqrt{\Lambda_{1}}R/2$ and $\sinh$, $\cosh$,
$\tanh$ if the argument is $\sqrt{|\Lambda_{2}|}R/2$. Gluing the solutions at
the interface $x=R$ yields

$A\cos=C\cosh(\sqrt{|\mu_{1}|}r/2)c_{11}+D\cosh(\sqrt{|\mu_{2}|}r/2)c_{21},$
$B\cosh=C\cosh(\sqrt{|\mu_{1}|}r/2)c_{12}+D\cosh(\sqrt{|\mu_{2}|}r/2)c_{22}.$

and

$-A\sqrt{\Lambda_{1}}\sin=-a(C\sqrt{|\mu_{1}|}\sinh(\sqrt{|\mu_{1}|}r/2)c_{11}+D\sqrt{|\mu_{2}|}\sinh(\sqrt{|\mu_{2}|}r/2)c_{21}),$
$B\sqrt{\Lambda_{1}}\sinh=-a(C\sqrt{|\mu_{1}|}\sinh(\sqrt{|\mu_{1}|}r/2)c_{12}+D\sqrt{|\mu_{2}|}\sinh(\sqrt{|\mu_{2}|}r/2)c_{22}).$

Equalising $A$ from the 1st equations yields

$\sqrt{\Lambda_{1}}\tan[C\cosh(\sqrt{|\mu_{1}|}r/2)c_{11}+D\cosh(\sqrt{|\mu_{2}|}r/2)c_{21}]$
$=a[C\sqrt{|\mu_{1}|}\sinh(\sqrt{|\mu_{1}|}r/2)c_{11}+D\sqrt{|\mu_{2}|}\sinh(\sqrt{|\mu_{2}|}r/2)c_{21}],$

Equalising $B$ from the 2nd equations yields

$\sqrt{|\Lambda_{2}|}\tanh[C\cosh(\sqrt{|\mu_{1}|}r/2)c_{12}+D\cosh(\sqrt{|\mu_{2}|}r/2)c_{22}]$
$=-a[C\sqrt{|\mu_{1}|}\sinh(\sqrt{|\mu_{1}|}r/2)c_{12}+D\sqrt{|\mu_{2}|}\sinh(\sqrt{|\mu_{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}[\sqrt{\Lambda_{1}}\tan-a\sqrt{|\mu_{1}|}\tanh(\sqrt{|\mu_{1}|}r/2)][\sqrt{|\Lambda_{2}|}\tanh+a\sqrt{|\mu_{2}|}\tanh(\sqrt{|\mu_{2}|}r/2)]$
$-c_{12}c_{21}[\sqrt{\Lambda_{1}}\tan-a\sqrt{|\mu_{2}|}\tanh(\sqrt{|\mu_{2}|}r/2)][\sqrt{|\Lambda_{2}|}\tanh+a\sqrt{|\mu_{1}|}\tanh(\sqrt{|\mu_{1}|}r/2)]=0,$

or as $\sqrt{\Lambda}_{1}\tan=aB/A$ with

$A=c_{11}c_{22}[\sqrt{|\Lambda_{2}|}\tanh+a\sqrt{|\mu_{2}|}\tanh(\sqrt{|\mu_{2}|}r/2)]$
$-c_{12}c_{21}[\sqrt{|\Lambda_{2}|}\tanh+a\sqrt{|\mu_{1}|}\tanh(\sqrt{|\mu_{1}|}r/2)],$
$B=c_{11}c_{22}\sqrt{|\mu_{1}|}\tanh(\sqrt{|\mu_{1}|}r/2)[\sqrt{|\Lambda_{2}|}\tanh+a\sqrt{|\mu_{2}|}\tanh(\sqrt{|\mu_{2}|}r/2)]$
$-c_{12}c_{21}\sqrt{|\mu_{2}|}\tanh(\sqrt{|\mu_{2}|}r/2)][\sqrt{|\Lambda_{2}|}\tanh+a\sqrt{|\mu_{1}|}\tanh(\sqrt{|\mu_{1}|}r/2)].$

Since $c_{12}c_{21}\leq 0$ and $c_{11}c_{22}\geq 0$ (and
$|\mu_{2}|\geq|\mu_{1}|)$, it follows that

$B\geq\det
c[a\sqrt{|\mu_{2}|}\tanh(\sqrt{|\mu_{2}|}r/2)+\sqrt{\Lambda_{2}}\tanh]\sqrt{|\mu_{1}|}\tanh(\sqrt{|\mu_{1}|}r/2),$
$0<A<\det
c[\sqrt{|\Lambda_{2}|}\tanh+a\sqrt{|\mu_{2}|}\tanh(\sqrt{|\mu_{2}|}r/2)],$

so that

$aB/A\geq a\sqrt{|\mu_{1}|}\tanh(\sqrt{|\mu_{1}|}r/2).$ (59)

Hence

$\sqrt{\Lambda}_{1}\tan(\sqrt{\lambda_{1}}R/2)\geq
a\sqrt{|\mu_{1}|}\tanh(\sqrt{|\mu_{1}|}r/2).$

Consequently, if (31) 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 2.10 .

As is easily seen, if inequalities (34) hold for $N_{ben}(0)$ and $N_{nb}(0)$,
then they hold for $N_{ben}(E)$ and $N_{nb}(E)$ with all positive $E$. Hence
it is sufficient to check condition (30) for $E=0$.

For a matrix

$A=\left(\begin{aligned} &-m_{1}\quad b_{1}\\\ &\,\,b_{2}\quad-
m_{2}\end{aligned}\right),$

the eigenvalues are

$\Lambda_{1,2}=\frac{1}{2}\left[-(m_{1}+m_{2})\pm\sqrt{(m_{1}-m_{2})^{2}+4b_{1}b_{2}}\right].$

Thus $\Lambda_{2}<0$ and $\Lambda_{1}>0$.

Looking for eigenvectors $v_{1,2}$, let us choose $v_{1}^{1}=1,v_{2}^{2}=1$.
Then

$v_{1}^{2}=\frac{\Lambda_{1}+m_{1}}{b_{1}}=\frac{b_{2}}{\Lambda_{1}+m_{2}},$
$v_{2}^{1}=\frac{b_{1}}{\Lambda_{2}+m_{1}}=\frac{\Lambda_{2}+m_{2}}{b_{2}}.$

Thus

$V=\left(\begin{aligned} &1\quad v_{2}^{1}\\\ &v_{1}^{2}\quad
1\end{aligned}\right),\quad W=\left(\begin{aligned} &1\quad w_{2}^{1}\\\
&w_{1}^{2}\quad 1\end{aligned}\right).$

Consequently,

$C=WV^{-1}=\frac{1}{\det V}\left(\begin{aligned} &1-v_{2}^{1}w_{1}^{2}\quad
w_{2}^{1}-v_{2}^{1}\\\ &-v_{1}^{2}+w_{1}^{2}\quad
1-v_{1}^{2}w_{2}^{1}\end{aligned}\right).$

Assuming $m_{1}>m_{2}$, we see that

$v_{2}^{1}=-\frac{(m_{1}-m_{2})+\sqrt{(m_{1}-m_{2})^{2}+4b_{1}b_{2}}}{2b_{2}}<0,$
$v_{1}^{2}=\frac{(m_{1}-m_{2})+\sqrt{(m_{1}-m_{2})^{2}+4b_{1}b_{2}}}{2b_{1}}<0.$

Similarly, $w_{2}^{1}<0$ and $w_{1}^{2}>0$. They have the same expression, but
with $\tilde{m}_{j}$ and $\tilde{b}_{j}$ instead of $m_{j},b_{j}$. Hence
$c_{11},c_{22}\geq 1>0$. Moreover,

$c_{12}c_{21}=-(w_{1}^{2}-v_{1}^{2})(|w_{2}^{1}|-|v_{2}^{1}|).$

Under (32) and (33), it follows that, for any $\omega\in(0,1)$,

$w_{1}^{2}\geq v_{1}^{2},|w_{2}^{1}|\geq|v_{2}^{1}|,$

and thus $c_{12}c_{21}<0$, which completes the proof.

##  5 Conclusion

In this paper we have developed the diffusion model for the propagation of
ticks with variable discontinuous diffusion coefficients. We have 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 have 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 have 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 most rough 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. Also 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 models aiming at the approximate description of the propagation in a
certain direction. Of course, it would be natural to extend the results to a
more realistic (for ticks and many other species) two-dimensional diffusion
models. 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 analysed based on the Fokas method, see [7] and [29], or via the
method of estimating the heat content, see e.g. [35]. It would be also of
interest to see whether the results of Theorem 2.9 have any reasonable
extensions beyond the restrictive assumptions (30) and whether one can get
some multi-dimensional analogs. Also 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. But 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. 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 stick to the diffusion models. Further investigations are
needed to see, whether more appropriate approximations can be achieved by
using other models, like the telegraph equations, or presently very popular
sub-diffusions or super-diffusions.

Notice 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 amount of studies were
devoted, see e.g. [19] and references therein.

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

## References

  * [1] J. W. Baron and T. Galla. Dispersal-induced instability in complex ecosystems. Nature Communications 11: 6032 (2020). https://doi.org/10.1038/s41467-020-19824-4
  * [2] T. Caraco, G. Gardener, W. Maniatty, E. Deelman and B. K. Szymanski. Lyme Disease: Self regulation and Pathogen Invasion. J. Theoretical Biology 193 (1998), 561-575.
  * [3] T. Caraco, S. Glavanakov, G. Chen, J. Flaherty, TK Ohsumi, BK Szymanski (2002). Stage-structured infection transmission and a spatial epidemic: A model for Lyme disease. American Naturalist (2002) 160: 348 - 359.
  * [4] J. P. Carvalho dos Santos, L. C. Cardoso, E. Monteiro and N. H. T. Lemes. A Fractional-Order Epidemic Model for Bovine Babesiosis Disease and Tick Populations. Abstract and Applied Analysis Volume 2015, Article ID 729894.
  * [5] Ch. L. Fefferman. The uncertainty principle. Bulletin of the Americal Math. Society 9:2 (1983), 129-206.
  * [6] N. A. Filippova (Ed.) Taiga tick Ixodes persulcatus Schulze (Acarina, Ixodidae): morphology, systematics, ecology, medical significance. Moscow, Nauka, 1985 (In Russian).
  * [7] A. S. Fokas. A Unified Approach to Boundary Value Problems. SIAM, Philadelphia, 2008.
  * [8] H. D. Gaff, L. J. Gross. Modeling Tick-Borne Disease: A Metapopulation Model. Bulletin of Mathematical Biology 69 (2007), 265–288.
  * [9] H. D. Gaff, E. Schaefer and S. Lenhart. Use of optimal control models to predict treatment time for managing tick-borne disease. Journal of Biological Dynamics, 5(5) (2011), 517-530.
  * [10] H. D. Gaff, A. White, K. Leas, P. Kelman, J, C. Squire, D. L. Livingston, G. A. Sullivan, E. W. Baker, D. E. Sonenshine. TickBot: A novel robotic device for controlling tick populations in the natural environment. Ticks and Tick-borne Diseases 6:2 (2015), 146 - 151.
  * [11] S. A. Gourley, N. F. Britton, M. A. J. Chaplain and H. M. Byrne. Mechanisms for stabilisation and destabilisation of systems of reaction-diffusion equations. J. Math. Biol. 34 (1996), 857 - 877.
  * [12] Harris, W. G. and E. C. Burns. Predation on the lone star tick by the imported fire ant. Environ. Entomol. 1 (1972), 362-365.
  * [13] E. E. Holmes. Are diffusion models too simple? A comparison with telegraph models of invasion. American Naturalist 142 (1993), 779 - 795.
  * [14] E. E. Holmes, M. A. Lewis, J. E. Banks and R. R. Veit. 1994\. Partial differential equations in ecology: spatial interactions and population dynamics. Ecology 75 (1994), 17 - 29.
  * [15] A. Kashkynbayev and D. Koptleuova. Global dynamics of tick-borne diseases Mathematical Biosciences Engineering (MBE), 17(4) (2020), 4064–4079.
  * [16] H. Kierstead and L. B. Slobodkin. The Size of Water Masses Containing Plankton Bloom. Journal of Marine Research, 12 (1953), 141 - 147.
  * [17] U. S. Korotkov. Life cycle of the taiga tick Ixodes Persulcatus in the conifer forests of the bottom land of the East Sayan mountain range. Parasotologya 48:1 (2014), 20-36 (in Russuan).
  * [18] L. D. Landau and E. M. Lifshitz. Quantum Mechanics. Pergamon Press. Second Edition 1965.
  * [19] Yu. Makarova, D. Han, S. Molchanov and E. Yarovaya. Branching random walks with immigration. Lyapunov stability. Markov Process. Related Fields 25:4 (2019), 683 - 708.
  * [20] M. Maliyoni, F. Chirove, H. D. Gaff, K. S. Govinder. A Stochastic Tick-Borne Disease Model: Exploring the Probability of Pathogen Persistence. Bull Math Biol 79 (2017), 1999 - 2021.
  * [21] M. Mimura and J. D. Murray. On a Diffusive Prey-Predator Model which Exhibits Patchiness. J. theor. Biol. 75 (1978), 249 - 262.
  * [22] G. A. Mount, D.G. Halle, D. R. Barnard and E. Daniels. New Version of LSTSIM for Computer Simulation of Amblyomma americanum (Acari: Ixodidae) Population Dynamics. J. Med. Entomol. 30: 5 (1993), 843 - 857.
  * [23] G. A. Mount, D.G. Halle and E. Daniels. Simulation of Blacklegged Tick (Acari: Ixodidae) Population Dynamics and Transmission of Borrelia burgdorferi. J. Med. Entomol. 34:4 (1997), 461 - 484.
  * [24] J. D. Murray. Mathematical Biology II: Spatial models and Biomedical Applications. Third Edition, Springer 2003.
  * [25] H. G. Mwambi, J. Baumgärtner and K. P. Hadeler. Development of a stage-structured analytical population model for strategic decision making: the case of ticks and tick-borne diseases. Riv. Mat. Univ. Parma (6) 3 (2000), 157 - 169.
  * [26] K. Nah and J. Wu. Long-term transmission dynamics of tick-borne diseases involving seasonal variation and co-feeding transmission. Journal of Biological Dynamics 15:1 (2021), 269-286.
  * [27] A. Okubo. Critical patch size for plankton and patchiness. In: S. A. Levin and T. G. Hallam (Eds.) Mathematical ecology. Lecture Notes in Biomathematics, 54, p. 456-477, Springer 1986.
  * [28] R. S. Ostfeld. The Ecology of Lyme-Disease Risk: Complex interactions between seemingly unconnected phenomena determine risks of exposure to this expanding disease. American Scientist 85:4 (1997), 338-346.
  * [29] N.E. Sheils and B. Deconinck. Initial-to-Interface Maps for the Heat Equation on Composite Domains. Studies in Appied Mathematics 137:1 (2016), 140-154.
  * [30] H. Shu, W. Xu, X.-Sh. Wang and J. Wu. Complex dynamics in a delay differential equation with two delays in tick growth with diapause. J. Differential Equations 269(2020), 10937 - 10963.
  * [31] J. G. Skellam. Random Dispersal in theoretical biology. Bulletin of Mathematical Biology 53: 1/2 (1991), 135 - 165. Reprinted from Biometrika 38 (1951), 196-218.
  * [32] C. Sweezy. Relating different conditions for the positivity of the Schrödinger operator. Rocky Mountain J. of Maths. 23:1 (1993), 353-366.
  * [33] J. Switkes, B. Nannyonga, J.Y.T. Mugisha, J. Nakakawa. A mathematical model for Crimean-Congo haemorrhagic fever: tick-borne dynamics with conferred host immunity, Journal of Biological Dynamics (2016) 10:1, 59 - 70.
  * [34] M. Tosato, K. Nah and J. Wu. Are host control strategies effective to eradicate tick-borne diseases (TBD)? Journal of Theoretical Biology 508 (2021), 110483.
  * [35] M. van den Berg and K. Gittins. Uniform bounds for the heat content of open sets in Euclidean space. Differential Geom. Appl., 40 (2015), 67 - 85.
  * [36] O. A. Vshivkova. Mathematical simulation of ant influence on the ixode ticks number in Euroasian ecosystems. Izvestia Samarskogo nauchnogo centra Russian Academy of Science 11:197 (2009), 1631-1633 (In Russian).
  * [37] O. A. Vshivkova, A. S. Komarov, P. V. Frolov, R.G. Khlebopros. The role of the heterogeneity of the habitat under the control of the ixode ticks number: cellular-automaton model. Control Scienes 4 (2013), 57-63 (In Russian).
  * [38] A. White, E. Schaefer, Ch. Wright Thompson, Ch. M. Kribs and H. Gaff. Dynamics of two pathogens in a single tick population. Letters in Biomathematics 6:1 (2019), 50–66.
  * [39] Xue Zhang, Bei Sun, Yijun Lou. Dynamics of a periodic tick-borne disease model with co-feeding and multiple patches. Journal of Mathematical Biology (2021) 82:27.
  * [40] Xiao-Qiang Zhao. Global dynamics of a reaction and diffusion model for Lyme disease. J. Math. Biol. (2012), 65: 787 - 808.