On the Bounds for a Two-Dimensional Birth-Death Process with Catastrophes

Anna Sinitcina 1,†, Yacov Satin 1,†, Alexander Zeifman 2,*,† ID , Galina Shilova 1,†, Alexander Sipin 1,†, Ksenia Kiseleva 1,†, Tatyana Panfilova 1,†, Anastasia Kryukova 1,†, Irina Gudkova 3,† and Elena Fokicheva 1,† 1 Faculty of Applied Mathematics, Computer Technologies and Physics, Vologda State University, 160000 Vologda, Russia; a_korotysheva@mail.ru (An.S.); yacovi@mail.ru (Y.S.); shgn@mail.ru (G.S.); cac1909@mail.ru (Al.S.); ksushakiseleva@mail.ru (K.K.); ptl-70@mail.ru (T.P.); krukovanastya25@mail.ru (A.K.); eafokicheva2007@yandex.ru (E.F.) 2 Faculty of Applied Mathematics, Computer Technologies and Physics, Vologda State University, Institute of Informatics Problems, Federal Research Center “Computer Science and Control” of the Russian Academy of Sciences, Vologda Research Center of the Russian Academy of SciencesSciences, 160000 Vologda, Russia 3 Applied Probability and Informatics Department, Peoples’ Friendship University of Russia (RUDN University), Institute of Informatics Problems, Federal Research Center “Computer Science and Control” of the Russian Academy of Sciences, 117198 Moskva, Russia; gudkova_ia@rudn.university * Correspondence: a_zeifman@mail.ru † These authors contributed equally to this work.


Introduction
There is a large number of papers devoted to the research of continuous-time Markov chains and models with possible catastrophes, see for instance , and the references therein.Such models are widely used in queueing theory and biology, particularly, for simulations in hight-performance computing.In some recent papers, the authors deal with more or less special birth-death processes with additional transitions from and to origin [9][10][11][12][13][18][19][20].In [22], a general class of Markovian queueing models with possible catastrophes is analyzed and some bounds on the rate of convergence are obtained.Here we consider a more specific but important model of a two-dimensional birth-death process with possible catastrophes and obtain the upper bounds on the rate of convergence in some weighted norms and the corresponding perturbation bounds.
Ergodicity bounds in l-1 norm (associated with total variation) for such processes can be obtained quite easily due to the possibility of catastrophes, i.e., transitions to zero from any other state.Obtaining the estimates in weighted norms that guarantee the convergence of the corresponding mathematical expectations as well as the construction of the corresponding limiting characteristics are more complex problems.
In addition, we consider in detail two examples with 1-periodic intensities and various types of death (service) rates.The bounds on the rate of convergence and the behavior of the corresponding mathematical expectations are obtained for each example.
Our results seem to be interesting for both queueing theory and biology applications.Let X(t) = (X 1 (t), X 2 (t)) be two-dimensional birth-death-catastrophe process (where X i (t) is the corresponding number of particles of type i, i = 1, 2) such that in the interval (t, t + h) the following transitions are possible with order h: birth of a particle of type i, death of a particle of type i, and catastrophe (or transition to the zero state 0 = (0, 0)).
Denote by λ 1,i,j (t), λ 2,i,j (t), µ 1,i,j (t), µ 2,i,j (t), and by ξ i,j (t) corresponding birth, death, and catastrophe rates for the process.Namely, λ 1,i,j (t) is the rate of transition from state (i, j) to state (i + 1, j) at the moment t, λ 2,i,j (t) is the rate of transition from state (i, j) to state (i, j + 1), µ 1,i,j (t) is the rate of transition from state (i, j) to state (i − 1, j), µ 2,i,j (t) is the rate of transition from state (i, j) to state (i, j − 1), and finally, ξ i,j (t) is the rate of transition from state (i, j) to state (0, 0) at the moment t.
The transition rate diagram associated with the process is presented in Figure 1.Suppose that all intensities are nonnegative and locally integrable on [0; ∞) as functions of t.Moreover, we also suppose that the condition of boundedness hold for any i, j and almost all t ≥ 0. We renumber the states of two-dimensional process X(t) = (X 1 (t), X 2 (t)) (0,0), (0,1), (1,0), (0,2), (1,1), (2,0), . . .by increasing the sum of coordinates, and in the case of the same sum, by increasing the first coordinate.Hence we obtain one-dimensional vector p(t) = (p 0 (t), p 1 (t), . . . ) T of state probabilities in a new numeration, and therefore, we can rewrite the forward Kolmogorov system in the following form: where A(t) = a ij (t) is the corresponding transposed intensity matrix: , and a ii (t) = − ∑ i a i,j (t).
Throughout the paper by • we denote the l 1 -norm, i. e., x = ∑ |x i |, and Let Ω be a set all stochastic vectors, i.e., l 1 vectors are with nonnegative coordinates and unit norm.Hence the assumption (1) implies the bound A(t) ≤ 2L for almost all t ≥ 0. Therefore, the operator function A(t) from l 1 into itself is bounded for almost all t ≥ 0 and locally integrable on [0; ∞).Therefore, we can consider the forward Kolmogorov system as a differential equation in the space l 1 with bounded operator.
It is well known, see [23], that the Cauchy problem for such a differential equation has a unique solution for an arbitrary initial condition, and p(s) ∈ Ω implies p(t) ∈ Ω for t ≥ s ≥ 0.
We have where U(t, s) is the Cauchy operator of Equation (2).

Bounds in l 1 Norm
Consider the first equation in forward Kolmogorov system and rewrite it in the following form: where ξ(t) = inf i,j ξ i,j (t).
Then one can estimate the logarithmic norm of B(t) in the space of sequences l 1 (see [24]): Then for all 0 ≤ s ≤ t we have Therefore, the following statement is correct (see details in [16,18]).
Theorem 1.Let the intensities of catastrophes be essential, that is Then the process X(t) is weakly ergodic and the following bound of the rate of convergence holds: for all initial conditions p * (0), p * * (0) and any t ≥ 0.
Consider now the "perturbed" process X = X(t), t ≥ 0, adding a dash on top for all corresponding characteristics.
Put Â(t) = A(t) − Ā(t), and assume that the perturbations are "uniformly small", i.e., for almost all t ≥ 0 the following inequality is correct Consider the stability bounds of the process X(t) under these perturbations.In addition, we assume that the process is exponentially ergodic, that is, that for some positive M, a and for all s, t, 0 ≤ s ≤ t the following inequality holds Then from Theorem 1: Here we apply the approach proposed in [25] for a stationary case and generalized for a nonstationary situation in [15,16].
We have Therefore we obtain It implies the following statement.
Theorem 2. If the condition ( 12) is fulfilled and the perturbations are uniformly small: for almost all t ≥ 0. Then the following bound holds: for any initial conditions p(0), p(0).

Corollary 1.
Let the intensities of the process be 1-periodic and instead of ( 12) we have the following inequality Then ( 17) is correct for where K = sup |t−s|≤1 t s ξ(τ)dτ < ∞.

Bounds in Weighted Norms
Consider the diagonal matrix D = diag(d 0 , d 1 , d 2 , d 3 , • • • ), with entries of the increasing sequence {d n }, where d 0 = 1, and the corresponding space of sequences l 1D = z = (p 0 , p 1 , p 2 , . ..)T such that Then one can estimate the logarithmic norm of operator B(t) in l 1D space.
According to the general approach, we obtain the matrix , where a ii (t) = − ∑ i a i,j (t).
Consider now the logarithmic norm Let us make the correspondence between the column number of matrix DB(t)D −1 and the number of zeros under the main diagonal in this column (till the first nonzero element).Then we obtain the arithmetic progression {a i }: We compose the sequence {b i } of the number of identical entries of the third line: 2, 3, 4, 5, 6, • • • .Note that ∑ N i=1 b i is equal to the last a k , corresponding to the number of zeros N in the k-th column.Then the sum of the first N elements of sequence {b i } is approximately equal to the number of the column a n = n: Knowing the column number n, one can find the formula for the number of zeros N under the main diagonal in this column till the first nonzero element.We note that the number of zeros over the diagonal till µ 1,i,j is one less.
If N is not an integer, we must take the nearest right to N an integer.One can see that columns 2, 5, 9, 14, ... (these columns correspond to sums ∑ j i=1 b i and integer N) contain death rates µ 2,i,j (t) only, and columns 3, 6, 10, 15, ... contain death rates µ 1,i,j (t) only, and all other columns contain the both death intensities.
Consider the following quantities: , in other cases: Then the following algorithm helps us to correlate the number n and pair (i, j): Therefore, for all 0 ≤ s ≤ t we have the bound for the corresponding Cauchy operator: and the following statement.
Theorem 3. Let for some sequence {d i } we have the condition Then the process X(t) is weakly ergodic and the following bound of the rate of convergence is correct: for any initial conditions p * (0), p * * (0) and for all t ≥ 0.
Mathematical expectations for both processes X 1 (t) and X 2 (t) can be obtained using formulas: , that is the number of all particles at the moment t.
Then one has for the mathematical expectation (the mean) of this process the following equality: We note that for W = inf i≥0 d i i the next inequality holds Denote by the conditional expected number of all particles in the system at instant t, provided that initially (at instant t = 0) k particles of both types were present in the system.Corollary 2. Let the condition (23) hold and there is a sequence {d i } such that W > 0, then the process N(t) has the limiting mean φ N (t) = E N (t, 0), and for any j and all t ≥ 0 the following bound of the rate of convergence is correct: Applying in addition the condition that the perturbations of the intensity matrix are small enough in the corresponding norm, that is Â(t) 1D ≤ ε, one can also obtain B(t) 1D ≤ ε.
We assume here that the process X(t) is exponentially ergodic in l 1D -norm, that is for some positive M 1 , a 1 and for all s, t, 0 ≤ s ≤ t the following inequality holds: Here we apply the approach from [18].One can rewrite the original system for the unperturbed process in the form: Then and Therefore, in any norm for any initial conditions we have the correct bound: Then we have the following inequality for the logarithmic norm: On the other hand, one can obtain the estimation using inequality (29): for any initial condition p(0).Moreover, f(τ) 1D ≤ ε.
Then using bound (33), we have Therefore, the following statement is correct.
Theorem 4. Let inequalities ( 12) and (29) hold for any initial condition p(0) ∈ l 1D and for all t ≥ 0, then we have and Corollary 3. Let in addition the sequence be increasing fast enough, such that W > 0, then for any j, t ≥ 0 we have

Example 2. Let now
The mean for the process N(t) on the interval t ∈ [0, 3] for different initial conditions are shown in Figures 9-12, and the bounds for the limiting perturbed mean is shown in Figure 13.Remark 1.These graphs give us the additional information on the considered examples.Namely, Figures 2 and 8 show the bounding on the rate of convergence, see Equation ( 21) for Examples 1, 2 respectively, in Figures 3-6 and 9-12 one can see the mathematical expectation of the number of all particles at the moment t until the stationary behaviour.Finally, the limiting behaviour of the limiting mathematical expectation of the number of all particles for the perturbed process is shown in Figures 7 and 13.

Figure 2 .
Figure 2. The values of several α n (t) for Example 1.

Figure 8 .
Figure 8.The values of several α n (t) for Example 2.