key: cord-0254158-s6hfr392
authors: Makarova, Iu.; Balashova, D.; Molchanov, S.; Yarovaya, E.
title: Branching Random Walks with Two Types of Particles on Multidimensional Lattices
date: 2022-02-06
journal: nan
DOI: 10.3390/math10060867 10.3390/math10060867
sha: f58117a59e2bc5e9f6b48476b98493b849f4cf90
doc_id: 254158
cord_uid: s6hfr392

We consider a continuous-time branching random walk on a multidimensional lattice with two types of particles and an infinite number of initial particles. The main results are devoted to the study of the generating function and the limiting behavior of the moments of subpopulations generated by a single particle of each type. We assume that particle types differ from each other not only by the laws of branching, as in multi-type branching processes, but also by the laws of walking. For a critical branching process at each lattice point and recurrent random walk of particles, the effect of limit spatial clustering of particles over the lattice is studied. A model illustrating epidemic propagation is also considered. In this model, we consider two types of particles: infected and immunity generated. Initially, there is an infected particle that can infect others. Here, for the local number of particles of each type at a lattice point, we study the moments and their limiting behavior. Additionally, the effect of intermittency of the infected particles is studied for a supercritical branching process at each lattice point. Simulations are presented to demonstrate the effect of limit clustering for the epidemiological model.

The branching random walk (BRW) is one of the widely used tools for describing the processes associated with the birth, death, migration, and immigration of particles [19, 6, 2, 3] . BRWs occur in population dynamics [11] and have numerous applications, e.g., in genetics [10] and demography [14] .

The continuous-time BRWs presented in this paper are the two-type branching processes with walking of particles which takes place on the multidimensional lattice Z d , d ∈ N. We will mainly study the distribution of subpopulations generated by a single particle of each type. We assume that each particle can produce not only particles of the same type, but also particles of a different type. We also assume that particles cannot change their type. However, later in Section 6, we lift this condition for a particular case.

The problem of multi-type processes is one of the most interesting and really complicated problems in the theory of random processes. Historically, these processes were apparently first considered by Sevastyanov in [15] . He considered both discrete and continuous-time branching processes with a finite number of types and studied the limit distribution of particles under different conditions. Nowadays this problem is also studied in detail by various research groups. In [4, 5] , for example, the authors consider a more complicated problem in which the number of types is not finite but countable. They consider the subclass of Galton-Watson processes called lower Hessenberg branching processes. In contrast to our studies, they investigated processes with discrete time, and the random walk was considered in a strip. Some works of Vatutin with coauthors [18, 17] are devoted to the problem of multi-type branching processes with discrete time and finite number of types in the random environment, but without walking of particles.

The structure of the paper is as follows. In Section 2, we describe a two-type BRW on Z d with infinitely many initial particles of both types. Here, we also define the main objects of the study. In Section 3, we study the first moments for subpopulations generated by a single particle of each type, and find their asymptotic behavior at the sites of Z d . To this end, we first obtain differential equations for the generating functions of subpopulations generated by a single particle of each type in Lemma 3.1 from Section 3.1. In order to find solutions for the corresponding equations, we turn to the equations for the Fourier transforms of the corresponding moments in Section 3.2 and show that the corresponding equations can be solved explicitly. This allows us to obtain explicit solutions (41) and (43) for the Fourier transforms of the corresponding moments. The results are then applied in Section 3.3 to find the asymptotics of the solutions for the Fourier transforms of the first moments of the subpopulations in the case of finite variance of the jumps. In Section 4, we study the second moments for the subpopulations. In Section 5, we study the particle clustering effect for BRWs under additional assumptions. Here, we assume that a two-type branching process is critical at every point on Z d and the distribution of an underlying random walk jumps has light tails. In Section 6, we discard one of the assumptions we imposed on our model in and instead assume that the particles of the first type can occasionally change their type to the second. This model can illustrate the situation related to epidemic spread, especially the spread of COVID-19 around the world. We refer to the first type of particles as infected and the second type of particles as immunity generated against COVID-19. Here, we assume that infected particles change their type after a short period of time so that they build up immunity. However, we assume that only one particle can do this after a short period of time. In Section 7, the algorithm for modeling the processes studied in Sections 5 and 6 is presented and examined using the Python programming language.

Here, we consider a population model with two types of particles. Let N i (t, y), with i = 1, 2, be the number of particles of type i at time t > 0 at the site y ∈ Z d , d 1. Then, the total population at the point y ∈ Z d at time t > 0 can be represented as the following column-vector whose components are non-negative integers:

N (t, y) = [N 1 (t, y), N 2 (t, y)] T .

(1)

We assume that N i (0, x) = l i for i = 1, 2 and all x ∈ Z d . We assume that the evolution of particles of each type consists of several possibilities. First, a particle of type i, i = 1, 2, can die with mortality rate µ i 0. Second, each particle of type i can produce new particles of either type. We denote by β i (k, l) 0, k + l 2, the rate at which a particle of type i produces k particles of type i = 1 and l particles of type i = 2. Then, we define the corresponding branching generation function (without particle death) for i = 1, 2, see, e.g., [15] :

Remark 2.1. In our notation, µ i = β i (0, 0) (i = 1, 2), and we do not consider the case when a particle of the type i = 1, 2 can transform to a particle of type j = 1, 2, j = i, and hence β 1 (0, 1) = β 2 (1, 0) = 0. Additionally, we assume that

where β 1 (1, 0) and β 2 (0, 1) denote the cases when nothing happens with the particles.

Remember that particles can jump between points on the lattice. We assume that the probability of a jump from a point x to a point x+v during the small period dt is equal to κ i a i (x, x+v) dt+ o(dt), i = 1, 2. Here, κ i > 0 is the diffusion coefficient. In what follows, we consider a symmetric random walk, i.e., the case where a i (x, y) = a i (y, x). Moreover, we assume that the random walk is homogeneous in space: a i (x, x + v) = a i (v) and irreducible such that span{v : a i (v) > 0} = Z d . Moreover, a i (0) = −1, v a i (v) = 0.

Then, the migration operator has the form

Let us introduce the subpopulations, which can be represented as the following column-vectors:

x, y) = [n 11 (t, x, y), n 12 (t, x, y)] T , n 2 (t, x, y) = [n 21 (t, x, y), n 22 (t, x, y)] T .

Here, n i (t, x, y) is the vector of particles at the point y, generated by a single particle of type i which at time moment t = 0 was at the site x ∈ Z d . Its components n ij (t, x, y) are the numbers of particles at the point y of type j, generated by a single particle of type i at x at the moment t = 0. Note that n ij (0, x, y) = δ i (j)δ x (y),

where δ u (v) is the Kronecker function on Z d (or R), that is if u, v ∈ Z d (or R)

Remark 2.2. In the BRW under consideration we assume that both random walk and branching process are "homogeneous". Namely, we assume that underlying random walk for each type of particles i = 1, 2 is homogeneous in space, so that a i (x, y) = a i (x − y, 0) = a i (x − y). At the same time branching process (which includes death and birth of particles) is also "homogeneous" due to the fact that all intensities µ i , β i (k, l), k + l 2, i = 1, 2 are constant and only depend on the type of particles (and independent of the lattice points). Such a "homogeneity" leads to the simplifying the relations which describe the evolution of considered BRW. First of all, we conclude that for all t 0 the probability P n ij (t, x, y) = k equals to P n ij (t, x − y, 0) = k for all k ∈ Z + , so that P n ij (t, x, y) = k ≡ P n ij (t, x − y, 0) = k , t 0.

To prove this equality, we consider the process n ij (t, x, y), which starts at some lattice point x, so that n ij (0, x, y) = δ i (j)δ x (y). Then, for each trajectory

which describes the transition of a particle from point x to the point y, there exists the "trajectory with a shift of y"

which describes the transition of a particle from point x − y to the point 0. At the same time, due to the homogeneity in space of random walk, all transition intensities for both trajectories are equal (a i (x, y) = a i (x − y, 0) = a i (x − y)), and because of the "branching homogeneity", all branching intensities are equal at every lattice point. As n ij (0, x − y, 0) = n ij (0, x, y), then for all t 0 we can conclude that Equation (6) is true. From (6), we get that for all values which can be obtained from n ij (t, x, y) we have the same relations. In particular, for En ij (t, x, y) we get

Finally, note that the similar relation [13] holds for transition probabilities p i (t, x, y), i = 1, 2 (definition will be given later)

The proof of the previous relation can be also obtained from the representation (47) from Section 3.3.

Remark 2.3. From Equations (6) and (8) obtained in Remark 2.2 we conclude that to investigate considered BRW, which starts at the lattice point x it is sufficient to consider the case x = 0. That can simplify the future narration. Now, using notation from Equation (4), we obtain the following representation of the total population specified by Equation (1):

where n i,l (t, x, y) is the subpopulation generated by the l-th particle at the point x at the time t = 0. Note that both internal series in Equation (9) do not depend on the order of enumeration of particles. The components of the vector N (t, y) are

where i = 1, 2 Given z = (z 1 , z 2 ), let us introduce the generating function

This generating function specifies the evolution of a single particle of type i = 1, 2. Let us consider what can happen to this particle (later we can use it to obtain a differential equation for the generating functions). First, the initial particle can die at a point x with probability µ i dt+o(dt) (then the subpopulation of this particle disappears). Second, this particle can produce k particles of type 1 and l particles of type 2 with probability β i (k, l) dt + o(dt). Third, the particle can jump from a point x to a point x + v with probability κ i a i (v) dt + o(dt). Finally, nothing can happen to a particle during time dt. From this, we get Lemma 2.1. The generating functions Φ i (t, x, y; z), i = 1, 2, specified by Equation (11), satisfy the differential equation

x, y; z))

Proof. Given an i = 1, 2, consider the generating function Φ i (t, x, y; z) at the time moment t + dt:

x, y; z)

Then,

Therefore,

x, y; z)).

Here, according to Equation (2), we have

and hence

x, y; z))

x, y; z)).

The initial condition for the latter equation follows from Equation (5):

So, we obtain the desired Equations (12) and (13) , which completes the proof of Lemma 2.1.

Remark 2.4. If we assume

for some c 0 > 0, then the Carleman condition is hold [21] , which guarantees that for each i = 1, 2 the function F i (z 1 , z 2 ) from Equation (2) is an analytic function in the strip |z i − 1| < δ 0 for some δ 0 > 0 [16] .

Recall that the goal of our article is to study the moments of the random variables n ij (t, x, y), i, j = 1, 2. In this section, we will consider the first moments. For this purpose, in Section 3.1, in Lemma 3.1, we obtain differential equations for the generating functions of subpopulations generated by a single particle of each type. In order to find solutions to the corresponding equations, we turn to the equations for the Fourier transforms of the corresponding moments in Section 3.2 and show that the corresponding equations can be solved explicitly. This allows us to obtain explicit solutions (41) and (43) for the Fourier transforms of the corresponding moments. The results are then applied in Section 3.3 to find the asymptotics of the solutions for the Fourier transforms of the first moments of subpopulations in the case of finite variance of the jumps.

Define m (1) ij (t, x, y) = En ij (t, x, y) and prove the following lemma playing important role in what follows. ij (t, x, y) satisfy the differential equation

Proof. Differentiating Equation (11) with respect to z j , j = 1, 2, we get

from which, by taking z = (z 1 , z 2 ) = (1, 1), we obtain

Now, differentiating Equation (12) over z j we can write

x, y; z) .

Again, by taking z = (z 1 , z 2 ) = (1, 1) in the above formula and applying Equation (17), we find that the left-hand part of the last equation takes the form:

while the right-hand part of the same equation equals

1j (t, x, y) + lm (1) 2j (t, x, y).

By combining Equation (18) and (19) , we obtain

1j (t, x, y) + lm (1) 2j (t, x, y)).

The initial condition for the latter equation can be found from Equation (5):

Lemma 3.1 is proved.

Remark 3.1. From Lemma 3.1 and the general theory of differential equations (in Banach spaces), for any i, j = 1, 2 one can easily obtain the inequality

(see, e.g., the proof of similar facts in [19, 20] ), but in order not to overload the exposition, its elementary proof is given in Section 3.2.

Nevertheless, let us explain the main ideas of the corresponding proof. Equation (15) with initial condition (16) can be treated as a linear differential equation in a Banach space whose right-hand side (for each t and y) is a linear bounded operator acting in any of the spaces l p (Z d ), p ≥ 1. Since in this case the initial condition m (1) ij (0, x, y) for each y as a function of the variable x also belongs to each of the spaces l p (Z d ), p ≥ 1, then, as shown, for example, in [19, 20] , m (1) ij (t, x, y) (for each t and y) as a function of the variable x also belongs to each of the spaces l p (Z d ), p ≥ 1, and is thus bounded.

In Lemma 3.1, we have obtained the differential equations for the subpopulations generated by a single particle of each type. Now, we want to obtain the differential equation for the full population N (t, y).

Define m (1) i (t, x, y) = En i (t, x, y), i = 1, 2, and rewrite the Equations (15) from Lemma 3.1 in the following form:

2 (t, x, y)

with the initial conditions m (1)

Let us denote n(t, x, y) = [n 1 (t, x, y), n 2 (t, x, y)] T and m (1) (t, x, y) = En(t, x, y). Then, the pair of equations, Equations (20) and (21) , can be rewritten in a more compact form:

2 (t, ·, y))(x)

where V is the matrix

The above calculus let us the opportunity to get the equation for the full population at the site y ∈ Z d . Using the representation of N i (t, y), i = 1, 2, in Equation (10) we obtain for m (1) 

2 (t, y)] T the following formula:

2i,m (t, x, y).

Taking the partial derivative over the parameter t for each component of m (1) (t, y) we derive from the above formula the equation:

where m (1) i,l (t, x, y) := En i,l (t, x, y). Formula (22) describes how the behavior of the full population depends on the behavior of each subpopulation. Later on, we will study the behavior of subpopulations in more details.

In this section, we will find an explicit form of the solutions of the differential equations obtained in Lemma 3.1. To find these solutions, we will use the discrete Fourier transform. For simplicity, we recall that the Fourier transform f (θ) of a function f (u) is defined as

where (·, ·) is the dot product in R d × R d , while the inverse Fourier transform is of the form

By applying the Fourier transform (23) to Equations (20) and (21), we obtain the equations ∂ m

2 (t, θ, y)

∂ m 

with the initial conditions m

To simplify formulas (25) and (26), let us introduce the following notations:

With the usage of these notations, Equations (25) and (26) can be represented in a more compact form:

To get a solution for this last system of differential equations, let us recall some facts from the theory of two-dimensional linear differential equations and perform some auxiliary calculations.

Remark 3.2. Represent Equations (31) and (32) arising in our treatment in a conventional form of a system of linear differential equations with two variables (see details in [7] ):

assuming that a, b, c and d here are some numerical parameters. In order to "keep the connection" with Equations (31) and (32) and not consider options unnecessary in the future, we will assume throughout this remark that b, c 0.

As is known (see, e.g., [7] or some other handbook on the theory of differential equations), the behavior of solutions of Equations (33) and (34) is completely determined, in a sense, by the roots of characteristic equation of the matrix of coefficients standing in the right-hand side of Equations (33) and (34):

These roots are as follows

Note that under the assumption b, c 0, the discriminant D is non-negative, and therefore the roots In this case, λ 1 = λ 2 coincide with each other and moreover λ 1 = λ 2 = a = d. Then (see, e.g., [7] ), the solution u(t) of Equations (33) and (34) is a linear combination of the functions e λt and te λt :

The solution v(t) can be expressed likewise. Let D = 0; this can be if and only if a = d or b = 0 and c = 0.

In this case, λ 1 = λ 2 and (see, e.g., [7] ) the solution u(t) of Equations (33) and (34) is a linear combination of the functions e λ1t and e λ2t :

The solution v(t) can be expressed likewise. Let us write out the precise forms of the solutions u(t) and v(t) of Equations (33) and (34); they will be needed in the further analysis. Consider the following combinations of the parameters b and c:

which exhaust all possible combinations of these parameters under condition b, c 0. The fact that the first two of these conditions intersect does not interfere with further considerations. We also mention that the case b = c = 0 is covered by both of the first two cases.

Case b = 0, c 0. Here, u(t) can be found directly from Equation (33):

To find v(t) it suffices to substitute the obtained expression for u(t) into Equation (34) and to solve the resulted non-homogeneous linear differential equation:

The value of the integral in the right-hand side of the obtained equality is different depending on whether the equality a = d holds. Direct evaluation shows that

Case b 0, c = 0. This case is treated similarly to the previous one, and we get:

Case b > 0, c > 0. In this case, both roots λ 1,2 of Equation (35) are different, and moreover λ 1 > λ 2 . In order to find the solutions u(t) and v(t) of Equations (33) and (34) let us first take t = 0 in Equation (36). Then, we obtain the following equation for the initial condition u(0):

Further, find bv(t) from Equation (33):

Using the obtained expression we will get the equation for bv(0):

By solving the resulting system of Equations (37) and (38), we get

Here, the last equation can be simplified by noting that

As a result, we obtain:

Now, we are able to write out the solutions of Equations (31) and (32). For this, it suffices to note that although in the reasoning of Remark 3.2 it was implicitly assumed that the functions u(t) and v(t) are scalar, but in fact this assumption was never used anywhere, and the functions u(t) and v(t) may be assumed vector-valued, for example, such as m (1) i (t, θ, y) in Equations (31) and (32).

One should also pay attention to the fact that in Equations (31) and (32), in contrast to Equations (33) and (34), the parameters a and d are actually functions of the variable θ, that is, a = a(θ) and d = d(θ), and then the values λ 1 , λ 2 and D are also functions of the variable θ:

and

Considering the above, we can write out the solutions m 

2 (t, θ, y) of Equations (31) and (32) using the appropriate initial conditions.

2 (0, θ, y).

Finally, it needs to be remembered that each of the functions m 2 (t, θ, y) for the obtained above three cases, we will obtain three cases of formulas for m 

Remark 3.3. The cases b = 0, c > 0 or b > 0, c = 0 can describe the situation when particles of one type cannot produce offsprings of both types. This can have the real-life interpretation: we have some species which have both "male" and "female" individuals and the "male" individuals cannot produce offspring.

Remark 3.4. The attentive reader will notice that our constructions are redundant in a sense. In the middle of this section, we made an effort to go from equations for the functions m (1) ij (t, x, y), i, j = 1, 2, to more general equations for the functions m ij (t, x, y), i, j = 1, 2 (or rather to their Fourier images m (1) ij (t, θ, y), i, j = 1, 2). We emphasize once again that from a technical point of view, this method of research is redundant, however, in our opinion, it contributes to a deeper understanding of the "nature of things" when analyzing the behavior of the functions m ij (t, x, y), i, j = 1, 2. Later, we will find the asymptotic behavior of each subpopulation m (1) ij (t, x, y) in one particular case.

In the previous section, in Equations (41)-(43), we found the solutions for the Fourier transform of the first moments of the subpopulations m ij (t, θ, y), i, j = 1, 2. In this section, we obtain their asymptotic behavior in one particular case that is natural for applications.

Remark 3.5. Consider the parabolic problem

where operators L i , i = 1, 2, are defined in (3). By applying the discrete Fourier transform (23) to Equation (44), we find that the Fourier image p(t, θ, y) of the function p(t, x, y) satisfies the Cauchy problem

whose solution can be found explicitly:

Applying the inverse Fourier transform to Equation (46), we obtain the solution of Equation (44):

Besides, from here we can see that p(t, x, y) depends only on x − y, which gives an alternative proof to the corresponding assertion from Remark 2.2.

Now, turn to the problem of finding m 

where | · | is the vector norm in R d . As was demonstrated, e.g., in [19] under condition (48), the solution of the parabolic problem (44) has, for each x, y ∈ Z d , the following asymptotics:

where

is a constant depending on the lattice dimension. For a more detailed description of asymptotics (49), including the form of the reminder term, see [12] . Let us now apply the above reasoning to Equations (31) and (32). Note that in the case where the migration operators L 1 and L 2 defined by Equation (3) coincide, i.e., L 1 = L 2 , we can refine the representation of Equation (39) for λ 1 (θ) and λ 2 (θ) by using Equation (40), which yields:

Replacing a(θ), b, c and d(θ) in Equation (51) by their values given by Equations (27)-(30) we obtain the following representations of C 1 and C 2 :

Let us denote

Then, in the case when κ 1 = κ 2 = κ, b = 0, c > 0 (or b > 0, c = 0) due to Equations (27) and (30), the following relation holds:

that is, the difference a(θ) − d(θ) does not depend on θ. This means that either a(θ) − d(θ) = 0 for all values of θ or a(θ) − d(θ) = 0 also for all values of θ. Moreover,

Consequently, in the case r 1 = r 2 we have not only that

Thus, from Equations (41)-(43) we have for t → ∞ (we prefer to consider the case b = c = 0 separately from other cases):

In this section, we will study the behavior of the second moments of the number of subpopulations.

To do this, we will essentially use the technique developed in the previous section, so we will omit some technical details.

Let us denote m

x, y) and let the estimate Equation (14) be true. Our goal in this section is to obtain differential equations for m (2) ij (t, x, y), i, j = 1, 2, which are similar to those obtained for m (1) ij (t, x, y), i, j = 1, 2, in Section 3.1. By taking the partial derivatives of the functions Φ i (t, x, y; z) in Equation (11) over z 1 and z 2 we can get the following equations:

Then, by fixing z 1 = z 2 = 1 in the last equation, we obtain

Now, by differentiating both sides of Equation (12) from Lemma 2.1 over z j twice, we obtain:

x, y; z)) .

Taking here z = (1, 1) and using the notation

we obtain the following representation of the left-hand side of Equation (52):

while the right-hand side of the same equation equals to:

By equating Equation (53) with (54) we get

1j (t, x, y)] 2 + lm 

Finally, to obtain the differential equations for the second moments m (2) ij (t, x, y), we add the term ∂ t m (1) ij (t, x, y) to each side of Equation (55). Then, we substitute the term m (1) ij (t, x, y) on the right side of the obtained expression by its representation (15) . Then, the left side of the resulting equation takes the form ∂ t m (2) ij (t, x, y), while the right side is equal to

2j (t, x, y)

1j (t, x, y) + lm (2) 2j (t, x, y)

Thus, we have proved the following lemma. ij (t, x, y), i, j = 1, 2, satisfy the differential equations

1j (t, x, y)

1j (t, x, y)] 2 + lm 

with the initial condition m

ij (0, x, y) = δ i (j)δ x (y).

Remark 4.1. Similar considerations as in Remark 3.1 show that in this case Equation (57) can be treated as a linear differential equation in a Banach space whose right-hand side (for each t and y) is a linear bounded operator acting in any of the spaces l p (Z d ), p ≥ 1. Therefore, for the same reasons as in Remark 3.1, we obtain that m (2) ij (t, x, y) (for each t and y) as a function of the variable x belongs to each of the spaces l p (Z d ), p ≥ 1, and is thus bounded.

So, we have obtained the differential equations for the second moments m (2) ij (t, x, y). In the next section, we will find the solutions for the obtained equations.

In this section, (as in Section 3.2), we will consider the equations for m 2j (t, x, y)), j = 1, 2. Then, using the notation (27)-(30) from Section 3, we obtain ∂ m (2) 1j (t, θ, y) ∂t = a(θ) m

1j (t, θ, y) + bm (2) 2j (t, θ, y) + f (j)

∂ m (2) 2j (t, θ, y) ∂t = c m (2) 1j (t, θ, y) + d(θ)m

2j (t, θ, y) + f 

where m (2) ij (0, θ, y) = δ i (j)e i(θ,y) , i = 1, 2,

and

1j (t, θ, y)

2j (t, θ, y) .

Here,

i.e., (F * G)(t, θ, y) is the convolution of the functions F (t, θ, y) and G(t, θ, y) with respect to the variable θ.

In what follows we will need the explicit form of the solution of the following linear differential equation:

This solution can be readily obtained by the method of variation of parameters, also known as the method of variation of constants:

Case b = c = 0. Here, the functions f ∂ m (2) 1j (t, θ, y) ∂t = a(θ) m The solution of the first equation due to Formula (63) is clearly as follows:

1j (t, θ, y) = m (2) 1j (0, θ, y)e a(θ)t = δ 1 (j)e i(θ,y) e a(θ)t , where the first equality follows from Equation (63), whereas the second equality follows from Equation (61).

To solve the second equation, we again use Formula (63) assuming x(t) = m (2) 2j (t, θ, y), k = d(θ) and f (t) = c m (2) 1j (t, θ, y) + f Then, Equations (59) also take the "triangle" form ∂ m (2) 1j (t, θ, y) ∂t = a(θ) m (2) 1j (t, θ, y) + b m (2) 2j (t, θ, y) + f 

2j (t, θ, y).

2j (0, θ, y)e d(θ)t = δ 2 (j)e i(θ,y) e d(θ)t . where again the first equality follows from Equation (63) whereas the second equality follows from Equation (61).

To solve the first equation we apply the formula (63) with x(t) = m 

12 (0, θ, y) = 0; m Case b > 0, c > 0. To address this case, we first recall the explicit form of the solution to the following linear differential equation:

where A is a matrix (in our problem A is a 2 × 2 matrix) with time-independent (constant) entries and f (t) is a column-vector function. The solution of Equation (64) can be easily obtained by the method of variation of parameters, see, e.g., [7] or any other textbook on linear differential equations:

where the matrix-function U (t) is the so-called "fundamental solution" of Equation (64). Let U (t) = exp{At}. Then It is known [7] that U (t) can be expressed as U (t) = exp{At}. However, for us, the following representation for U (t) will be more useful:

where the vector-functions

are solutions of the homogeneous system

satisfying the initial conditions, respectively,

The components u ij (t) of the solutions u 1 (t) and u 2 (t) can be computed by using calculations from Remark 3.2 (case b > 0, c > 0). By doing the needed computations, we obtain:

where λ 1,2 (θ) are specified by Equation (39). Therefore, from Equations (61) and (65) we obtain the following solutions for Equations (59) and (60):

2 (s, θ, y) e λ2(θ)(t−s) ds

2 (s, θ, y) e λ1(θ)(t−s) ds

2 (s, θ, y) e λ2(θ)(t−s) ds

where the functions f 

In this section, we consider BRWs with two types of particles satisfying the condition that the particle reproduction law at each lattice point is described by an irreducible critical two-type branching process and that the underlying random walks have finite variances of the jumps. We show that for particles of both types with the underlying recurrent random walks on Z d , a phenomenon of clustering of particles can be observed over long times, implying that the majority of particles are concentrated in some particular areas. We generalize the study started in [1] for BRW with one type of particles.

In this section, based on the results for two-type critical branching processes we show that the probability of degeneracy of the subpopulation tends to 1 for the underlying recurrent random walk. We also show that, at the same time, subpopulations that are not degenerate exhibit linear growth in t at infinity. Let us introduce some notation. Denote by D = (d ij ) the matrix with the elements

where F i (z 1 , z 2 ) are the generating functions defined in Equation (2). We also define the densities of second factorial moments of F i (z 1 , z 2 ) (cf. Equation (4) in [15] (Ch. 4, § 7)) as

, i, j, k = 1, 2, and assume that condition (14) holds, so that d ij and b 

where the sum on the right-hand side is finite due to Remark 3.1. Then, the matrix D(t, x) := (d ij (t, x)) with elements

ij (t, x, y), i, j = 1, 2, is well-defined, i.e., its elements d ij (t, x) are finite. Moreover, the relations (66) prove that the quantity d ij (t, x) does not indeed depend on the spatial coordinate x, i.e.,

Then, according to [15] (Ch. 4, § 7, Th. 5), we have

where r is the Perron root (see [15] (Ch. 4, § 5, Def. 6)) of the matrix D and r 1 is some quantity satisfying r 1 < r. We denote by

the left and right eigenvectors, respectively, corresponding to the eigenvalue r of D. 

be the number of particles in a subpopulation at time t generated by a particle of the i-th type provided that at the initial moment of time the particle was at the point x.

Remark 5.1. Evaluate the quantity n i (t, x), i = 1, 2, at the time moment t + dt. Let n ij (t, x), j = 1, 2 be the number of offsprings of type j generated by a single particle of type i, which at the time moment t = 0 was located at the point x ∈ Z d , so that n i (t, x) = n i1 (t, x) + n i2 (t, x). Assume G i (t, x; z) = Ez ni(t,x) . Then, by using the Kolmogorov forward equation, we obtain the following relations:

From these relations, it can be seen that the behavior of the process n i (t, x) depends only on its "branching component" and the evolution of the process coincides with the evolution of the branching process with continuous time treated in [15] . For this reason, we apply the results of [15] in the following.

Remark 5.2. Note that from Remark 2.2, we have for all k ∈ Z + :

Recall that the branching process under consideration is assumed to be critical and irreducible. In this case, [15] (Ch. 6, § 3, Th. 4) implies that the probability of non-degeneration of a subpopulation has the following asymptotic behavior for all x ∈ Z d as t → ∞:

where c i is a constant. Thus, the probability of degeneration P n i (t, x) = 0 of the subpopulation tends to 1 for all x ∈ Z d as t → ∞. Now, we will estimate the conditional mathematical expectation

which is the main object of the study in Section 5.1. By the definition of conditional expectation we have for all

where I{A} is the indicator of the set A. Note that from Equation (69) it follows that

x, y) = 0 for j = 1, 2.

Then, with the usage of the formula of total probability we have

Thus, from Equation (67) we obtain

At the same time, substituting r = 0 in Equation (68) we obtain d ij (t) = u i v j + o(1), whence, denoting C ij := uivj ci = const and using Equation (70) we get

In virtue of Equation (70) we have that, in the case when t → ∞, the probability of degeneration of the subpopulation P n i (t, x) = 0 tends to 1. At the same time, due to Equation (71) those subpopulations that are not degenerate have a linear growth in t at infinity.

In this section, we study the effect of clustering at each point for an irreducible critical branching process under the condition that the tail of a random walk is superexponentially light, i.e., for each λ ∈ R d , i = 1, 2, the following condition holds:

By p i (t, x, y), we denote the transition probability of the random walk on Z d defined by L i , i = 1, 2, see Equation (3) . From [9] (Ch. 3, § 2) it follows that p i (t, x, y) is the solution of the Cauchy problem ∂p i (t, x, y) ∂t = (L i p i (t, ·, y))(x), p i (0, x, y) = δ x (y).

Note that here, as was shown in Remark 3.5, p i (t, x, y) = p i (t, x − y, 0) = p i (t, 0, y − x), which follows from the property of spatial homogeneity of the process under consideration. Denote y − x = s, then from [13] (Eq. (4.7)) we have for s = O( √ t) the following equality:

where

For d = 1, the one-dimensional matrix B i is as follows:

Consider the probability that a particle located at point 0 ∈ Z d will jump no further than a distance C √ t, C > 0 is some constant. Then,

Note that in the last equality, under the sign of integral is the function f (τ ) such that

Function f (τ ) is the density function of a random variable, which has normal distribution with mean 0 and variance b i t. Therefore, choosing appropriate constant C > 0, we can make the quantity

to be arbitrarily close to 1. Hence, for every ε > 0 there exists C > 0 such that

Thus, as t → ∞ a particle will move at a distance of no more than C √ t with probability arbitrarily close to 1.

Turn to the case when Z d , d = 2. Then, from Equation (72) takes the form

As in the previous case, consider the probability that a particle located at point 0 ∈ Z d will jump no further than a distance C √ t. Consequently,

Note that in the last equality, under the sign of integral is the function f (τ 1 , τ 2 ) such that

is the density function of a random variable which has two-dimensional normal distribution with mean vector (0, 0) and covariance matrix B i t.

Therefore, choosing an appropriate constant C we can get

arbitrarily close to 1. Hence, due to Equation (74) for every ε > 0, there exists C > 0, such that

Thus, as t → ∞, a particle will move with probability arbitrarily close to 1 over a distance no greater than C √ t. This result for the lattice dimension d = 2 is similar to that obtained in Equation (73) for the lattice dimension d = 1.

Let us now consider the situation where there is one particle of type i at each point at the initial time. We denote the set of odd positive integers by N 1 and the set of even positive integers by N 2 . Let us consider a given particle at time t and all its progenitors up to the initial time.

We consider a particle at time t and a sequence (evolutionary lineage ) K consisting of all its m progenitors (from the initial particle to the immediate parent) and the particle itself, K = (k 1 , . . . , k m , k m+1 ), m ≥ 0. If m > 0 (i.e., the particle is not included in the set of initial particles on the lattice), then we select s from the sequence of indices [2, . . . , m + 1] such that type(k s ) = type(k s−1 ), where type(k) denotes the type of the particle k. We denote the sequence of selected indices by S = (s 1 , . . . , s n ), n ≤ m + 1. If the sequence S turns out to be empty (i.e., no type changes were observed in the evolutionary lineage considered), then we add to it the index s 1 := m + 1, then n = 1. We denote by h(k) the lifetime of the particle k and construct the sequence τ = (τ 1 , . . . , τ n ), where τ 1 =

Assuming that the evolutionary lineage started with a particle of type 1, we obtain that in the time intervals (0, τ 1 ), (τ 2 , τ 3 ), . . . the particles of this evolution lineage walk on the lattice under the action L 1 and on the time intervals (τ 1 , τ 2 ), (τ 3 , τ 4 ), . . . under the action of L 2 .

We denote by p(τ, x, y) the probability for a particle to move from a point x to y on Z d in time τ . Due to the Kolmogorov-Chapman equation, for n ≥ 2 we obtain

where s(i) = 1 for i ∈ N 1 and s(i) = 2 for i ∈ N 2 . This representation will be needed in the following lemma.

Lemma 5.1. Let t 1 := i∈N1 τ i be the total time spent by a particle in the first state and t 2 := i∈N2 τ i be the total time spent by the same particle in the second state. Then,

Proof. Let us show first that p(τ, x, y) = p((τ 1 , . . . , τ n−2 + τ n , τ n−1 ), x, y).

For the proof, due to Equation (75), it is enough to consider the following sequence of relations:

× p s(n−1) (τ n−1 , x n−2 , y) = p((τ 1 , . . . , τ n−2 + τ n , τ n−1 ), x, y).

Consecutively, applying the Formula (77) n − 2 times, we obtain p(τ, x, y) = p((τ 1 , . . . , τ n−2 + τ n , τ n−1 ), x, y) = p((τ 1 , . . . , τ n−3 + τ n−1 , τ n−2 + τ n ), x, y) = . . . = p((τ 1 + τ 3 + . . . + τ 2[(n−1)/2]+1 , τ 2 + τ 4 + . . . + τ 2[n/2] ), x, y), whence the assertion of the lemma follows. Now, we will apply Equation (76) from Lemma 5.1 for understanding how far the particles can go from the initial position of their initial progenitor by some time t when t → ∞.

In case of t → ∞ and t i → ∞, i = 1, 2, for every ε > 0, there exists C i > 0 such that

In case of t → ∞ and t i < C, i = 1, 2, and for every ε > 0, we have:

Thus, for t → ∞, ∀ε > 0, ∀τ : i τ i = t,

For d = 1, the distance between the start points of subpopulations that did not degenerate by the time t → ∞ has a geometric distribution with an average value of t ci + o(t) Equation (70) and non-degenerate subpopulations have particles at a distance from the initial particle of the order of no more than √ t with a probability arbitrarily close to 1, see Equation (73). Thus, particle clusters with length of order √ t are separated by empty intervals with length of order t.

Let us turn to the case d = 2. Choose two functions, ν(t) and f (t), such that ν(t) → ∞ and

. Now, consider the square of the lattice with side tν(t)e ν(t) ci and divide it into cells with side tν(t) ci ; then the number of cells will be e ν(t) . We call a cell degenerate at time t if it does not contain the starting points of populations that do not degenerate at time t. Then, the probability that for t → ∞ all subpopulations of cell degenerate is

with some constant C > 0. The probability of the existence of a cell, whose all subpopulations of the initial particles are degenerated, is

Non-degenerate subpopulations have particles at a distance from the initial particle of the order of no more than √ t tν(t) ci . Therefore, by the time t → ∞, we get particle-free circles with a radius of the order of tν(t) ci at a distance of the order of tν(t)e ν(t) ci . Thus, we have proved that both in the case of dimension d = 1 and d = 2, the effect of clustering of particle subpopulations takes place.

One of the assumptions of the model we presented in Section 2 was the fact that particles cannot change their type over time (see Remark 2.1). Here, we will consider an example where particles of the first type can become particles of the second type. This example can describe the distribution of a virus.

In Section 6.1, we will describe a new model of BRW with two types of particles, where the particles can change their types. We will use the designations from Section 2. In Section 6.2, we study the first moments for the number of particles of type i = 1, 2 at each lattice point. In Section 6.3, we obtain the solutions for the second moment for the number of particles of the first type at each lattice point in the more general case. In Section 6.4, we study the effect of intermittency in the simplest case for the number of particles of the first type. In Section 6.5, we obtain the differential equation for the second moment for the number of particles of the second type at each lattice point and find its asymptotic behavior as t → ∞.

Consider a new model of BRW with two types of particles. Here, we will study the behavior of the processes N i (t, x), i = 1, 2 defined in (1). We call the particles of the first type infected and the particles of the second type particles with immunity. Let us denote by r the intensity to build up immunity for an infected particle during the small time dt. This means that the particle can change type with probability r dt + o(dt). Moreover, we assume that there was only one infected particle on the lattice at time t = 0. Without limiting generality, we can assume that this initial particle was at the origin. Then, N 1 (0, x) = δ 0 (x), N 2 (0, x) ≡ 0 for all x ∈ Z d . Let b n , n 2 be the intensity to infect n − 1 new particles. Here we assume that there are enough healthy particles at each point of the lattice to get sick. We are also interested in studying the moments of the number of particles of both types. In the previous notation, we say that b n = β 1 (n, 0), β 2 (k, l) ≡ 0 for all k, l : k + l 2.

In what follows, due to the fact that particles can change their types, we use the forward Kolmogorov equations approach to obtain the differential equations for the moments of N i (t, x), i = 1, 2. The derivation of the forward Kolmogorov equations is based on the following representation:

where ξ(dt, x) and ψ(dt, x) are discrete random variables with the following distributions:

We are going to study the first two moments for the random variables N i (t, x), i = 1, 2. In the next section, we pay attention to the first moments.

In this section, we consider the first moments for N i (t, x), i = 1, 2. We obtain the differential equations for them and find their explicit solutions in terms of the Fourier transform (23). We also obtain their asymptotic behavior in certain cases. Define the first moments R i (t, x) := EN i (t, x).

Note that R i (0, x) = δ 1 (i)δ 0 (x). Let F t be the sigma-algebra of events up to and including t. Note that ξ(dt, x) (and ψ(dt, y)) and F t are independent.

Derive the differential equations for these functions:

Then, as dt → 0 the differential equation for R 1 (t, x) is

The same technique helps to find the differential equation for R 2 (t, x):

From this we get as dt → 0:

Firstly, solve equation Equation (78). Write again the equation for this function:

To solve the equation apply the discrete Fourier transform Equation (23). Then,

The solution has the form: Therefore,

We find out the solution of Equation (79). Substituting the solution for R 1 (t, x) from Equation (80) in Equation (79), we obtain:

To obtain the solution of differential Equation (81), we consider two cases:

Case β − µ 1 − r = −µ 2 . Apply the discrete Fourier transform (23) and the variation of constants methods to solve the Equation (81). If κ 1 a 1 (θ) = κ 2 a 2 (θ), then R 2 (t, θ) = rte (κ2 a2(θ)−µ2)t .

Consequently,

If κ 1 a 1 (θ) − κ 2 a 2 (θ) =: d > 0, then

Then,

Here, the solution of differential Equation (81) has the form:

Remark 6.2. Find out the asymptotic behavior for the first moments R i (t, x), i = 1, 2 in a particular case. Assume that generators L i , i = 1, 2 defined in Equation (3) 

where γ d is specified by Equation (50). With the usage of (82) we get for R i (t, θ), i = 1, 2 obtained above we have, for each x ∈ Z d , as t → ∞:

when β − µ 1 − r = 0 and µ 2 = 0, and

when A = β − µ 1 − r = 0 and µ 2 = 0 So, we have found the first moments for both types of a particle. In the next sections, we are going to get the explicit form of the second moments.

Here, we will find out the asymptotic behavior for the second moments of the random variable  N 1 (t, x) .

To derive the second moment for N 1 (t, x), we consider a more general problem. Let N 1 (t, x, y) be the number of particles of the first type at time t at point y ∈ Z d generated by the single particle of the first type located at time t = 0 at the site x ∈ Z d . The initial condition for N 1 (t, x, y) is N 1 (0, x, y) = δ x (y). If we use the designations from Section 6.1, if N 1 (t, x) = N 1 (t, 0, x) . Here, we will use the method of backward Kolmogorov equations. Define the generating function for the random variable N 1 (t, x, y) as

where z ∈ R, t 0. From this, we get the following lemma.

Lemma 6.1. The generating function F (t, x, y; z) specified by Equation (83) satisfies the differential equation:

x, y; z));

F (0, x, y; z) = e −zδx(y) .

(84)

Proof. Consider the generating function F (t, x, y; z) at the time moment t + dt. Then,

x, y; z) + (L 1,x F (t, ·, y; z))(x) dt + f F (t, x, y; z) dt

where f (s) = µ 1 + s(−µ 1 − n 2 b n ) + n 2 b n s n and

Then,

x, y; z) = (L 1,x F (t, ·, y; z))(x) dt + f F (t, x, y; z) dt

Therefore, as dt → 0 ∂F (t, x, y; z) ∂t = (L 1,x F (t, ·, y; z))(x) + f F (t, x, y; z) + r(1 − F (t, x, y; z)).

The initial condition for the latter equation follows from Equation (83):

F (0, x, y; z) = Ee −zN1(0,x,y) = Ee −zδx(y) = e −zδx(y) .

In what follows, we are going to study the effect of intermittency in one area of x and y when |y − x| = O( √ t) as t → ∞. Denote by p(t, x, y) the solution of the following Cauchy problem ∂p(t, x, y) ∂t = (L 1 p(t, ·, y))(x), p(0, x, y) = δ x (y).

Then, the representation from Equation (85) has the form

Using Duhamel's principle and Equation (88), from Equation (86) we obtain

Remark 6.6. In cases when underlying random walk has infinite variance of jumps, so that relation (see definition in [8] (Equation (1.2))) when ν := β − µ 1 − r > 0, the results were obtained in [8] (Th.1.2).

Now, consider the case where ν > 0 and the underlying random walk has superexponentially light tails of a random walk such that for all λ ∈ R d z∈Z d e (λ,z) a 1 (z) < ∞.

z k z j a 1 (z), k, j = 1, . . . , d.

Note that for all x, y ∈ Z d , t > 0

Then, for t → ∞ p(t, 0, 0) = C t d/2 + o(t −d/2 ) and 0 < p(t, 0, 0) < 1, we have

For t → ∞, |x − y| = O( √ t) and ν > 0 from Equation (88) in the numerator of the last summand. Then, we can continue the estimation = β (2) p(t, x, y) t 0 e νs p(s, 0, 0)ds e νt p 2 (t, x, y)

t 0 e νs p(s, 0, 0)ds e νt p(t, x, y)

Thus, the random variable N 1 (t, x, y) is non-intermittent for |x − y| = O( √ t).

In this section, we will set up the differential equation for the second order correlation function for N 2 (t, x) and determine the asymptotic behavior in the special case. Define the following correlation function for the second type particles

To obtain the differential equation for R 22 (t, x, y), we consider this random variable at time moment t + dt. Unlike the differential equation for N 1 (t, x), here we consider two cases: x = y and x = y.

Firstly, consider the case when x = y. Here, we have for R 22 (t + dt, x, x):

Then, as dt → 0 we obtain the differential equation for R 22 (t, x, x):

Now, consider the case when x = y. Then,

x) + ψ(dt, x))(N 2 (t, y) + ψ(dt, y))|F t ] = R 22 (t, x, y)

x, y) dt + r R 12 (t, x, y)

Therefore, as dt → 0 we have the differential equation for R 22 (t, x, y) when x = y:

Above, we have defined function R 12 (t, x, y) = EN 1 (t, x)N 2 (t, y). This is unknown function, consequently, we need to obtain the differential for this function. Using the same technique as for R 22 (t, x, y) we get:

− δ x (y)rR 1 (t, x) dt + (L 1,x R 12 (t, ·, y))(x) + (L 2,y R 12 (t, ·, y))(x) dt

Then, as dt → 0 the differential equation for R 12 (t, x, y) is ∂R 12 (t, x, y) ∂t = (L 1,x R 12 (t, ·, y))(x) + (L 2,y R 12 (t, ·, y))(x)

− δ x (y)rR 1 (t, x);

To get the behavior of this function, we also need to have the differential equation for it. Then, consider R 11 (t + dt, x, y):

If dt → 0, we obtain R 11 (t, x, y): ∂R 11 (t, x, y) ∂t = (L 1,x R 11 (t, ·, y))(x) + (L 1,y R 11 (t, ·, y))(x)

In the future calculus, we need the following lemma.

Lemma 6.2. Let G(t, x, y) := R 1 (t, x)R 2 (t, y) and K(t, x, y) = R 1 (t, x)R 1 (t, y). Then, functions G(t, x, y) and K(t, x, y) satisfy the following differential equations ∂G(t, x, y) ∂t = (L 1,x G(t, ·, y))(x) + (L 2,y G(t, ·, y))(x) + (β − µ 1 − r − µ 2 )G(t, x, y) + rK(t, x, y); G(0, x, y) = 0.

∂K(t, x, y) ∂t = (L 1,x K(t, ·, y))(x) + (L 1,y K(t, ·, y))(x) + 2(β − µ 1 − r)K(t, x, y);

Proof. Note that

Then, with the usage of Equations (78) and (79), we have G(t, x, y) ∂t = R 2 (t, y) ∂R 1 (t, x) ∂t + R 1 (t, x) ∂R 2 (t, y) ∂t = R 2 (t, y) (β − µ 1 − r)R 1 (t, x) + (L 1 R 1 (t, ·))(x) + R 1 (t, x) (L 2 R 2 (t, ·))(y) − µ 2 R 2 (t, y) + rR 1 (t, y) = (L 1,x G(t, ·, y))(x) + (L 2,y G(t, ·, y))(x) + (β − µ 1 − r − µ 2 )G(t, x, y) + rK(t, x, y).

Similarly, we can obtain the differential equation for K(t, x, y). Notice that K(t, x, y) ∂t = R 1 (t, y) ∂R 1 (t, x) ∂t + R 1 (t, x) ∂R 1 (t, y) ∂t .

Consequently, using Equation (78), we get K(t, x, y) ∂t = R 1 (t, y) ∂R 1 (t, x) ∂t + R 1 (t, x) ∂R 1 (t, y) ∂t = R 1 (t, y) (β − µ 1 − r)R 1 (t, x) + (L 1 R 1 (t, ·))(x) + R 1 (t, x) (β − µ 1 − r)R 1 (t, y) + (L 1 R 1 (t, ·))(y) = (L 1,x K(t, ·, y))(x) + (L 1,y K(t, ·, y))(x) + 2(β − µ 1 − r)K(t, x, y).

The initial conditions follow from G(0, x, y) = R 1 (0, x)R 2 (0, y) = 0, K(0, x, y) = R 1 (0, x)R 1 (0, y) = δ 0 (x)δ 0 (y).

Find out the asymptotic behavior of the second moment R 22 (t, x, y) in the particular case where we consider Z d . Let the generators of the random walk be identical for both types of particles, so that L 1 = L 2 =: L. Assume that the random walk with generator L has finite variance of jumps (see (48)). For A = β − µ 1 − r > 0 and µ 2 = 0, it was found in Section 6.2 that, as t → ∞, for each x ∈ Z d ,

re At At d/2 . Note that from Lemma 6.2 and Equation (91), we obtain ∂(K(t, x, y) − R 11 (t, x, y)) ∂t = (L 1,x (K(t, ·, y) − R 11 (t, ·, y)))(x) + (L 1,y (K(t, ·, y) − R 11 (t, ·, y)))(x) + 2(β − µ 1 − r)(K(t, x, y) − R 11 (t, x, y)) + F (R 1 (t, x)),

where F (R 1 (t, x)) is the function which linearly depends on R 1 (t, x). Then, from the representation of R 1 (t, x) as t → ∞ we have, for each x ∈ Z d ,

where C(·) function which can be obtained from Equation (91). Then, as t → ∞ K(t, x, y) − R 11 (t, x, y) = Ce At t d/2 , C − constant.

The same technique helps us to find out that G(t, x, y) − R 12 (t, x, y) = C e At t d/2 , t → ∞ C − constant.

Due to the homogeneity of space, we can consider the following values: R 22 (t, u) := R 22 (t, 0, y− x).

Using the above results, rewrite the equation for the second correlation function:

∂R 22 (t, u) ∂t = 2(LR 22 (t, ·))(u) + r 2 γ 2 1 (u)e 2At t d + δ 0 (u) rγ 1 (u)e At t d/2 − a 2 (u) rγ 2 (u)e At √ 2πt ;

R 22 (0, u) = 0.

Divide the last equation into two equations and find the solutions R re At √ 2πt − a 2 (u) re At t d/2 , R

22 (0, u) = 0.

Then, we will have R 22 (t, u) = R Here, for a fixed space coordinate, we do not have the intermittency effect:

for each x ∈ Z d .

This section is devoted to process modeling using the Python programming language. We consider the state of the system as an array whose elements are lists of the form [i, x, t 1 , t 2 ], where i characterizes the type of a particle and x = (x 1 , . . . , x d ) is its spatial coordinate, t 1 is the time of its entry into a given position (it was born at x at this time or jumped at x at this time), t 2 is the time of its exit from this position (it died at x at this time or jumped out of x at this time). Recall that we perceive all events related to the reproduction of offspring, including the degeneration from one type to another, as the death of the parent particle with the production of k descendants of the first type and l of the second type. Initialization. We set the characteristics of BRW, i = 1, 2: 

In this work, we considered a continuous-time branching random walk with two types of particles.

The main results were devoted to the study of the limiting behavior of the moments of subpopulations generated by a single particle of each type. In particular, in Section 3, we have derived the differential equations for the first moments of subpopulations and found their solutions, which allows us to find exact expressions for their asymptotics. Similar results for the second moments were obtained in Section 4. In Section 5, we have shown that for particles of both types with the underlying recurrent random walks on Z d a phenomenon of clustering of particles can be observed over long times, which means that the majority of particles are concentrated in some particular areas. The obtained results were then applied in Section 6 to study epidemic propagation. In this model, we considered two types of particles: infected and immunity generated. At the beginning, there is an infected particle that can infect others. Here, for the local number of particles of each type at a lattice point, we study the moments and their limiting behavior. Additionally, the effect of intermittency of the infected particles was studied for a supercritical branching process at each lattice point. To demonstrate the effect of limit clustering for the epidemiological model, we provided the results of a numerical simulation in Section 7.

We would like to emphasize that the present work is primarily devoted to a theoretical analysis of the effects that occur in branching random walks with two types of particles. We have tried to illustrate the obtained theoretical results by a numerical simulation. However, this simulation should not be considered as a full-fledged numerical analysis of the considered situation. Therefore, the question of developing real programs (in Python, R or any other language) for the numerical analysis of the considered processes arises quite naturally. In this paper, the authors did not undertake such a task, mainly because the fact that the computational aspects of modeling multidimensional processes are a special science that can hardly be treated professionally and completely in one or a few sections of even such an extensive article as ours. Possibly, further special studies will be devoted to it. 

R is the array consisting of a finite number of lists [i, x, 0, 0] characterizing types i and particle positions x = (x 1 , . . . , x d ) at the initial moment of time

A i = (a i (x, y)) are matrices of the random walk intensities, by which the generators (3) are determined

We select one of the elements [i, x, t 1 , t 2 ] of the array R such that t 2 < T . The particle spends exponential time dt at the current position x, after which it does one of the following: 1. with probability µ i /(κ i + µ i + k+l 2 β i (k, l) + r i ) dies

)/(κ i + µ i + k+l 2 β i (k, l) + r i ) divides into k + l particles

/(κ i + µ i + k+l 2 β i (k, l) + r i ) jumps from position x to position x + z = x, then we append

+ µ 1 + k+l 2 β 1 (k, l) + r) turns into a particle of the second type

After the process is completed, the entire history of particles in different states is located in arrays R and H. To find out the number, type and spatial configuration of particles at time t, we select those elements [i, x, t 1 , t 2 ] of the array H for which t 1 t < t 2 . Simulations. Suppose d = 1 and at time t = 0 there are 300 initial particles on segment [0, 300] ∈ Z. The random walk for the particles of the first type has intensities a 1 (z) = 1/2 for simulation with parameters µ 1 = 0.25, β 1 (2, 0) = 0.125, β 1 (1, 1) = 0.125, µ 2 = 0.375, β 2 (0, 2) = 0.125, β 2 (1, 1) = 0.25, all other birth/death intensities equal to 0. It demonstrates the clustering effect in the case of the critical branching law described in Section 5

Structure of the particle population for a branching random walk with a critical reproduction law

Stochastic differential equation with jumps for multi-type continuous state and continuous time branching processes with immigration, ALEA Lat

On aggregation of multitype Galton-Watson branching processes with immigration

Extinction in lower hessenberg branching processes with countably many types

The probabilities of extinction in a branching random walk on a strip

Spread of a catalytic branching random walk on a multidimensional lattice

Sbornik zadach po differentsial'nym upravleniyam. uchebnoe posobie, Moskva: Nauka

Intermittency for branching walks with heavy tails

Vvedenie v teoriyu sluchaȋnykh protsessov, Izdat

On two-type branching random walks and their applications for genetic modelling

Branching random walks with immigration. Lyapunov stability, Markov Process

Limit theorems for the Green function of the lattice Laplacian under large deviations for a random walk

Large deviations for a symmetric branching random walk on a multidimensional lattice

Spatial models of population processes, Modern problems of stochastic analysis and statistics

Sevast'yanov, Vetvyashchiesya protsessy, Izdat

Counterexamples in probability

The survival probability for a class of multitype subcritical branching processes in random environment

Multi-type subcritical branching processes in a random environment

Branching random walks in a heterogeneous environment

Spectral properties of evolutionary operators in branching random walk models

On conditions for a probability distribution to be uniquely determined by its moments

Later, we also will use the following notation:(L i,y Ψ(t, ·, y))(x) = κ i v =0 a i (v)[Ψ(t, x, y + v) − Ψ(t, x, y)], i = 1, 2.Let M 1 (t, x, y) = EN 1 (t, x, y) be the first moment of N 1 (t, x, y). Note thatThen, from Equation (84) we can derive the differential equation for the first moment taking the partial derivative of both sides of (84). Omitting the calculus we obtainAs above, with the usage of the discrete Fourier transform (23), we can find the solution for this equation:Then, using the inverse Fourier transform (24), we will obtain:Remark 6.3. Notice that from the obtained representation for M 1 (t, x, y), we get that the first moment M 1 (t, x, y) is a function that depends on the difference of the considered sites on the lattice, so thatFrom Equation (84) (by taking partial derivative over parameter z twice and substituting z = 0) we can derive the differential equation for the second moment, which can be defined asConsequently, omitting the calculuswhereApply the discrete Fourier transform (23) to this equation we obtain:Remark 6.4. As in Section 6.2, we consider the asymptotic behavior of M 2 (t, x, y) in case when random walk for the particles of the first type has finite variance of jumps, so that v a 1 (v)|v| 2 < ∞.

Let M 2,h (t, x, y) be the solution of homogeneous equation and M 2,p (t, x, y) be the particular one. Then, M 2 (t, x, y) = M 2,h (t, x, y) + M 2,p (t, x, y). Assume that β − µ 1 − r = 0. Then, previous equation takes the formFrom Remark (3.5) and Equation (82), we get thatThen, as t → ∞, f (t) is the solution of Equation (87) and M 2,Similarly, for d 2, we can find that, for each x, y ∈ Z d ,Consequently, as M 2 (t, x, y) = M 2,h (t, x, y) + M 2,p (t, x, y) we obtain that, for each x, y ∈ Z d ,We have thus obtained the asymptotic behavior of the first two moments of the random variable N 1 (t, x, y). In the next section, we will examine the effect of intermittency (see definition in the next section ) for the random variable N 1 (t, x, y) using M i (t, x, y), i = 1, 2.

In Sections 6.2 and 6.3, we obtained the solutions for the first two moments for the random variable N 1 (t, x, y). Here, we will study the effect of intermittency in the simplest case for the number of particles of the first type. Introduce the following definition (see, for example, [8] ). where x ∈ Ω(t), and Ω(t) is a non-decreasing family of sets.Remark 6.5. We are going to consider the effect of intermittency for random variable N 1 (t, x, y). In our designations N 1 (t, y − x) = N 1 (t, 0, y − x).