A general mathematical method for predicting spatio-temporal correlations emerging from agent-based models

Agent-based models are used to study complex phenomena in many fields of science. While simulating agent-based models is often straightforward, predicting their behaviour mathematically has remained a key challenge. Recently developed mathematical methods allow the prediction of the emerging spatial patterns for a general class of agent-based models, whereas the prediction of spatio-temporal pattern has been thus far achieved only for special cases. We present a general and mathematically rigorous methodology that allows deriving the spatio-temporal correlation structure for a general class of individual-based models. To do so, we define an auxiliary model, in which each agent type of the primary model expands to three types, called the original, the past and the new agents. In this way, the auxiliary model keeps track of both the initial and current state of the primary model, and hence the spatio-temporal correlations of the primary model can be derived from the spatial correlations of the auxiliary model. We illustrate the agreement between analytical predictions and agent-based simulations using two example models from theoretical ecology. In particular, we show that the methodology is able to correctly predict the dynamical behaviour of a host–parasite model that shows spatially localized oscillations.

Generator L R,C,P describes the following random event. Let, at a moment of time t ≥ 0 the system is given by a configuration γ t . The event is that, within a small time-interval [t, t + δt], a group of reactants {x i 1 1 , . . . , x iR R } disappears from the configuration and a group of products {z l 1 1 , . . . , z lP P } will become a part of the configuration, so that z 1 ∈ Λ 1 , . . . , z P ∈ Λ P for some disjoint bounded subsets Λ 1 , . . . , Λ P of R d . Thus, γ t+δt = γ t \ {x i 1 1 , . . . , x iR R } ∪ {z l 1 1 , . . . , z lP P }. The probability of this event is then δt · δt = 0. The influence on catalysts y j 1 1 , . . . , y jC C reflects the interaction between agents. The catalysts remain unchanged within the event, however, they influence the probability of the event.
We include also the case when, for some m, n, a reactant x m im and a product z n ln are such that x m = z n , whereas i m = l n , i.e. when an agent keeps its position with changing its type only. In this the corresponding integral w.r.t. dz n is omitted. We can also formally treat this as like the function r includes the factor δ(x m − z n ); henceforth δ(x) is the Dirac delta-function.
The dynamics of γ t if defined then through the differential equation: (1.5) which should be satisfied for a large class of functions F .

Spatial correlation functions and cumulants
Definition 1.1. For each i ∈ I N , a function k i t (x) ≥ 0 is said to be the first order spatial correlation function of type i (for the distribution of γ t ), if for any function g 1 (x) ≥ 0, Henceforth, E[·] denotes the expected value of a random quantity (w.r.t. the distribution of γ t ).
The function k i t (x) is also called the density of agents of type i, since, taking g 1 (x) = 1 1 Λ (x) for some bounded subset Λ of R d , where we get from (1.6) that Henceforth, |η| denotes number of points in a finite subset η of R d .
Definition 1.2. For each i, j ∈ I N , a function k i,j t (x 1 , x 2 ) = k j,i t (x 2 , x 1 ) ≥ 0 is said to be the second-order spatial correlation function (between agents of the types i and j), if, for any symmetric function g 2 (x 1 , x 2 ) ≥ 0, (1.7) Remark 1.3. Recall, we assume that agents cannot occupy the same position, hence, for i = j, γ i t and γ j t are disjoint, and thus the restriction x 1 = x 2 in (1.7) is redundant then.
Combining (1.7) with (1.6), we can also write, for all i, j ∈ I N , Henceforth, 1 1 i=j denotes the Kronecker delta (that is 1 if i = j and 0 otherwise). Substituting to (1.8) the symmetric function where Λ 1 , Λ 2 are bounded subsets of R d , we get (1.10) One can also consider the centralized spatial moment that is the expectation of the product of centralized random quantities |γ i t ∩ Λ 1 | − E[|γ i t ∩ Λ 1 |], i ∈ I N (called so because the expectation of each such quantity is 0): Definition 1.4. The function (1.11) is called the second order spatial cumulant between types i and j.
We have hence (1.12) We going to define now a general spatial correlation function.
Remark 1.7. If N = 1, so that i 1 = . . . = i n = s 1 , then we will normally use the notation: k In particular, k (1) t (x) := k s 1 t (x). We can also rewrite then (1.13) as follows,  t (x 1 , . . . , x n ); dx 1 . . . dx n . and the spatial correlation function k (n) t is also called the n-th order spatial factorial moment. Note that k (n) t is a symmetric function.
To define n-th order spatial cumulants, we note that [i 1 , . . . , i n ] is called a multiset, i.e. a collection of n (possibly repeating) elements from I N .
Definition 1.8. We set u i t (x) := k i t (x), i ∈ I N , and define spatial cumulants through the equality Remark 1.9. To see that (1.14) indeed defines spatial cumulant u t for given spatial correlation functions k t , note that the right hand side (r.h.s. henceforth) of (1.14) contains the term with m = 1 which is just u i 1 ,...,in Differential equations for spatial correlation functions can be obtained from (1.5) by using the definition (1.13). Namely, we take as F in (1.5) the integrand in the left hand side (l.h.s. henceforth) of (1.13), i.e. F = F n , where F (γ t ) = F n (γ i 1 t , . . . , γ in t ) = where g n is a symmetric function such that, for some bounded subset Λ of R d , g n (x 1 , . . . , x n ) = 0 if only x m / ∈ Λ for some 1 ≤ m ≤ n. Differentiating both part of (1.13), we will get then from (1.5) that where, recall, F = F (g n ).
The next step is to represent where g m are also symmetric functions depending on g n , and types j 1 , . . . , j m ∈ I N ; they all depend on the particular form of the operator L; in the case of an RCPgenerator, g m and types depend on the rate (1.2). Note that then the summation in m is finite. We will get then, by (1.13), Since F n depends on g n linearly and LF depends on F linearly, we have that L m,n g n := g m depend on g n linearly as well. By considering a dual operator L n,m := ( L m,n ) * , we will get that Since g n was arbitrary, we get then Considering an infinite vector k t of all functions k i 1 ,...,in t indexed by n ≥ 1 and by different multisets [i 1 , . . . , i n ] of types, we can treat the r.h.s. of (1.15) as the action of an infinite matrix L , whose entries are operators L n,m .
Stress that, typically, m can take values larger than n, so that the system of linear differential equations (1.15) is not closed and cannot be solved analytically nor numerically.
The explicit form for the action of L in case of the RCP-generator L, given by (1.3), can be found in [1,Supplementary Note 1].
The differential equations for spatial cumulants can be obtained by substituting (1.14) into (1.15). The equations will have a similar form with, however, nonlinear operators Q n,m . For their explicit form, in the case of L given by (1.4), we also refer to [1, Supplementary Note 1].

Beyond mean-field expansion for spatial dynamics
Equation (1.15) has initial conditions at, say, time t = 0: k i 1 ,...,in 0 (x 1 , . . . , x n ). The important class of such initial conditions are product functions where q i 1 0 , . . . , q in 0 are nonnegative functions on R d . Spatial correlation function (1.16) corresponds to the Poisson distribution of configurations. The characteristic feature of the Poisson distribution is that random numbers are independent for all disjoint bounded subsets Λ 1 , . . . , Λ n of R d ; in particular, all corresponding spatial cumulants of an order more than 1 are equal to 0, cf. (1.12). The Poisson distribution is also called chaotic because of the mentioned independence.
In most cases, however, the solution to (1.15) with the initial condition (1.16) does not have a product structure. The idea of the mean-field approximation (with a small parameter ε > 0) is to find a modification L ε of the Markov operator L in (1.5), such that the solution k i 1 ,...,in would be approximately (up to certain order of ε) equal to a product function. Hence the distribution of γ ε,t would be approximately chaotic, in a certain sense. This is called the propagation of chaos in statistical physics. The realization of the scaling procedure for the RCP-generator is as follows. We assume that r in (1.2) is given through combinations of various kernels of the form (1.18) where v, w ∈ R d , k, m ∈ I N , a k,m ≥ 0 is a function on R d . Here v k , w m are some agents among reactants x i 1 1 , . . . , x iR R , catalysts y j 1 1 , . . . , y jC C or products z l 1 1 , . . . , z lP P . We consider L ε given by (1.4) where r ε,R,C,P has the same structure as r R,C,P , however, the kernels a(v k , w m ), given previously by (1.18), are replaced now by i.e. the scaled kernels have the same full integral but a scaled (expanded) shape.
Next, we consider the initial condition to the corresponding equation (1.17), as follows (1), (1.20) where lim ε→0 o(1) = 0 (in particular, one can consider the initial condition without that o(1) at all). The statement is that then the solution to (1.17) has the property where q i t , i ∈ I N , solve a system of (nonlinear) differential equations where q t is the vector of all q s 1 t , . . . , q s N t ; with certain (nonlinear) mappings H s 1 q , . . . , H s N q . For the exact form of H i q , i ∈ I N , we refer to [1,Supplementary Note 1]. By (1.21), the cumulants of all orders bigger than 1 corresponding to the function k i 1 ,...,in ε,t (x 1 , . . . , x n ), n ≥ 2, through expansion (1.14) are equal to o(1). In particular, cf. (1.11), The term o(1) in (1.21) depends in a non-trivial way on both time t, variables x 1 , . . . , x n and types i 1 , . . . , i n . To partially reveal this dependence, one needs the next term of the expansion. It was shown in [1,Supplementary Note 1], that Here q i t satisfies (1.22) and p i t , g i,j t satisfy certain linear differential equations • p t is a vector of all p i t , i ∈ I N ; and g t is the vector of all g i,j t , i, j ∈ I N ; • mappings H i p [q t ] and H i,j g [q t ] depends on q i t in a nonlinear (in general) way.
One can also get then from (1.23) the following enhancement of (1.21) for n = 2: Space-homogeneous case Consider the special case, where, initially, the density does not depend on space and the pair-correlation is translation invariant, namely: (1.28) Then if the operator L ε has the form (1.4) with r R,C,P in (1.3) replaced by r ε,R,C,P which is a combination of pair-interaction kernels as above, then, for all t ≥ 0, where q i t , i ∈ I N , satisfy the system (1.22) of ordinary differential equation, and (see [1,Supplementary Note 1,formula (241)]) the equation (1.25) can be rewritten in terms of the Fourier transform of functions g i,j t (x), defined by where x · ξ denotes the standard dot-product in R d and i 2 = −1. Namely, g i,j t (ξ) satisfies the following differential equation, for each ξ ∈ R d , Here, similarly to above, the result is a function of ξ; i ∈ I N , in general, nonlinearly. When i, j run over I N , the system of equations (1.30) can be read as a linear nonhomogeneous system of (ordinary) differential equation (considered independently for each value of ξ). Since all q i t , i ∈ I N , are known, one can solve (1.30) explicitly.

Spatiotemporal characteristics 2.1 Spatiotemporal correlations
By (1.13), spatial correlation function k i 1 ,...,in t (x 1 , . . . , x n ) characterizes the probability to find n agents of types i 1 , . . . , i n of the configuration γ t in vicinities of positions x 1 , . . . , x n ∈ R d , respectively. It is naturally also important to characterize the similar probability when agents appear in those vicinities at different moments of time.
We restrict ourselves to two moments of time only: t ≥ 0 and t + ∆t for some ∆t ≥ 0. Moreover, we consider the second order spatiotemporal correlations only. Namely, we are interested to find, for each i, j ∈ I N , a function k i,j t,∆t (x 1 , x 2 ) ≥ 0, such that, for each symmetric It it worth noting that, as we will see below, To obtain k i,j t,∆t , we proceed as follows. Recall, we consider dynamics of For each i ∈ I N , we consider auxiliary dynamics of three configurations γ iO ∆t , γ i+ ∆t , γ i− ∆t , which will be created at the moment t and will have 'own local time' ∆t ≥ 0. Namely • γ iO ∆t contains all agents of the considered system which were present in the system at time t having type i; didn't change their positions nor types within the time interval [t, t+∆t], and hence they are still present in the system at time t + ∆t having type i; • γ i+ ∆t contains all agents of the considered system which appeared in the system within the time interval [t, t + ∆t] having type i; didn't change their positions nor types after that, and hence they are still present at the system at time t + ∆t having type i; • γ i− ∆t contains all agents of the considered system which were present in the system at time t having type i; do not present in the system at time t + ∆t: namely, each of such agent, within the time interval [t, t + ∆t], either disappeared from the system or changed its position and/or type.
We have hence Therefore, our full auxiliary dynamics is i.e. it contains agents of 3N types (recall that a type is just a label). We set (2.6) Recall also that there are two notations to represent configurations, see (1.1), hence, we can also write where, for A ∈ {O, +, −}, Next, by (2.4), we can rewrite the l.h.s. of (2.1) as follows and using (1.7) and Remark 1.3 for (2.5), one can continue Therefore, by (2.1), however, by (2.8), in general, (2.2) holds.

Dynamics of the auxiliary model
To study the dynamics of γ ∆t , we need to consider a modification of (1.5): with an appropriate modification L of L given by (1.4). Namely, we need that • a reactant x iO does not just disappear from γ O ∆t , but changes its type becoming • a reactant x i+ just disappears from γ + ∆t (since then this agent did not exist at time t and will not exist at time t + ∆t); • products may be of the type '+' only; a product z i+ just appears in γ i+ ∆t ; • the event should happen with the same rate as for L, but applied to the union of O-catalysts and +-catalysts (as they all present at the system on the time interval [t, t + ∆t]); • agents of the type '−' do not perform own dynamics (hence they appear because of the transformation from O-reactants only).
Let R, C, P be fixed and L R,C,P be given by (1.3). It means that there are still R reactants, some of them are O-reactants (denote their number by r, so that 0 ≤ r ≤ R), the rest are +-reactants, namely, there are R − r ≥ 0 of +-reactants. Next, there are P products, all are +-products by the above. Finally, C catalysts should be chosen from γ O ∆t ∪ γ + ∆t ; let, similarly to reactants, there be c O-catalysts (0 ≤ c ≤ C) and hence C − c ≥ 0 +-catalysts.
As a result, we will get Naturally, we also set Remark 2.2. It may be also convenient to use another style of writing. Namely, we will interpret now γ ∆t as the union: rather than as the the tuple (2.7). Then, for

Beyond mean-field expansion for spatiotemporal dynamics
Let L ε be given by (1.4) with r R,C,P in (1.3) replaced by (1.19) as it is described in Subsection 1.3. Let, initially, γ ε,0 be distributed so that the corresponding correlation functions has the form (1.20) for certain fixed collection of functions q i 0 , i ∈ I N . Let γ ε,t and γ ε,t+∆t be the corresponding random configurations at times t and t + ∆t, respectively.
We consider the corresponding auxiliary dynamics distributed according to the generator L ε obtained from L ε by an analogy to that done in Subsection 2.2, in particular, by (2.3), . , x n ), ı 1 , . . . , ı n ∈ I N , be the corresponding system of correlation functions.
Recall that, cf. Remark 2.1, for i, j ∈ I N , Space-homogeneous case Consider again the special case where (1.28) holds. By (2.3), the auxiliary (and scaled by ε) dynamics will inherit that property as well: for all ı,  ∈ I N , k ı ε,∆t will not depend on a space coordinate, and k ı,  ε,∆t (x) (and hence u ı,  ε,∆t (x), g ı,  ∆t (x)) will depend on one space coordinate only. Rewriting the formulas above, one gets and where, initially, (2.35) By using (1.30) with i, j ∈ I N replaced by ı,  ∈ I N , we get differential equations for all g ı,  ∆t (ξ) = g iA,jB ∆t (ξ), A, B ∈ {O, +, −}, with coefficients dependent on q ı t . Solving the obtained system of differential equations, one can find g i,j t,∆t (ξ) by (2.32). However, we are interested to simplify the computations by finding a differential equation on g i,j t,∆t (ξ). We are going to formulate now a conjecture which can be verified for various models (in particular, for the considered below). We are going to prove it in a forthcoming paper.
To formulate the conjecture, we consider auxiliary functions on R d : By (1.30) and (1.22), we get that the vector h t = (h i,j t ) i,j∈I N satisfies a nonhomogeneous system of linear differential equations: where, similarly to above, A i,j [q t ](·) is a multilinear mapping, calculated here at the vector h t , and B i,j [q t ] is a function; both depend on q i t , i ∈ I N , nonlinearly (in general).

Conjecture. Consider another auxiliary functions on
recall that q iO ∆t depends on t. Then the vector h t = (h i,j t,∆t ) i,j∈I N satisfies a homogeneous system of linear differential equations: The system of linear equations (2.39), can be solved in matrix form (or, rather, tensor form, as vector h t,∆t is two-dimensional). Note that the initial condition to (2.39), when ∆t = 0, can be obtained, by (2.38), (2.33), (2.35), (2.32), as follows: Next, if i = j, one has to find q iO ∆t (that can be often done explicitly), and get g i,j t,∆t (ξ) from (2.38). Finally, one has to take the inverse Fourier transform, to obtain g i,j t,∆t (x); the latter, of course, can be done only numerically. As a result, one gets an approximate value of u i,j ε,t,∆t (x) from (2.31).
3 Case study 1: Spatial and stochastic logistic model

Spatial characteristics
We consider agents of one type, i.e. N = 1. Let L be given through some of three operators, cf. 1.4: Here m > 0 is a constant, and a ± (x) ≥ 0 are kernels such that We denote also We will always assume that A ± > 0, i.e. it is not the case that a ± (x) = 0 for almost all (a.a. henceforth) x ∈ R d . Operator L 1 describes that any catalyst at x ∈ γ may create a product at y ∈ R d (send an off-spring to y) according to the dispersion kernel a + ; L 2 describes that any reactant at x ∈ γ may disappear with an density independent mortality m; L 3 describes that any reactant at x ∈ γ may also disappear because of competition with catalysts at y ∈ γ \ {x} given through the competition kernel a − .
Following the scheme above, we consider L ε with a ± (x − y) above replaced by ε d a ± (εx − εy); next, we consider the dynamics of γ ε,t defined by (1.5) ε,t (x 1 ) and k (2) ε,t (x 1 , x 2 ) be the corresponding first-and second-order correlation functions, and let u be the corresponding first-and second-order cumulants. Consider the space homogeneous case. Then, by e.g. [5], Function q t , cf. (1.22), satisfies the mean-field equation which can be solved explicitly: ( 3.4) If we assume, additionally, that we obtain, see [3], an equation for the Fourier transform of g t (x), namely Equation (3.3) has two stationary solutions q t ≡ 0 and q t ≡ q * , where We will always assume that i.e. that q * > 0; otherwise, by (3.4), lim t→∞ q t = 0, i.e. the population would extinct.
Under (SL 1 ), we have lim By (3.1), J t is integrable; as a result, (3.5) holds with g 0 replaced by g t . In particular, cf. (1.29), According to (2.36), we define also By (3.6) and (3.3), is just a multiplication operator given by, cf. (3.7), for a function f : R d → R.
We will consider below the stationary regime, when t → ∞. We define, for a.a.
We will assume, additionally to (SL 1 ), that there exists α > 0, such that Since, for an integrable function f ≥ 0, one has, by (3.15) and (SL 1 ), the following sufficient condition for (SL 2 ): It was shown in [3], that if (SL 1 ), (SL 2 ), (3.10) and (3.5) hold, then there exists Surely g * (ξ) is just the stationary solution to (3.6), i.e. it satisfies (3.6) with the left hand side replaced by 0. It was also shown in [3] that the inverse Fourier transform g * (x) of g * (ξ) is just the pointwise (and even uniform) limit of g t (x) as t → ∞.

Spatiotemporal characteristics
We consider the auxiliary dynamics of γ t = (γ O t , γ + t , γ − t ) described by the generator where L 1 , L 2 , L 3 are defined according to the rules postulated in Subsection 2.2.
i.e. both O-reactants and +-reactants may produce +-products. Next, in the counterparts of L 2 and L 3 , O-reactants become −-products, whereas +-reactants just disappear. Therefore, and since there are both O-and +-catalysts, we have Following the general scheme, we consider L ε with a ± (x − y) above replaced by ε d a ± (εx − εy); next, we consider the dynamics of γ ε,t defined by (2.18) ∆t , A, B ∈ {O, +, −} be the corresponding functions from the beyond mean-field expansion.
One can now solve, for each ξ ∈ R d , a linear ordinary differential equation (3.24) with the initial condition where h t (ξ) is given by (3.12). Then, one can get g t,∆t (ξ) from (3.21); to this end, one needs q O ∆t . The latter function satisfies the following differential equation (again, see Subsection 3.3 for details): where we used (3.23). As a result, we will get the following statement.
where q t and q t+∆t can be obtained from (3.4). If, additionally, (3.10) and (SL 2 ) hold, then g t,∆t (ξ) is an integrable function, and one can apply the inverse Fourier transform to it, to get g t,∆t (x) for a.a. x. Then, for all t ≥ 0 and a.a. x ∈ R d ,

28)
and g ∞,∆t (ξ) is an integrable function. Let g ∞,∆t (x) be its inverse Fourier transform. If, additionally, g 0 (x) and g 0 (ξ) are both integrable, then, for all ∆t ≥ 0, the following limit holds uniformly in a.a. x ∈ R d :

Derivation of equations
In this Subsection, we are going to derive equations (3.22) and (3.25). We will partially use the Model Constructor toolbox presented in [1]. Firstly, we express L 1 , L 2 , L 3 given above as sums of model components in the terminology of [1,Supplementary Note 2]: Here L 11 represents the Birth component: L 12 represents the BirthToAnotherType component: L 21 represents the DensityIndependentDeath component: L 22 represents the ChangeInType component: L 31 represents DeathByCompetition: L 32 represents DeathByExternalFactor: L 33 and L 34 require a new component called ChangeInTypeByFacilitation (which is defined below): The Model Constructor is written on Wolfram Language and requires Wolfram Mathematica R v10 or later. The Model Constructor packages are available at [2] and should be installed before running the following code. Next, we define the ChangeInTypeByFacilitation model component needed for L 33 and L 34 above. It describes the event when an agent at a position x 1 changes own type from s 2 to s 1 . The event happened because of interaction of the agent with each of other agents of a type s 3 placed at a position x 2 . The interaction is defined through a kernel a(x 1 − x 2 ). In particular, s 3 may be equal to s 1 as it is needed for L 33 .
In [5] We define now the AuxiliaryProcess which includes all model components corresponding to operators L ij above. Here the agent types 1, 2, 3 correspond to O, +, −, respectively. We sum up now the right hand sides of the differential equations for the needed g A,B ∆t (ξ), A, B ∈ {O, +, −}, and q O ∆t , cf. (3.20), (3.21). In the notations of the Model Constructor: We use the following code: We are going to verify now (3.22); to this end, we define the expression for h := h t,∆t (ξ), cf. (3.20), (3.21): In[19]:= (*Define h*) h=g [1,1,ξ]+g [1,2,ξ]+g [1,3,ξ]+g [2,3,ξ]+q [1]; Finally, we equate Ch with the obtained sum of the right hand sides of the equation, and find C, that is nothing but 1 2  Here q[2] = q + ∆t and also, by the very definition (1.29) of the Fourier transform: Therefore, the found expression for C coincides with the factor before h t,∆t (ξ) in the right hand side of (3.22). The second found alternative just means that h = 0, i.e. that h t,∆t (ξ) ≡ 0 also solves (3.22), that is trivial. Hence, (3.22

Numerics for the stationary regime on plane
We consider the 2-dimensional case: d = 2, and radially symmetric kernels with equal Gaussian shapes: Note that, indeed, R 2 β(|x|) dx = 1, and also, see Example 3.1, the assumption (3.17) holds, hence (SL 2 ) holds. We will need the following lemma.
where β is given by (3.30). Then Moreover, if f is a function such that the function g(ξ) := f ( c(ξ)), ξ ∈ R 2 is integrable, then the inverse Fourier transform g(x)of g(ξ) can be found by the formula

32)
where J 0 is the Bessel function of the first kind.
The latter integral can be calculated numerically. We use the following Wolfram Mathematica code (where dt = ∆t and r = |x|): The simulations were done with ε = 1 2 . The covariance between numbers of agents in two areas satisfies (2.25); note that we are actually interested in the covariance between 'small' areas (a local characteristic), so we may assume that they are disjoint, hence, the second summand in the right hand side of (2.25) is redundant. Next, the value of u ε,t,∆t (x) can be approximated by the formula (2.31), with d = 2 and ε = 1 2 , in our case. Therefore, we are interested in We plot now graphs for 1 4 g ∆t, r 2 with ∆t ∈ {0, 1, 2}, r ∈ [0, 10]:

Spatial characteristics
We consider now agents of two types, called hosts and parasites: Let L be given through some of three operators, cf. 1.4: Here L 1 describes an independent birth process of hosts: any H-catalyst sends an off-spring which is an H-product, according to a dispersion kernel a + ≥ 0: next, hosts may die because of the competition with other hosts (for resources), according to a competition kernel a − ≥ 0: next, parasites may die with a constant mortality rate m > 0: finally, any host may be transformed to a parasite (keeping the position) because of interaction with the existing parasites, according to a kernel b ≥ 0: We will assume that (3.1) holds for both a ± and for b, we define A ± through (3.2), and set, similarly, We consider L ε by replacing a ± (x − y) and b(x − y) by ε d a ± (εx − εy) and ε d b(εx − εy), respectively. We consider the space homogeneous case. Then, by the general scheme described in Section 1, (4.1) We define also, cf. (2.36), Differential equations for q A t and h AB t (ξ) = g BA t (ξ), A, B ∈ {H, P } are derived in Subsection 4.3 below. We will show that Next, we define and consider the matrix The second and third rows will correspond to h HP t (ξ) and h P H t (ξ) which are equal. Hence, we consider also the matrix with swapped second and third rows (4.6) Finally, we define the vector-function henceforth, the superscript T denotes the transpose vector. Note that Then, we will show in Subsection 4.3 that, for the vector we have, cf. (2.37), We consider now the stationary regime when t → ∞. The only pair of non-zero stationary solutions of (4.3) is Therefore, the condition ensures that (q H * , q P * ) is a the only pair of positive stationary solutions to (4.3). Proposition 4.1. Let (HP 1 ) hold. Then, for any q H 0 > 0, q P 0 > 0, then q H t and q P t oscillate around q H * and q P * , respectively, with a decreasing amplitude (damping oscillation).
We define also the following analogue of (3.14): for x, ξ ∈ R d , we set (4.14) We will assume that there exist α > 0, such that The following proposition provides simple sufficient conditions for (HP 1 )-(HP 3 ).  3. An example when (4.15) holds is the case (3.19) of an equal shape c of kernels a + = A + c and a − = A − c, provided that (HP 1 ) holds. Indeed, then

Proposition 4.2. Suppose that
In the Appendix below, we consider also how the condition (4.16) can be relaxed to still get (HP 3 ).

Spatiotemporal characteristics
Similarly to Subsection 3.2, we consider the auxiliary dynamics of described by the generator where L 1 , L 2 , L 3 , L 4 are defined according to Subsection 2.2. Namely, now we have, similarly to the corresponding operators in spatial logistic model The situation with L 4 is more complicated: when HO-reactant of the spatial dynamics becomes P +-product (keeping the position), it should be also transformed to H−-product, according to the general scheme.
Hence, formally, one could write: However, the first summand in the latter expression does not satisfy the basic requirement that different agents should have different positions: there are agents x H− and x P + simultaneously. 1 To overcome this, we consider a formal modification the host x HO transforms to the parasite z P + distributed in space according to the kernel given by the Dirac δ-function δ(x − z).
Therefore, we define Considering the first summand of the operator L 4 as such with a regular (integrable) kernel δ(x − z) one apply the technique described above and derive the corresponding differential equations for correlation functions, cumulants, and for the beyond mean-field expansion. We should then 'replace' the regular kernel by the real δ-function. This latter includes two steps. Firstly, in differential equations in terms of the Fourier transform, we replace δ(ξ) by 1.
Secondly, one has to distinguish the terms which, in real coordinates (before passing to the Fourier transform), contained the δ-function, i.e.
If such term contained also the point x H− , then, after the integration w.r.t. z P + = z, the corresponding integral would disappear, and z will be replaced by x. The latter is however the position of the H−-agent. We will get as a result the δ-function between the positions of H− and P + agents. This can be written heuristically as follows: for a (regular) function f , It means that the pair correlation between H− and P + will include now a δfunction. By (2.17) (considered in the case i, j = I 2 = {H, P }), pair cumulant u H−,P + ∆t = u P +,P − ∆t in the r.h.s. appears for u HP t,∆t in the l.h.s. only. Hence, in the space homogeneous case, u HP t,∆t (x) is (the only) non-integrable pair cumulant. This is the effect of using L 4 instead of L 4 . Since the dynamics we consider is linear, one has u HP t,∆t (x) = u HP t,∆t (x) + u HP t,∆t,sing δ(x), (4.23) so, to get the answer, that is the corresponding cumulant u HP t,∆t (x) for the model generated by L 4 , one has to get rid of the singular term. We will do this in the beyond mean-field approach below. 2 Following the general scheme, we consider L ε with a ± (x − y) and b(x − y) above replaced by ε d a ± (εx − εy) and ε d b(εx − εy), respectively; next, we consider the dynamics of γ ε,t defined by (2.18) with L replaced by L ε . Let q X ∆t , g X,Y ∆t , X, Y ∈ {HO, H+, H−, P O, P +, P −} be the corresponding functions from the beyond mean-field expansion. We consider, cf. (2.32), for all A, B ∈ {H, P }, (4.24) and, cf. (2.38), (4.25) Recall that, g HP t,∆t = g P H t,∆t and hence h HP t,∆t = h P H t,∆t . By (4.23), g HP t,∆t (x) is (the only) non-regular function which should include the Dirac δ function, as a result g HP t,∆t will be non-integrable at infinity. By (4.23), if we subtract from g HP t,∆t its limit at infinity, the result will be nothing but the Fourier transform of g HP t,∆t (x), so that, cf. (2.31), AB ∈ {HH, HP , P H, P P }. (4.26) We will get hence the values of the spatiotemporal cumulants for the initial model. In Subsection 4.3 below, we will show that d d∆t (4.28) where D t,∆t (ξ) is the following matrix: and hence, cf. (4.4), . (4.29) Therefore, by (4.5), we can represent A t+∆t (ξ) as a block matrix, namely, (4.30) where 0 denotes 2 × 2 matrix of zeros, and hence, denoting, cf. (4.9), we get, by (4.27), (4.28), (4.31) cf. also (4.10). Note that, by (4.29), cf. (4.8), A t+∆t (ξ) = A[q H t+∆t , q P t+∆t ](ξ). Comparing this with (4.10), we see that the conjecture holds.
An explicit formula for g HP t,∆t (∞) is provided in the Appendix below, see (A.33). As a result, for each AB ∈ {HH, HP , P H, P P }, one can apply the inverse Fourier transform to an integrable function g AB t,∆t (ξ), to get, for a.a. x ∈ R d , the needed for (4.26) function g AB t,∆t (x).
As a result, we obtained the desired leading terms g AB ∞,∆t (x) in the beyond meanfield expansion (4.26), for each AB ∈ {HH, HP , P H, P P }, of the spatiotemporal cumulant of the considered model in the stationary regime, i.e. when t in (4.26) is replaced by ∞.

Derivation of equations
In this Subsection, we are going to derive, in particular, equations (4.3), (4.10), and (4.31). Similarly to Subsection 4.3, we will partially use the Model Constructor toolbox from [1]. Firstly, we note that the operators L 1 , L 2 , L 3 , L 4 defined above represent Birth, DeathByCompetition, DensityIndependentDeath and Infection components, respectively.
Next, we express L 1 , L 2 , L 3 , L 4 given above as sums of model components in the terminology of [1,Supplementary Note 2]. Namely Here L 11 represents the Birth component: L 12 represents the BirthToAnotherType component: L 21 represents the DeathByCompetition component: L 22 represents the DeathByExternalFactor component: L 23 and L 24 both represent the ChangeInTypeByFacilitation component introduced in Subsection 4.3: L 31 represents the DensityIndependentDeath component: L 32 represents the ChangeInType component: L 41 also represents the ChangeInTypeByFacilitation component: L 42 represents the Infection component (which is, actually, a partial case of the ChangeInTypeByFacilitation component): Finally, L 43 and L 44 both represent the ChangeInTypeAndBirthByFacilitation component defined below: We are going to describe the Wolfram Mathematica code we used. First six lines (In [1][2][3][4][5][6]) are the same as in Subsection 3.3: we load libraries, set-up internal variables, define the ChangeInTypeByFacilitation and Relax components.
Next, we define the ChangeInTypeAndBirthByFacilitation model component needed for L 43 and L 44 above. It describes the event when an agent at a position x 2 changes own type from s 3 to s 2 and, simultaneously, it sends an off-spring of a type s 1 to a position x 1 . The off-spring is sent through a kernel d(x 1 − x 2 ) (which is the Dirac δ-function δ(x 1 − x 2 ) in operators L 43 and L 44 ). The event happened because of interaction of the agent with each of other agents of a type s 4 placed at a position x 3 . The interaction is given through a kernel b(x 2 − x 3 ). In particular, s 4 may be equal to s 1 as it is needed for L 44 .

In[7]:=
To simplify the representation of the calculations below, we introduce a replacement rule to replace values of the Fourier transform at the origin by the corresponding integral, i.e.
we include here also the replacement δ ≡ 1: In the Model Constructor toolbox, HQfALL and also HGfALL are the functions providing the r.h.s. of the differential equations for functions q and g, respectively.
Firstly, we create a collection rule: ,-q [4],-q [5], q [1],q [2],q [3],q [4],q [5]}; Next, we collect the terms in the r.h.s. of the equations obtained in In[14] according to the rule above; and then one can replace the sums of g by the corresponding notations for h AB t,∆t . Then we expand, and thereafter we collect again, now to distinguish the coefficients before functions h:

True
We are going now to define the functions in the stationary regime, using the replacement according to (4.12). Firstly, we consider the stationary version (the limit as t → ∞) of the matrix (4.5): Finally, we are going to verify the corrections obtained in (4.34)-(4.35) to get integrable functions. We were actually interested at the values of the obtained solutions h AB ∞,∆t (ξ) as |ξ| → ∞. By the Riemann-Lebesgue lemma, the Fourier transforms of a ± , b converge to zero at infinity.

Analysis and numerics for the stationary regime
By Theorem 4.7, the eigenvalues of the matrix E * (ξ) have negative real parts for each ξ ∈ R d . Then, by (4.36)-(4.37), functions h AB ∞,t,∆t (ξ), A, B ∈ {H, P } converge (for each ξ) to 0 as ∆t → ∞. Similarly to Proposition 4.1, this convergence to 0 can be monotone or oscillating with damping, depending on whether the eigenvalues of E * (ξ) are real negative or complex with negative real parts, respectively. By (4.34) and (4.35), functions g AB ∞,∆t (ξ), AB ∈ {HH, HP , P H, P P } have the same properties.
By the proof of Theorem 4.7, for a fixed ξ ∈ R d , the eigenvalues of E * (ξ) are real negative iff is non-negative; otherwise, they are complex with negative real parts. The character of convergence hence depends on an interplay between behavior of the Fourier transforms a + (ξ), a − (ξ) and b(ξ) in different zones of ξ ∈ R d , and is especially non-trivial if the Fourier transforms take negative values, that may be the case e.g. for the Gaussian-like kernels considered in (A.18) in the Appendix.
One can make, however, several general observations about the mentioned convergence.

Since
we conclude that the convergence is oscillating for ξ = 0 if and only if (4.13) holds. By continuity and the strict inequality in (4.13), the oscillations will take place for ξ at a neighbourhood of the origin.
2. Directly from (4.39), one gets that the convergence is monotone for all ξ ∈ R d such that b(ξ) ≤ 0. We assume now that this is not the case and (4.40) 3. To have oscillating convergence for all large enough |ξ|, one needs with necessity, by the Riemann-Lebesgue lemma, that As a result, to have oscillating in neighbourhoods of the origin and infinity, we require, with necessity that both (4.40), (4.41), and (4.13) hold. Note also that (4.13) under assumption (4.41) can be easily rewritten as follow: We assume now that both (4.40)-(4.42) hold; moreover, we also consider the case of equal shapes for all kernels: (4.43) Then, since, by (3.16), c(ξ) ≤ 1, one can easily rewrite (4.39) as follows: because of (4.42). Therefore, (4.40)-(4.42) imply that that all functions g AB ∞,∆t (ξ), AB ∈ {HH, HP , P H, P P }, converge to 0 as ∆t → ∞ with damping oscillations for all ξ ∈ R d .
Note that the opposite to (4.42) inequality does not imply that the convergence is monotone for all ξ ∈ R d . Moreover, one can not show analytically that the oscillations will be preserved for the inverse Fourier transforms g AB ∞,∆t (x), AB ∈ {HH, HP , P H, P P }.
The result of function gReal is the vector of g AB ∞,∆t (x), AB ∈ {HH, HP , P H, P P }; here valt = ∆t and valx = |x|. The result will be used to find the corresponding cumulants by (4.26): the simulations were done for ε = 1 2 , hence, the cumulants measured in simulations should be approximated as follows: AB ∈ {HH, HP , P H, P P }.
We are going to discuss now the results of the numerical calculations of the inverse Fourier transforms.
1. For a fixed space variable, e.g. for x = 0, the graphs of the g-functions as functions of ∆t ∈ [0, 20] are shown on Figure 2a. They demonstrate the damping oscillations. As we can see, values of g HH ∞,∆t (0) and g P P ∞,∆t (0) are pretty close, as well as values of g HP ∞,∆t,reg (0) and −g P H ∞,∆t (0). We observe also the initial negative correlation between finding a host at the current position of a parasite for some positive time interval.
2. Next, in the considered case of monotone kernels, cumulants converge to 0 in |x| monotonically, see Figure 2b for a fixed ∆t = 7. In other words, for each |x| ≥ r 0 (perhaps, staring with some r 0 > 0), we will see a picture similar to that on Figure 2a, however, the corresponding cumulants have smaller amplitudes.
3. We consider now the dependence of solutions on A + when A + = m. We keep other parameters in (4.48) and compare the graphs of Figure 2a with those for A + = 2 > m and A + = 1 2 < m. We present the comparison on Figure 3, for e.g. |x| = 5. We can observe that, with the growth of A + , the oscillation became more frequent (the period becomes smaller) and the amplitude becomes higher.
4. Finally, recall that the parameters in (4.48) satisfy (4.42), that is, for A + = m, nothing but (4.13). We consider the case when (4.42) fails. If we keep A + = B = m = 1, then (4.42) fails iff A − ≥ 2( √ 2 − 1) ≈ 0.83. We take A − = 0.85, and then, see Figure 4, the convergence of g-functions in ∆t to 0 becomes monotone, without visible damping oscillations. Note that, however, we can not prove this analytically, so the fluctuation may still happen for large values of ∆t and/or |x|. Appendix: Mathematical proofs and discussions Proof of Theorem 3.2 Step 1. By (2.40), h t,0 (ξ) = g t (ξ) + q t .
Then, the solution to (3.24) has the form By the flow property of solutions to (3.3), q t+τ = q τ , where We have then d dτ and integrating from τ = 0 to τ = ∆t, we get Replacing back q τ by q t+τ , we obtain Substituting into (A.1) and using that m + q * A − = A + , one gets Then, one can rewrite: To find g t,∆t (ξ) from (3.21), one has to solve (3.25). By (2.35), q O ∆t = q t . Therefore, by (A.2), As a result, (3.21) implies (3.26).
Next, we have and by (A.9), (A.5), Finally, rewriting back and using that we obtain from (A.8), that Hence, by the dominated convergence theorem, (A.11) holds with I 1 replaced by each of I 2 , I 3 , I 4 . Therefore, g t,∆t converges to g ∞,∆t in L 1 (R d ) as t → ∞. Since g t,∆t − g ∞,∆t is the inverse Fourier transform of g t,∆t − g ∞,∆t , we obtain from an analogue of (3.16) that g t,∆t converges to g ∞,∆t in L ∞ (R d ).

Proof of Lemma 3.3
Firstly, for the Fourier transform defined as in (1.29), formula (3.31) follows from e.g. [4,Example 2.2.9,Proposition 2.2.11]. Next, we note that the Fourier transform of a (radially) symmetric function is also (radially) symmetric, and it coincides with the inverse Fourier transform of that function. Therefore, by e.g. [4,Back Matters B.5], Making the substitution s = 2πr, one gets the desired formula (3.32).

Proof of Proposition 4.1
We set f (x, y) = x(A + − A − x − By), g(x, y) = y(Bx − m). Consider the Jacobian Then, it is straightforward to check that Since under (HP 1 ), we conclude that both eigenvalues of j(q H * , q P * ) have negative real parts, hence (4.12) holds.
The monotone convergence in (4.12) takes place when the eigenvalues of J are real (and hence negative), otherwise there will be damping oscillations. The eigenvalues of J are real iff tr(J) 2 ≥ 4 det(J), that is equivalent to 2 , and hence the statement is proved.
with small enough s > 0. Indeed, it is straightforward to check that b(ξ) < 0 for |ξ| > d(s) for certain continuous d(s) with d(0) = 0. Therefore, taking s > 0 small enough we ensure that b(ξ) < 0 for some ξ ∈ Λ given by (A.17). Note also that a small s here corresponds to a large length scale of the kernel b(x) describing the influence of parasites on hosts.
operators on X) and also provided that the operator (matrix) C * is invertible. We have that C * (ξ) is given by (4.18) and Th show that C * is invertible in X it is evidently enough to show that the function det C * (ξ) is separated from 0. We have By (HP 2 ) and (A.13), the proves the needed. As a result, by [3,Lemma 3.1], and moreover g AB t → g AB * in both L 1 (R d ) and L ∞ (X) as t → ∞, A, B ∈ {H, P }. The convergence in L 1 implies, by (SL 2 ), the convergence in L ∞ (R d ) for the inverse Fourier transforms, that fulfills the proof.