A Many-Body Field Theory Approach to Stochastic Models in Population Biology

Background Many models used in theoretical ecology, or mathematical epidemiology are stochastic, and may also be spatially-explicit. Techniques from quantum field theory have been used before in reaction-diffusion systems, principally to investigate their critical behavior. Here we argue that they make many calculations easier and are a possible starting point for new approximations. Methodology We review the many-body field formalism for Markov processes and illustrate how to apply it to a ‘Brownian bug’ population model, and to an epidemic model. We show how the master equation and the moment hierarchy can both be written in particularly compact forms. The introduction of functional methods allows the systematic computation of the effective action, which gives the dynamics of mean quantities. We obtain the 1-loop approximation to the effective action for general (space-) translation invariant systems, and thus approximations to the non-equilibrium dynamics of the mean fields. Conclusions The master equations for spatial stochastic systems normally take a neater form in the many-body field formalism. One can write down the dynamics for generating functional of physically-relevant moments, equivalent to the whole moment hierarchy. The 1-loop dynamics of the mean fields are the same as those of a particular moment-closure.

S1 Sketch-proof of moment-hierarchy dynamical law.
In the main text, we mentioned that a generating functional Here we sketch a proof of this. First note that because the reference state is coherent, we have Also, it is not difficult to see that commuting an a † through a's to its left acts like differentiation. In multi-index notation therefore: where the derivative acts to the left. This provides a neat way of calculating the dynamics of a specific normal-ordered moment. In the formal limit, one can replace the partial derivatives with functional derivatives.
Consider further the behaviour of a piece of the hamiltonian which can be represented ← − ∂ p a q . Then X n ≥0 1 n! J n a n← − ∂ p a q = X n ≥p 1 n! J n n! n − p! a n−p a q = J p X n ≥0 1 n! J n a q+n = J p " ∂ ∂J « q X n ≥0 1 n! J n a n which together with the above, proves the result.

S2 Coherent-state path integrals and the effective action
The material in this section is standard to field theory. It is included here to set conventions and for the benefit of readers who do not have a physics background. For more details, see e.g. Refs. [6,11] of the main text.
Write φ for a field where φα (vaguely discretised) refers to site α. Then take our coherent state definition as These states have the properties

S2.1 Coherent state path integral
For illustration, we consider where the M δt k give equally spaced intervals which partition [0, t]. Before each k−th element in this product, we insert (in slightly abusive notation -dropping subscripts to the αs where we can get away with it). Then One then uses e δtH(a † ,a)=:e δtH(a † ,a) :+O(δt 2 ) (10) and the properties of coherent states above to obtain (ignoring O(δt 2 )) Taking various 'limits', we write this formally as and similarly for other 'observables'.

S2.2 Path integral representations of interesting quantities
We willconsider initial states which have poisson-distributed numbers of bugs uniformly in space. This gives a multiple of a coherent state: The final state will always be |1 for us, which is coherent. So for this state, e.g.
etc. More general initial probability distributions can be considered by expanding them in terms of coherent states.

S2.3 The shift trick
In [7], Cardy and Täuber point out one can make a shift φ * = 1+φ, which changes the term in the exponential to Then we have such things as For the bug case then . We will write this more concisely as:

S2.4 Relation of coherent states to probability distributions
A coherent state corresponds to poisson sprinkling in each infinitesimal area. In the translation invariant case, this means a poisson choice of total number sprinkled uniformly across the area. For low densities, this state is a good approximation to a fixed initial number of bugs.

S2.5 Calculation of Γ
Introduce a counting parameter , which will in fact count the number of loops by its power. From the definitions we have We have substitued for J from its definition as the derivative Γ and shifted the variable of integration. Integer subscripts indicate functional derivatives, Q is the part of the action quadratic in the shifted variable φ, and R is the functional taylor series for the shifted action from third order up. Writing ∆ = Γ − S, this reads where we have made the rescaling φ/ 1/2 → φ. (Rescalings of the 'measure' do not concern us because they respresent an irrelevant constant shift in Γ. More careful treatments consider only the ratio of such integrals as meaningful: see, e.g. [11].) This yields an equations which we can solve iteratively for Γ, order by order in : Although it appears that this series contains non-integer powers of , this is not the case because only even moments of a gaussian integral are not zero.
The lowest term in this series reads: where capitals have been used for functional determinants and traces and we have dropped irrelevant constants.

S3 The general multimode fluctuation integral
This section extends the mutlimode fluctuation integral to non-hermitian hamiltonians with sources; a case which does not appear to have been treated previously. It is this case which is needed for effective action calculations for non-quantum situations, such as those we treat. We closely follow the techniques of [14].

S3.1 The action
We are interested in the coherent state path integral integral where the integral has zero limits and The discretized version of the action is (indices refering to time-slices) with appropriate . Writing a = x + iy,ā = x − iy and using where convenient that ξ0 = ξn+1 = 0, the action reads Here, following [14], we have introduced the coordinate system where x and y live in R k for a k-mode system. Using this, the action can be written compactly as This allows the identification

S3.2 The integrations
Repeated use of the formula allows us to see that In, the discretized version of I, is given by where the M j are defined recursively by

S3.3 Solving the recursions
One shows by induction that The following properties of the gamma's will be used: and Dropping the subscripts for the moment: Note that in the recursion, the L's carry factors of Γ3 and Γ T 3 , only Γ1 and Γ2 terms survive this, and both yield Γ2. This shows the form given above is correct. We wish to compute the recursion satisfied by X to order however. It is not hard to see that up to order we have Hitting these with the L's one finds that: The continuum version of this is precisely It only remains to calculate Note both Γ1 and Γ2 are traceless. Using det(1 + A) = 1 + trA + ((trA) 2 − trA 2 ) + . . ., and keeping only order , the above becomes We then have

S3.5 The u part
From the above relations, it is clear that the form u k = Y k ⊗ µ1 + Z k ⊗ µ2 holds. To order then: which means which becomes in the continuous limit: Up to order we have Which means that the u contribution becomes

S3.6 Summary
We have shown that (up to an irrelevant constant factor): where, with trivial initial conditions In the main text, we have written Y rather than Z to avoid confusion with generating functionals.

S4 Explicit computations for the two models considered
In this section, there are more details of the calculations pertaining to the specific models we consider.

S4.1 Non-spatial bugs
For our action f , g and ω are given by Varying this yields the zero loop equations plus a correction (note that φ ≡ 0 after variation): Note that the dynamical equations inφ are satisfied by the trivial solution. Further, ∂φf (ψs) = −V giving us the system For comparison, the central moment-closed equations read with the first and second moment (n and N2) having initial conditions N0 and N 2 0 + N0 respectively. N3 is the third moment, posited to satisfy the closure relation in Eqn.56.

S4.2 Non-spatial SIR
From the SIR-action above, one deduces that where the field letters now refer to the mean fields. The 1-loop effective action is: After variation, it is clear that the the barred fields are non-dynamical, and stay zero (as they should). Given this, and writing The differential equations for X become: coupled to the equations of motion: Note that only the first three of the equations in A to F are relevant to the correction term 2B, in the same was that only 3 of the second moments are relevant in moment-closure. For comparison, the zero third central-moment equations are with the covariance CSI and variances VI and VS satisfying the initial conditions 0, I0 and S0 respectively. The covariance plays the role of the correction to the rate equations for S, I and R.

S4.3 Spatial bugs
We want to write the action in terms of its fourier modes. For a field ψ write where d is the number of spatial dimensions. We have where we have used Also, φn 1 Vn 2 φn 3 e 2πi/L(x.n 1 +n 2 .(x−y)+n 3 .y) = X n 1 ,n 2 ,n 3 The quadratic part of the bugs action is (in condensed notation) Where we have introduced Assume X−n = Xn, equivalent to X ∈ R. Summarizing, this gives where we have dropped the m-subscript. The multimode result is unchanged by the L d prefactor, and each mode is in fact uncoupled from the rest, giving the effective action 1-loop term in the form Similar arguments apply to the one loop case, which means we need only consideṙ

S4.3.1 Converting back to real space
Starting witḣ Xn (74) one can write the correction in terms of real-space again. First note that Vn 1 Xn 2 e −2πix.n/L e 2πiy.n 2 /L e 2πi(x−y).n 1 /L = L d VnXn We therefore have The correction to the tree-level equations of motion for φ is then This compares with a moment-closure correxion of the form where for 3rd-cumulant-zero closure These are the same, with the identification C = 2X/L d .

S4.4 Spatial SIR
The quadratic part of the action is given by where In fourier space (81) As above, we will convert all the −n in this equation to simply n.
Again, everything separates out by mode. The one loop correction is of the form 2L d P n trfnXn witḣ Dropping the ns and writing  (87) where for these calculations The correction antisymmetrically applied to theṠ andİ equation is then 2L d P n βVnBn.

S4.4.1 Converting back to real space
The real-space equations are We expand the integration variables a andā around paths a0 andā0, with the correct initial and final conditions, respectively, to give δā,δa] (93) where SQ[δā, δa] is the quadratic part of the action and Sint[δā, δa] is everything else. The idea is to insert a counting parameter in front of Sint and expand perturbatively in this: One truncates this at a given order, e.g. 2: and then asserts a principal of minimal sensitivity (PMS). Namely, that

S5.2 The mean fields
This transition probability is a special case of a generating functional with J = K = 0. We can carry out the above procedure for more generally. Then we can calculate the mean fields as for a truncation at order n.

S5.3 Relating to number states
As defined above, we are referring to the transition probability between Poisson distributions: For large n these things are not so different, but one can also use the completeness relation to write CHECK!

S5.4 The non-spatial SIR model
The original action is which means that which means that The extra terms in . . . Q are, at order (using the fact that the expectation of an odd number of fields is zero) and at order 2 The integrand is +ā0 (1)  In our case

S5.5 Moments
Calculation of the functions L(t) and q(t, t ) requires the computation of the moments of . . . Q up to order 8. These can be calculated from the generating functional I, as above: Henceforth, we will confine ourselves to computation of mean fields as these are simpler. The equations admit a solution withā ≡ 0. This can be seen as follows. The δS/δa terms are all at least linear in barred fields, and so the only forcing is from the L and q terms. If the barred fields are zero, then f ≡ 0, so that The first of these implies that and hence that higher derivatives of Y with respect to K are zero This further implies that δ i+j I δJ i δK j˛0 = 0 (i < j) and this leaves no forcing on the barred dynamics, allowing the solution a0 ≡ 0. The terms of the integrand Eqn. 110 which are are relevant are: +b0 (1) which agrees with Eqns. 120 and 121 with the understanding that θ(0) = 0. One could obtain this result directly, if more delicately, by considering the discretization which underlies the differential equations.