A unified framework for analysis of individual-based models in ecology and beyond

Individual-based models, ‘IBMs’, describe naturally the dynamics of interacting organisms or social or financial agents. They are considered too complex for mathematical analysis, but computer simulations of them cannot give the general insights required. Here, we resolve this problem with a general mathematical framework for IBMs containing interactions of an unlimited level of complexity, and derive equations that reliably approximate the effects of space and stochasticity. We provide software, specified in an accessible and intuitive graphical way, so any researcher can obtain analytical and simulation results for any particular IBM without algebraic manipulation. We illustrate the framework with examples from movement ecology, conservation biology, and evolutionary ecology. This framework will provide unprecedented insights into a hitherto intractable panoply of complex models across many scientific fields.

T here are many systems in biology 1 , the social sciences 2 and other disciplines comprising populations of individuals that move and interact with each other, but where the properties of focal interest arise at the whole-population level. In ecology and evolution 1,3,4 , for instance, large-scale patterns, such as coexistence, biodiversity, and the evolution of novel types emerge from individual-level processes, such as births, deaths and inter-agent interactions. Historically, however, many fundamental ideas in ecology were developed using differential equation models that treat populations as continuous numbers.
The mathematical transparency of these models means that we understand their predictions completely and can propose general principles, e.g., regarding the number of resources needed by competing species 5 or the processes that destabilise host-parasite dynamics 6 . However, it has long been recognised that these simple models are inadequate for describing many biological phenomena 7 . For example, the conditions for competitors to coexist, and for predator-prey cycles to occur, are altered when space and individual discreteness are considered 8 . Space and stochasticity can similarly play an important role in the dynamics of chemical species 9 . Individual-based (or agent-) based models 10 faithfully capture the discrete and spatial nature of population dynamics, but these are usually studied by computer simulation 1 , which only tells us about a limited set of parameter values and not the general model behaviour. They therefore do not yield the sort of general understanding that would be given by mathematical analysis, which encapsulates behaviour across all possible parameter combinations 8 .
To use individual-based models to develop general principles, of the sort derived from classic differential equation models, requires a method for analysing them mathematically. Individualbased models can be formulated by spatiotemporal point processes 11 , where individuals (or, more generally, entity types, such as juveniles or infected individuals, see Fig. 1a) are created, destroyed and move at rates that can depend on the positions of other individuals in the system, see Fig. 1b. In principle, the dynamics of such systems are described exactly by equations for the time evolution of spatial moments, representing mean density of individuals (first-order moment), spatial covariance (secondorder moment) and so on (Fig. 1c). However, in practice there are two obstacles to using these to predict model behaviour. First, they need to be laboriously derived separately for each model, and for each order of moment 12 . Second, the moment equations form an unclosed hierarchy, with the dynamics of each moment depending on the higher order moments, so if the equations for all orders of moment were known it would still not be possible to solve them-even numerically 13 . 'Moment closure' is a widely used approximation scheme that closes the hierarchy by an adhoc assumption relating moments of different order 12 , but this is an uncontrolled approximation that is not guaranteed to perform well in any particular limit 8,14 . A more reliable alternative based on a perturbation expansion has been proposed, which gives asymptotically exact results when agents interact over large enough scales, but the algebraic burden for this method remains prohibitive because it requires arduous derivations for each particular model 11,13 (Fig. 1d).
Here, we overcome these difficulties by formulating a unified theoretical framework for a wide class of systems, which allows us to derive general analytical results. The results for any particular model within this general class can be obtained directly without further analysis. We make analytical results available to the nonspecialist by providing software, which generates mathematical expressions describing a model that the user specifies in an accessible and intuitive way. We first describe our framework and the mathematical results leading to it, and then give three applications that illustrate how our method allows us to answer questions that are not addressable by a simulation approach alone.

Results
The framework. We begin by classifying the participants in demographic processes into three types of individuals (borrowing terminology from chemistry): (i) reactants (that are destroyed by the process); (ii) products (that are created by the process) and (iii) catalysts (that are unaffected by the process, but whose presence affects the rate at which it occurs). For example: a death event has a single reactant and no products or catalysts; dioecious birth has two catalysts (the mother and the father) and one product (the offspring) and movement can be represented by a reactant at the initial position and a product at the final position. This representation can describe processes with an arbitrarily high degree of complexity: the number of reactants, catalysts and products that can participate in any event is unlimited, as is the number of entity types within the system. We could, for instance, model a population where individuals can consume a food item, thereby increasing their energy level, provided they are in the shelter of a tree and that there is another individual helping to capture the prey. In this case, there would be two catalysts (tree and helper), two reactants (original individual and food item) and one product (an individual with a higher energy level than the original), and different entity types would be used to represent trees, food items and individuals of the focal species with different energy levels. This general classification allows us to derive our first main result: an exact expression for the moment equations to all orders for the general model containing processes with arbitrary sets of reactants, products and catalysts and interactions between them (see the Methods section Eq. (3)).
Then, using this exact expression for the general model, we apply a perturbation scheme 13,15 to derive a robust approximation for the dynamics of this general class of model. If the interactions between individuals take place over very large spatial scales, the demographic rates depend only on globally averaged population densities. In that case, the moment equations reduce to 'mean-field' equations-ordinary differential equations of the type used in classical population dynamics. If the interactions are more local, then stochastic and spatial fluctuations cause the model to deviate from the mean-field behaviour. We assume that spatial interactions between individuals depend on their separation x (in d-dimensional space) according to interaction 'kernels' of the form ϵ d f ðϵxÞ, where 1=ϵ is the typical length scale of interactions, and prove mathematically (Supplementary Note 1) that the mean densities and spatial covariance (including autocovariance) satisfy an expansion where q is the density computed by the mean-field model, p is the correction to the mean-field density due to spatial stochastic fluctuations and g is the dominant contribution to the spatial covariance that describes the degree to which individuals are aggregated or segregated in space (oðϵ d Þ denotes a term that, when divided by ϵ d , vanishes when ϵ ! 0). For simplicity, in the above expression we have assumed translational invariance so that the density does not depend on position and the spatial covariance depends only on spatial separation, but our underlying mathematical framework does not have this limitation, see Fig. 1e and Eq. (1) therein. Depending on the scientific question being asked, it may be enough to study the mean-field behaviour described by q alone, whereas p or g will need to be calculated if we are interested in the nature and consequences of spatial patterns that emerge from local interactions. This perturbation scheme has been applied to some specific ecological models before, and has been found to perform as well as or better than alternative moment closure approaches 13,15 . The framework described here allows us to apply this technique to a much wider class of models. Our second main result is an expression for differential equations for q, p and g for the general reactant-catalyst-product process ( Fig. 1e and Eq. (2) therein, and Eqs. (4), (5), (6) in the Methods section). This allows us to construct directly the perturbation expansion for any model containing any number of processes in our general reactant-catalyst-product class because, even though the interactions between individuals may be nonlinear, the differential equations for q, p and g contain sums of independent contributions from each process. Our results make mathematical analysis available for systems that would be prohibitively complex using other approaches, and, while Eqs. (4), (5), (6) may appear daunting, we have automated this process by writing Mathematica code that computes analytical expressions for general reactant-catalyst-product models (see Supplementary Note 2). To complement and verify the analytical results, we have also written computer code in C that simulates a very broad class of reactant-catalyst-product models (see Supplementary Note 2). This provides a unified framework for analysis and simulation (illustrated in Fig. 2), where the user wanting to study a particular model has only to convert the verbal model description into a graphical model a c d e b Fig. 1 Agent-based models in ecology and evolution: the problem, and our solution. Agent-based models consist of entities of different types (panel a) that can interact in space in many ways (panel b). At any time, the state of the model can be quantified by spatial moments, such as density or spatial correlations (panel c). Agent-based models are generally considered too complicated to analyse mathematically (panel d). We address this problem by defining a general class of models and processes, for which we provide computer code that yields analytical solutions and runs stochastic simulations (panel e) description, which can easily be translated to the syntax understood by the computer code. For simplicity, our software assumes a translational invariant initial condition, because there are many ways in which the underlying geometry could be spatially heterogeneous (see the Methods section for more discussion of this point). Below we provide a step-by-step description of the application of the framework to derive the underlying equations for a dynamic landscape patch occupancy metapopulation model. We then use three case studies drawn from different areas of ecology and evolution to illustrate how the framework of Fig. 2 allows researchers to obtain deep analytical insights into important problems.
Using the framework. The use of our framework follows the workflow illustrated in Fig. 2. We demonstrate this with an example for which the perturbation expansion has previously been derived via other methods: a dynamic landscape patch occupancy metapopulation 15 . For simplicity, we will here consider the case without landscape correlations (the limit ν ! 0, αν ¼ constant in ref. 15 ).
We begin by deriving the equations for the sub-model that describes the dynamics of the habitat patches themselves (without considering whether they are occupied), which is explicitly solvable. We assume that individual patches are created, at random, at rate per unit area r, and destroyed (at random) at rate μ. Within our framework, this can be represented by a spatiotemporal point process with a single-entity type (habitat patches, which we denote as species type '1'). The formal mathematical representation of this model is the generator ; where the generators L IM and L D have technical definitions in the Supplementary Note 1, section 1.3, and can be understood intuitively as representing an immigration process (patch creation) with intensity r, and a density-independent death process (patch destruction) with rate μ. The subscript '1' states the species to which the process applies (in this case, habitat patches are the only 'species'). These are two examples of the reactant-catalyst-product processes that, for convenience, we have defined in Supplementary Note 1 section 1.3, but note that we have derivations for the general reactant-catalyst-product process (Supplementary Note 1 section 1.4)-we could specify the model explicitly in terms of such general processes, but this would be somewhat more long-winded than using the above shorthand.  Fig. 2 Steps in our framework for constructing models and calculating their predictions. These are illustrated using our dynamic landscape metapopulation case study. The user wanting to study a particular model has only to convert the verbal model description (a) into a graphical model description (b), formulated in terms of a list of possible events and the rates at which they occur. An unambiguous mathematical statement of the model as a spatiotemporal point process is obtained from the operators L for these components 11 ; the superscripts denote the corresponding class of reactant-catalyst-product process, namely immigration (IM), infection (I), change of type (CT), density-independent death (D). The colour denotes the type of individual (blue: occupied patch; pink: unoccupied patch). A computer-readable statement of the components (c) can be used to generate code for simulating the model (d), or to derive the equations describing the perturbation expansion using computer symbolic algebra (e). The differential equations for the mean-field q, correction p and spatial covariance g in the perturbation expansion comprise a sum of terms contributed by each of the model components. These expressions can then be evaluated numerically (f), or analysed further by human or machine (g). The predictions from simulations and mathematical methods can then be compared (h) and used in many ways to gain insights into the focal problem (see Figs. 3,4 and 5) We will now show the steps for deriving the equations for this model using our  Supplementary Table 3 for the full list of pre-defined processes, and Supplementary Note 2 sections 2.1.3 and 2.3.4 for details on how to define specific reactant-catalyst-product processes). The first argument '1' is the species label, the second argument (r or mu) is the corresponding rate and the third process is a (redundant) overall prefactor for the rates, included for consistency with other processes.
The quantity H q (the right-hand side for the differential equation for the mean-field patch density) is given by the Mathematica command: HQfALL[{q, p, g},processes,1] >r -mu q [1] where {q,p,g} specifies the names of the variables Mathematica will use to denote the mean-field density, stochastic correction and correlation functions, processes is the vector defining the model given above, and '1' is the species label. To obtain the quantities H p , which requires an integration in Fourier space, and H g , we need to define the spatial dimension, and specify the name of the integration variable and the name for the frequency which is the argument of g: For HPfALL, the argument '1' refers to the species label (of which there is only one in this simple model), and the argument '2' refers to the spatial dimension (which needs to be specified because HPfALL will in general include an integral over space). In HGfALL, both arguments '1' denote species labels, of which there need to be two because this quantity relates to correlation functions. The argument 'k' in both functions specifies the name of the variable that denotes the frequency, which appears as an argument to the Fourier transforms. Thus, our framework has derived the following equations for the model: which are exact because this model is a simple immigration-death model, with no nonlinear processes (and hence no correction to mean-field) and no processes that introduce spatiotemporal correlations.
Adding metapopulation dynamics to this dynamic landscape, we now have two entity types: unoccupied patches (type 1), that can be colonised by occupied patches (type 2) (Fig. 2a, b). The dynamics are now represented by the following generator ; where the last two terms represent patch of either type being destroyed (with the same rate μ), and the first one represents the fact that patches are unoccupied (type 1) when they are created. The second and third terms describe the metapopulation dynamics. The second term represents colonisation, where occupied patches turning an unoccupied patch into an occupied patch with kernel c (i.e., the colonisation rate is cðrÞ if the two patches are separated by distance r), and has superscript 'I' (denoting 'Infection') because this is functionally equivalent to type 2 patches 'infecting' type 1 patches. The third term corresponds to spontaneous changes in type (superscript 'CT') from type 2 to type 1, representing extinction of occupied patches.
These processes are represented in Model Constructor by the Mathematica process variable where the argument to the Infection process represents the kernel c, and the argument $c represents its Fourier transform. The equations for the quantities q; g; p can again be obtained using the HQfALL, HGfALL, HPfALL functions. For example, dq 1 =dt, dp 2 =dt and dg 12 ðkÞ=dt are, respectively, given by the outputs of the following commands: These equations are identical to those derived in ref. 15 for the case of an uncorrelated landscape (see the Methods section for the full set of differential equations). Our Mathematica toolbox contains tools for readily computing expressions for these quantities at equilibrium, as well as evaluating them numerically (see Supplementary Note 2 section 2.3).
Simulation code is obtained by the 'Model simulator' toolbox (Supplementary Note 2 section 2.2). The model is specified using a very similar markup to that for 'Model Constructor', except that parameters are given numerical values rather than symbolic names, and kernels have to be assigned specific functional forms. For instance, when colonisation has a top hat kernel with strength cð0Þ ¼ 1:5 and length scale 3, and the other parameters are ðr; μ; eÞ ¼ ð1:1; 0:9; 0:5Þ, the syntax for Model Simulator is specified in file modelCaseStudy1.txt containing the text: Simulations of the model are run by passing, as a commandline argument, the name of the file containing this model definition to the simulation programme ppsimulator we provide (Supplementary Note 2 section 2.2).
Other kernel shapes are implemented within ppsimulator, but for the sake of computational efficiency kernels that represent interactions between individuals are truncated so as to be zero beyond a threshold distance. For example if 'tophat [1.5,3]' were replaced by 'truncatedGaussian [1.5,3]' then the colonisation strength would still becð0Þ ¼ 1:5, but the kernel would be a truncated Gaussian (i.e., zero beyond a certain threshold distance, by default three standard deviations and Gaussian otherwise) with standard deviation 3. For more details, see Supplementary Note 2 section 2.2.
Optimal landscape connectivity. Our first case study concerns the question in conservation biology of how to identify the most important habitat to conserve. A commonly used metric to assess the value of habitat patch to a network, and aid the design of nature reserves, is 'connectivity' 16 . A critical open question is how well connectivity predicts patch occupancy-that is, whether connectivity is a good way to identify valuable habitat. Connectivity is expressed as a weighted sum of the proximity to other habitat patches in the form 16 where h is a 'connectivity kernel', and x i and x j are positions of the ith patch and jth patch, respectively. Commonly, h is taken to be an exponential function h ¼ expðÀ ½distance between patches =λÞ, where λ is usually chosen to equal the species' average dispersal distance. While this is inspired from metapopulation ideas, there is no underlying theory for how well connectivity S predicts habitat occupancy, or whether this kernel is the optimal choice 17 . Our framework allows us to solve this problem and find the bestperforming connectivity measure.
We start with the dynamic landscape patch occupancy metapopulation model on a landscape where habitat patches are ephemeral (Figs. 2a, b, c; 3a) 15 . We showed above in 'Using the framework' how to use our 'model constructor' software to compute the quantities q, p and g for this model. We used this to derive an analytical expression for the correlation between patch occupancy and connectivity S for a general connectivity kernel h (see the Methods section), which can be expressed in terms of spatial moments and therefore in terms of our quantities q, p and g. Considering first an exponential connectivity kernel, we find that the standard choice of the length scale λ is not optimal as it gives a much lower correlation than when λ is 3-4 times larger (Fig. 3b). We further derived an expression for the optimal connectivity kernel, i.e., the one that maximises the x (2) x (1) Log(λ/l c ) Occupancy−connectivity correlation  Fig. 2b. Snapshot a illustrates that patches are more likely to be occupied (blue symbols) rather than unoccupied (pink symbols) when they are close to other patches b. Connectivity-occupancy correlation depends strongly on the ratio of the length scale λ of the connectivity kernel to the length scale l c of the colonisation kernel. Lines show analytical predictions and symbols are simulation results. Colonisation and connectivity kernels are chosen to have the same standard deviation when λ ¼ l c . The highest correlation occurs when λ is 3-4 times larger that l c . The optimal length scale depends only weakly on whether the landscape is static (solid lines, circles) or dynamic (dashed lines, triangles), or whether occupancy P 0 is high (P 0 ¼ 0:667, cyan) or low (P 0 ¼ 0:167, black). The correlation for the exponential connectivity kernel with the optimal choice of length scale is within 10-20% of the horizontal lines, which represent the occupancy-connectivity correlation for the optimal connectivity kernel shapeh ? . The choice λ ¼ λ Ã (vertical lines), where the connectivity kernel has the same standard deviation as the optimal connectivity kernel h Ã , gives a higher occupancy-connectivity correlation than λ ¼ l c . c Occupancy-connectivity correlation decreases when rate of patch turnover μ increases relative to patch extinction rate e. Lines are analytical approximation, and symbols are the results of simulations. Parameters r and μ are chosen so that patch density =1 and mean-field occupancy P 0 is kept constant as μ=e changes; colour denotes P 0 ¼ 0:167 (black); P 0 ¼ 0:333 (magenta); P 0 ¼ 0:667 (cyan). Error bars not shown as standard errors are smaller than the plotting symbols. Parameter values and kernels are given in the Methods section connectivity-occupancy correlation: h ? ðωÞ ¼c ðωÞ cð0Þcð0Þr þ μ 2 ð ÞÀμðe þ μÞcðωÞ : Here and throughout this paper, a tilde represents Fourier transform;h ? ðωÞ is the Fourier transform of the optimal connectivity kernel,cðωÞ is the Fourier transform of the colonisation kernel and ω is a frequency (the Fourier conjugate variable to space). The parameters r; μ and e are defined in the caption of Fig. 2 and have the same meanings as in 'Using the framework' above (i.e., rates of patch creation, destruction and extinction respectively). It is clear from this expression that, in general, this kernel will have a different shape from either the exponential function or the colonisation kernel. However, we find that an exponential connectivity kernel with the optimal choice of λ performs nearly as well (Fig. 3b). Thus, we show that an exponential connectivity kernel is a reasonable choice for landscape design, but we suggest that its ability to identify important habitat would be improved if its standard deviation were set to match that of the optimal kernel (Fig. 3). This means that the connectivity kernel should not just depend on the dispersal properties as is usually assumed 16 , but also on the metapopulation occupancy and landscape turnover as well (Fig. 3b). Moreover, we find that connectivity-occupancy correlation becomes weaker when the landscape becomes more dynamic (i.e., patches are created and destroyed more rapidly, Fig. 3c) 18 . Also, for static landscapes (left side of Fig. 3c), the correlation is stronger when the mean patch occupancy is lower, but this trend is reversed for dynamic landscapes (right side of Fig. 3c). Thus, connectivity is of less use to biological conservation when either occupancy or habitat turnover is high 19 .
Genetic similarity. Our second case study is in the field of population genetics. Molecular ecologists study patterns of genetic similarity and differentiation between populations or across space to infer the underlying processes that cause and maintain that genetic structure. A long-standing problem is to find tractable models where local density is regulated so that population dynamics are stable 20 . Many ingenious solutions to this problem include imposing constant local density artificially [21][22][23] . However, since the aim of these models is to relate pattern to process, a much preferable solution would be to do this via local density-dependent population dynamics. Such models are regarded as too 'hopelessly complicated' 24,25 to study using standard methods: even in the simplest models, it is highly nontrivial to obtain analytical expressions for genetic similarity, and the incorporation of more realistic processes (e.g., such as selection) poses even greater difficulty 24,25 . Our framework overcomes these obstacles, allowing us to consider local densitydependent population regulation explicitly, and to derive analytical expressions for how genetic similarity varies in space. We illustrate this with an example including local competition and limited dispersal, selection and mutation between alleles in neutral and non-neutral loci, and also adaptation to heterogeneous environments, which is important for example for tropical forests 22 .
The model involves two habitat types, and haploid individuals with four different genotypes (one selective and one neutral locus, two alleles at each). Habitat patches appear independently as a Poisson process with rate κ and vanish with rate μ. All individuals have density-independent mortality (rate m) and densitydependent mortality with the same interaction kernel c between all genotypes. Birth is density-independent and offspring are distributed relative to their parents with kernel d, but the fecundity depends on habitat availability (via the patch kernel r) and the match between the habitat type and the genotype of the individual via the function ϕ τ ¼ f 0 Á ð1 ± τÞ, where f 0 is a constant base fecundity and τ is the strength of selection. Therefore, each habitat patch a distance x away from the individual contributes rðxÞ Á f 0 Á ð1 þ τÞ to the fecundity of that individual if allele in selective locus matches the type of habitat, and rðxÞ Á f 0 Á ð1 À τÞ otherwise. Offspring inherits their parent's genotype (i.e., reproduction is clonal) unless there is a mutation, which happens independently at each allele at rate ν per birth and specified by the function ψ ν . All model components are shown in Fig. 4a, and a typical snapshot of the dynamics is illustrated in Fig. 4b.
We are interested in the similarity functions F n ðxÞ and F s ðxÞ at the neutral and selective loci, respectively, i.e., the probability that two individuals separated by distance x have the same allele at the locus in question. These can be computed from the probability densities to find specific individuals at this spatial separation. Using the following function i;j ðx i ; x j Þ that in turn can be expressed in terms of the quantities q, p and g in the perturbation expansion 11 , i.e., k ð2Þ i;j ðxÞ ¼ q i q j þ ϵ d ðg i;j ðxÞ þ q i p j þ p i q j Þ þ o ϵ d À Á : Using our framework, we obtained analytical expressions for how genetic similarity depends on the distance between individuals (Methods) which show good agreement with numerical simulations, Fig. 4c. Genetic similarity at the selective locus is always higher than at the neutral locus, and also extends further in space due to having different length scales l n and l s (Fig. 4c), where l s ; l n are, respectively, the length scales of genetic similarity at the selective and neutral loci. Our analytical expressions for genetic similarity (shown in the Methods section) enable us to understand what controls the differences in these length scales. The results become particularly simple under the plausible assumption that mutation is rare (ν ( 1, see the Methods section): is the ratio of the total density of individuals N e to the total density of habitat patches, M ¼ 4νq h f 0 e rð0Þ e dð0Þ is the per-capita creation rate of mutants in the population, μ is the turnover rate of habitat patches and r ¼ ð1 þ μ=MÞ À1=2 . We see that the ratio of these length scales is controlled by just two parameter combinations: μ=M, which quantifies how ephemeral the landscape is relative to the appearance of mutants, and τ 2 R 2ν , which will be numerically large when ν ( 1 unless selection is very weak or habitat patches vastly outnumber the organism. Similarity at the selective locus will resemble that at the neutral locus if mutation is common or selection weak, but will extend up to twice as far when habitat turnover is slow relative to mutation (Fig. 4d).
Optimal foraging. Our third case study originates from the field of movement ecology, in which a long-standing challenge is to understand the causes and consequences of movement behaviour on foraging efficiency 26,27 . In environments with patchily distributed resources, commonly observed behaviours involve slow foraging movements within patches, interrupted by fast exploratory movements between the patches 28 . Individual-based models are needed to account for stochasticity and resource heterogeneity in foraging 29 , but there are several processes at work so it is very difficult from simulations to understand the key evolutionary drivers and ecological consequences of variation in foraging behaviour. Our framework solves this problem by providing simple analytical expressions that quantify the ecological factors determining the optimal foraging strategy.
We used the graphical model description (Fig. 5a) to construct a model in which aggregates of resources are continuously generated at new locations, while existing resource units decay, resulting from the consumer's point of view in an unpredictable resource distribution. The model consists of resources (targets), which are generated in clusters as follows. Target 'generators' appear as a Poisson process with intensity b, and disappear at rate h. During their lifetime, target generators create targets at rate hλ at a distance from the generator determined by kernel r (normalised sorð0Þ ¼ 1). We assume h ! 1 such that hλ is constant, so that each generator instantaneously creates a Poisson distributed number of targets with mean $λ). Targets vanish at rate μ. The forager moves as a jump process with jump kernel c at rates m S and m F , respectively, when in slow or fast mode. Foragers consume targets at rate γ with kernel f . We assume the consumer has evolved to switch to a fast movement mode (at a rate we denote by α) when it does not encounter new resources, and to switch back to slow movement mode after encountering a resource unit (Fig. 5a, b).
Our mathematical formalism yields a simple closed-form expression for the mean consumption rate (Fig. 5c), ρ % ρ 0 þ ϵ d ρ 1 , where ρ 0 ¼ γβ=μ is the mean-field consumption rate, β ¼ λb; and the first-order correction to ρ is x (2) x (1) Distance between individuals, x Similarity function F(x)  Fig. 5c). This is because a consumer that remains continuously in the slow mode is not efficient in finding new resource aggregates, whereas a consumer that remains continuously in the fast mode misses the opportunity of elevated resource availability within aggregates. We also find a simple expression for the optimal switching rate: , λ is the mean number of resource items per cluster, and I ¼ I 1 =I 2 (0 < I < 1) quantifies the spatial scale of the resource detection process relative to the resource cluster size. Thus, although we started with a model with nine processes, the optimal foraging strategy is only determined by the two parameter combinations ρ 0 μ and λI (Fig. 5d)-for example, an increase in resource acquisition rate has the same effect on the switching rate as the same proportionate increase in resource production rate. It therefore pays to spend a long time in the slow foraging mode (small α Ã ) if the resources are highly aggregated (large λ) or if the consumer's resource detection efficiency (γ) is low, Fig. 5d. However, if fast movement requires more energy than the slow mode, we find that α Ã ¼ 0 when these energetic costs per unit time exceed a threshold (see the Methods section). This threshold value shows that it may be better to stay put and wait for resources to be generated locally, rather than to roam widely, when the resource consumption rate (ρ 0 ) is low, resource clusters contain few items (λ small), or the spatial scale of clusters and/or resource detection is large.

Discussion
Our framework is game-changing because (i) it makes mathematical approximations available for a very wide class of individual-based models, and (ii) this approximation is available to the non-specialist because it requires no laborious derivations on the part of the researcher. Since we have derived expressions for general reactant-catalyst-product model, and implemented these in Mathematica, a researcher interested in a particular model has only to convert their model into the markup understood by our software (see the Methods section and Fig. 2a-c) and the approximation is generated for them automatically. Without our method, mathematical results for a new model would only be available to a researcher with the skill and patience to spend weeks performing algebra without error. Our software also generates code for simulating the same model. The mathematical x (2) x ( 5 Optimal foraging model. a shows a graphical definition of the model and b shows a simulation snapshot, illustrating a forager (large cyan circle) that has found a resource aggregate, and has switched to the slow movement mode to consume resources (small magenta circles). c shows how mean resource consumption rate ρ varies from a totally stationary searcher (α ¼ 0; continuous horizontal line) through a searcher following the optimal rate α ? (vertical line) to a highly mobile searcher (α ! 1, dashed line). In c, the dots with error bars show simulation results, whereas the lines follow the analytical prediction of ρ % ρ 0 þ ϵ d ρ 1 (see Eq. (1)). The contour lines in d show the optimal switching rate α ? as a function of the rescaled detection rate γβ=μ 2 of the searcher and the clustering level λI of the targets. Error bars represent three standard errors. Parameters are shown in the Methods section validity of our approach is proved in Supplementary Note 1, and the validity of the approximation in these examples is shown by its close agreement with simulations in Figs. 3-5. Since our framework is based on a perturbation expansion, it is guaranteed to give a good approximation when the length scale of interaction kernels is large enough compared to the separation between individuals in the system 15 . One way to assess the accuracy of the approximation is to compare p to q: if the first-order correction p is comparable in magnitude to the mean-field density q, then approximation may not be very accurate for those parameter values because higher order corrections are likely to be important. Our examples illustrate the conceptual and methodological advantages of analytical methods over simulations. In the first case study, we have a general expression for the best-performing connectivity kernel, valid for any colonisation kernel, whereas simulations can only explore specified examples and functional forms. In the third case study, we obtained directly an expression for the optimal foraging strategy, whereas in a simulation study the optimal strategy would have to be searched for separately for each parameter combination. We found in the second and third examples that the behaviour depends on a small number of parameter combinations, which effectively reduces the dimension of the parameter space of the problem. This means that the model behaviour can be explored more economically, and explained more intuitively, than in a simulation study in which each set of parameter values needs to be investigated independently. While simulations will always play an important role in studying spatial and stochastic models, our mathematical framework will facilitate much deeper understanding by making analytical results readily available for models that were previously thought to be too complex.
There are already platforms that allow simulation of spatial, stochastic individual-based models specified in a user-friendly way 10,30 , but ours is the first framework that automatically performs mathematical analyses as well. We have formulated, analysed mathematically and written computer code for the general model containing any combination of reactant-catalyst-product processes. This general model can be used to describe a very large array of processes, so that a wide variety of models of interest in ecology and evolution are straightforward to construct and analyse using the unified framework presented here. Thus, our framework provides a long-sought analytical tool 8 with which many classical questions in ecology, evolution and other fields can be revisited to understand the role of space and stochasticity.

Methods
Derivation of the analytical framework. Here, we give a summary of the approach and key results of our analytical framework; full details of the calculations are given in the Supplementary Notes.
We begin by classifying the participants in demographic processes into three types of individuals: (i) reactants (that are destroyed by the process); (ii) products (that are created by the process) and (iii) catalysts (that are unaffected by the process but whose presence affects the rate at which it occurs). A given demographic process is therefore characterised by integer numbers denoted by P; R and C that give number of agents in the groups 'reactants', 'products' and 'catalysts'. The group 'products' is characterised by the set P consisting of pairs of the type of agent i and the location x (from dÀdimensional space x 2 R d ) assigned to the agent: P ¼ ffi 1 ; x i 1 g; fi 2 ; x i 2 g; ; fi P ; x i P gg: Similarly, reactants are described by R ¼ ffj 1 ; x j 1 g; fj 2 ; x j 2 g; ; fj R ; x j R gg; and catalysts are described by C ¼ ffk 1 ; x k 1 g; fk 2 ; x k 2 g; ; fk C ; x k C gg: We introduce the function r P; R; C ð Þ , which depends on types of agents and on distances between locations of agents used in P, R and C. The function r determines interactions between agents by specifying those coordinates and indices which are the same in P, R, C. For example, the 'Change in type' process where agents of type j 1 transform with rate r 0 into agents of type i 1 without changing their locations is described by: P ¼ ffi 1 ; x i 1 gg; R ¼ ffj 1 ; x j 1 gg, C ¼ +, and r P; R; C ð Þ¼r 0 δðx i 1 À x j 1 Þ. For the purpose of this paper, it is convenient to work with cumulants u instead of moments, as cumulants carry the same information as moments but lead to a more simple perturbation expansion 11 . We define the operator Q Δ that determines the full hierarchy of equations for cumulants: Here, all orders of cumulants are denoted by uðη; tÞ, where η denotes any finite number of points in dÀdimensional continuous space R d occupied by agents of a specified type. For example, for a one-point configuration for an agent of type m one has η ¼ fx; mg, and the expression uðη; tÞ gives a first-order cumulant u ð1Þ m ðx; tÞ which is the density of agents m at location x, denoted as D m ðx; tÞ ¼ u ð1Þ m ðx; tÞ. For a two-point configuration η ¼ fx; m; y; ng one obtains a second cumulant uðη; tÞ ¼ u ð2Þ m;nðx; y; tÞ that determines spatial covariance denoted Cov m;n ðx; y; tÞ ¼ u ð2Þ m;nðx; y; tÞ. Following ref. 11 , we made the following steps in the derivation (for detailed derivations see Supplementary Note 1): (1) using notions of locally finite configurations, transform an agent-based model description into the dynamics of spatial moments; (2) using the correspondence between correlation functions and cumulants, obtain operator Q Δ for the evolution of cumulants. We obtained: where we used the following operations D ðmÞ x m and D yðmÞ x m , that can be considered as related to creation and annihilation operators: where we used the operation Ã that is defined below, also see ref. 11 . First, let Γ 0 denotes the set of all finite subsets η of R d . For any functions u; v on Γ 0 , the Ã operation is defined as: where the symbol F denotes a disjoint union defined as X thus, the symbol Ã defines a convolution. Using these notations, the function ðexp Ã uÞðηÞ is The function exp ÃÀ1 u is the inverse with respect to the Ã convolution, i.e., as was shown in e.g., 31 , for any u with uð+Þ ¼ 0, there exist a function exp ÃÀ1 u such Þ Ãn ðηÞ: The system of Eq. (2) cannot be solved exactly, because an equation for a cumulant of any order depends on higher order cumulants. In our derivations we follow the method from 11 , which assumes that interaction between individuals is long ranged; details of derivations are presented in Supplementary Note 1, sections 1.1-1.3. Using this exact expression of the full hierarchy of moment equations for the general model, we apply the perturbation scheme 11,13 to derive a controlled approximation that gives asymptotically exact results when agents interact over large enough scales. It is assumed that spatial interactions between individuals depend on their separation x (in d-dimensional space) according to interaction 'kernels' of the form ϵ d f ðϵxÞ where 1=ϵ is the typical length scale of interactions. In models with long-ranged interactions (i.e., where typical scale of interaction is large), the parameter ϵ is a small parameter that can be used to develop a perturbation expansion. It is this small parameter ϵ that makes this approximation controlled: the smaller ϵ, the better leading terms in a perturbation expansion describe the exact solution.
We showed mathematically that for a model made up of a set general process defined by P, R, C and the interaction function r κ ðP; R; CÞ, the mean densities and spatial covariance (including autocovariance) satisfy the following expansion: density of species m ¼ q m ðϵx; tÞ þ ϵ d p m ðϵx; tÞ þ oðϵ d Þ; spatial covariance between species m and n ¼ ϵ d g m;n ðϵx 1 ; ϵx 2 ; tÞ þ oðϵ d Þ; where dq m ðx; tÞ dt ¼ X For a general (arbitrary) product-reactant-catalyst process κ, we derived the corresponding contributions H ðκÞ q m , H ðκÞ p m and H ðκÞ g m;n to the differential equations for q m ; p m and g m;n , see Eq. (2) in Fig. 1e. The expression H ðκÞ q m for the mean-field density q m ðx; tÞ of an agent of type m for an arbitrary spatially heterogeneous product-reactant-catalyst process is shown below: Expression for H Using expressions H ðκÞ q m ðx; tÞ, H ðκÞ p m ðx; tÞ and H ðκÞ g m;n ðx 1 ; x 2 ; tÞ, it is straightforward to show that values of quantities q m ðϵx; tÞ, ϵ d p m ðϵx; tÞ and ϵ d g m;n ðϵx 1 ; ϵx 2 ; tÞ calculated using the given kernels a are identical to the values of quantities q m ðx; tÞ, p m ðx; tÞ and g m;n ðx 1 ; x 2 ; tÞ calculated using the scaled kernels a ϵ , a ϵ ðx À yÞ ¼ ϵ d aðϵðx À yÞÞ.
Computer code. The analytical software 'Model Constructor' is developed in Mathematica 32 . The code first requires a user to specify the model in terms of products, reactants and catalysts, and the interactions. Then, using Fourier transform and a simplifying assumption of translational invariance, the code uses the input data to obtain the analytical expressions for H ðκÞ q m ðx; tÞ, H ðκÞ p m ðx; tÞ and H ðκÞ g m;n ðx 1 ; x 2 ; tÞ for a given system in 1D, 2D or 3D infinite space. Also, the code can write down these analytical expressions in real space and in Fourier space into .tex file. The analytical software has been checked to reproduce the results for 15 specific processes (see Supplementary Notes) that were derived analytically.
The expressions produced by Model Constructor assume translational invariance, i.e., that the agents do not interact with a spatially heterogeneous extrinsic environment, and the expressions are averaged over initial conditions that are translationally invariant. This is not due to a limitation of the underlying mathematical framework, as Eqs. (4), (5), (6) do not make this assumption. However, in the absence of translational invariance the equations for g and p involve two spatial co-ordinates rather than one, and as a result are much more challenging to solve both analytically and numerically. Any analytical solution would require exploiting whatever symmetries are present in the initial condition, which depends on the details of each individual case. While we are studying some specific cases without translational invariance (work in progress), software that allows spatial heterogeneity to be specified in a general way is beyond the scope of this paper.
The simulation software 'Model simulator' is a C-programme for simulating continuous-time point processes. Each point is associated with a coordinate and a discrete species attribute. Points are located either on 1D or 2D torus space. The user defines the set of processes and the initial configuration after which the simulator runs the Gillespie algorithm 33 in such a way that the information of point locations are taken into account, i.e., the system is not assumed to be wellmixed. The state of the configuration can be outputted at user-defined constant time intervals. Auxiliary R-functions are provided for calculating summary statistics and creating figures and animations based on the simulation. Input for the simulator is given by means of text files and few command line arguments. Output of the simulator is written in text files. The simulation software has been checked to replicate the results that we obtained earlier using more specific implementations (e.g. ref. 11 ) that were fully independent of the current implementation (e.g., coded in different programming languages).
In all systems studied in this paper results of both types of software, analytical and simulation, match each other in the way that the mathematical theory suggests. The detailed derivations of underlying mathematical expressions, and the detailed tutorials for analytical and simulation software are shown in Supplementary Notes.
Case studies. Model Constructor was used to obtain the analytical results used in the case studies, and for representative parameter values these were compared with simulation results obtained using Model Simulator. The toolboxes and files containing the markup for the case studies are included at Figshare, https://doi.org/ 10.6084/m9.figshare.9633161. We give here the essential details for obtaining the case study results used in the paper; further details are given in Supplementary Notes 3, 4 and 5.
Optimal landscape connectivity model. For case study 1, the model is defined verbally in Fig. 2a and the processes are listed in Fig. 2b.
Perturbation equations: The section 'Using the framework' in the main paper contains a full description of the steps required to derive the following expressions for q i and g ij (p i is not needed for the rest of the calculation): dq 2 dt ¼ Àeq 2 À μq 2 À q 1 q 2c ð0Þ dg 11 dt ¼ À2μg 11 þ 2eg 12 À 2g 11 q 2c ð0Þ À 2g 12 q 1c ðωÞ dg 12 dt ¼ ðg 11 Àg 12 Þq 2c ð0Þ þ ðg 22 Àg 12 Þðe À q 1c ð0ÞÞ À 2μg 12 À q 1 q 2c ðωÞ , and solving (which can be performed by Mathematica), we find the following equilibrium densities and correlation functions that will be needed later: . These averages were expressed as spatial moments (Supplementary Note 3 section 3.4), which were expanded using the perturbation expansion to leading order in ϵ d were retained. Expressions for q and g above were then substituted into the expression for R to give Proof that landscape turnover weakens connectivity-occcupancy correlation: From Eq. (7), we compute ∂R=∂ρ, keeping the mean patch occupancy P 0 and mean patch density Q 0 constant: This expression is always negative. Kernel that maximises correlation: We take a variational approach and writẽ hðωÞ ¼h ? ðωÞ þ νδðjωj À jω 1 jÞ, where δ is a Dirac Delta function. If the correlation is maximal when h ¼ h ? , then ∂R=∂ν ¼ 0 for all ω 1 , which leads tõ Variance and standard deviation of the optimal kernel: For any rotationally symmetric function MðxÞ with variance the Taylor expansion of the Fourier transform of M is where a is a number that depends on spatial dimension, but not on M. Thus, cðωÞ ¼cð0Þ 1 À aV c ω 2 þ oðω 2 Þ ð Þ , and expandingh ? gives Therefore, the variance of the optimal connectivity kernel h ? is Thus, the standard deviation of the optimal kernel is λ Ã ¼ l c ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðcð0Þr þ μ 2 Þ=ðcð0Þr À eμÞ p , where l c is the standard deviation of the colonisation kernel.
Genetic similarity model. Full operator L defining the model: Denoting two habitat types as species s 1 and s 2 , and denoting species with four different genotypes as s 3 ; s 4 ; s 5 and s 6 , the graphical definition of the model presented in Fig. 4a leads to the complete definition of the model by the following full operator L: where abbreviations have the following meaning: IM stands for Immigration, D for density-independent death, DE for death by external factor, BTF for birth to another type by facilitation (these processes are defined in Supplementary Note 1 section 1.3).
Probability densities: Using notations introduced in the main text, the probability density for any two individuals at locations y and x þ y is considered as i;i ′ ðy; x þ yÞ, where the sum is taken over all possible pairs ði; i′Þ of genotypes (there are four possible genotypes: fA1; A2; B1; B2g). Probability density for individuals at y and x þ y with allele j in neutral locus (i.e., j ¼ 1; 2) is defined as k Analytical expressions of similarity functions: We use Model Constructor to obtain expressions for q, p and g for the model, and hence expressions for the leading terms in the expansion for the two-point correlation functions (second spatial moments) k ð2Þ i;j . To first-order in ϵ d , and from now on assuming we are in spatial dimension d ¼ 2, the similarity functions F n ðxÞ and F s ðxÞ take the form where J 0 ðωÞ is Bessel function of the first kind of order zero and q Ã ¼ 2f 0 κ e dð0Þe rð0Þ À mμ 4μe cð0Þ is the mean-field density for any single genotype.
In the limit of very rare mutations: ν ! 0, the denominator in the integrals above is dominated by the behaviour for ω small. In small-ω limit in 2D one obtains e dðωÞ % e dð0Þ 1 À 2π 2 σ 2 ω 2 ð Þ þ oðω 3 Þ, where σ 2 is the variance of the kernel d. As a result, we have N e ¼ 4q Ã , and we have assumed x ( σ; K 0 is a modified Bessel function of the second kind of order zero. Similarly, when ν ! 0, the similarity at the selective locus is F s ðxÞ % F n ðxÞ þ τ 2 8νq h πσ 2 K 0 ðx=x n Þ À K 0 ðx=x s Þ μ=M ; x s ¼ x n ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þ μ=M p ; where q h ¼ κ=μ is the density of habitat patches of a single type, and M ¼ 4νq h f 0 e rð0Þ e dð0Þ is the rate with which an existing individual produces a new mutant. Derivation of length scales: We define the length scale for a similarity function as the mean distance over which it decays to its asymptotic value at large distances lim x!1 FðxÞ ¼ 1 2 À Á : According to this definition, using the small-ν expressions for F n and F s in terms of Bessel functions, we obtain the expression for the ratio of length scales presented in the main text. Parameters used in Fig. 4: b ν ¼ 0:01, κ ¼ 0:05, μ ¼ 1, τ ¼ 1, f 0 ¼ 1, m ¼ 1; rðxÞ ¼ 3:5 Tophat(x; 2); dðxÞ ¼ 3 Tophat(x; 2); cðxÞ ¼ 1:25 Tophat(x; 3=20), where Tophat(x; R) is a normalised tophat kernel with radius R. c dðxÞ ¼ rðxÞ ¼ Tophatðx; 6Þ, and cðxÞ ¼ Tophat(x; 20).
Optimal foraging model. The model is specified in section 'Optimal foraging' in the paper and Fig. 5a where the entity types are defined as follows: 1: resource target generators; 2: resource targets; 3: slow-moving consumer; 4: fast-moving consumer; the abbreviation BT stands for birth to another type, J for jump, CTC for change in type by consumption, other abbreviations were already explained above (the full definition is presented in Supplementary Note 1 section 1.3, see also Supplementary Tables 1-3). Parameter h determines the rate at which target generators disappear; we will study the case where clusters of targets are created simultaneously by taking the limit h ! 1. Since each target generator creates resources at rate hλ, and has mean lifetime 1=h, it generates a Poisson distributed number of resource targets with mean λ. This is a convenient method for generating clusters of individuals from an underlying process that generates individuals; a similar approach was used in ref. 15 to generate clusters of habitat patches.
Our interest is in the equilibrium rate ρ Ã by which each searcher consumes targets. This can be computed by solving the stationary state of the system under the condition that the density of searchers, denoted here by A, is constant. This yields ρ Ã ¼ ðβ=k Ã 2 À μÞk Ã 2 =A, where β ¼ λb is the appearance rate of targets, k Ã 2 is the stationary density of targets and the expression in brackets is the rate at which a randomly selected target disappears due to consumption. As we are not interested in resource competition among multiple searchers, we consider the case of a single searcher, obtained technically by taking the limit A ! 0. Note that, in the main text, we have set ϵ ¼ 1, i.e., the length scale of the kernels f and r is the true biological one and has not been rescaled by ϵ.
Resource consumption rate: Using Model constructor, the stationary density of targets in the mean-field approximation is k Ã 2 % q Ã 2 ¼ β=ðAγ þ μÞ, and thus the consumption rate of targets by a single searcher is ρ Ã % ρ 0 ¼ γβ=μ independently of the parameter α.
In the first-order approximation, the stationary density of targets is k Ã 2 % q Ã 2 þ ϵ d p Ã 2 . The general expression for p Ã 2 is given using 'The model constructor' toolbox: p Ã 2 ¼ Using 'The model constructor' toolbox to calculate the limit h ! 1 (after which the results no longer depend on h) and to extract the leading term in the limit A ! 0; which is achieved by the following Mathematica commands: