Territorial Dynamics and Stable Home Range Formation for Central Place Foragers

Uncovering the mechanisms behind territory formation is a fundamental problem in behavioural ecology. The broad nature of the underlying conspecific avoidance processes are well documented across a wide range of taxa. Scent marking in particular is common to a large range of terrestrial mammals and is known to be fundamental for communication. However, despite its importance, exact quantification of the time-scales over which scent cues and messages persist remains elusive. Recent work by the present authors has begun to shed light on this problem by modelling animals as random walkers with scent-mediated interaction processes. Territories emerge as dynamic objects that continually change shape and slowly move without settling to a fixed location. As a consequence, the utilisation distribution of such an animal results in a slowly increasing home range, as shown for urban foxes (Vulpes vulpes). For certain other species, however, home ranges reach a stable state. The present work shows that stable home ranges arise when, in addition to scent-mediated conspecific avoidance, each animal moves as a central place forager. That is, the animal's movement has a random aspect but is also biased towards a fixed location, such as a den or nest site. Dynamic territories emerge but the probability distribution of the territory border locations reaches a steady state, causing stable home ranges to emerge from the territorial dynamics. Approximate analytic expressions for the animal's probability density function are derived. A programme is given for using these expressions to quantify both the strength of the animal's movement bias towards the central place and the time-scale over which scent messages persist. Comparisons are made with previous theoretical work modelling central place foragers with conspecific avoidance. Some insights into the mechanisms behind allometric scaling laws of animal space use are also given.


Introduction
Understanding the mechanisms behind animal territoriality is of great importance to many areas of ecology [1], from conservation biology [2] to epidemiology [3] to predator-prey dynamics [4]. A species is called territorial if each animal, or group of animals, constructs and defends a region of space from conspecific neighbours or possible intruders. Maintaining a territory relies on the animal's ability to exclude conspecifics from the area it occupies. Since the animal needs to spend time moving inside its territory to carry out vital activities such as foraging, continuous monitoring of territory boundaries is not possible. Therefore many animals have evolved mechanisms whereby their territory is identified by visual, auditory or olfactory signals [5], thereby obviating the need for constant border patrolling.
In this paper we focus on a model where the signals are olfactory (scent marks). It is based on an agent-based model of socalled territorial random walkers, first introduced in [6], where animals are modelled as lattice random walkers that deposit scent as they move. The scent is only active for a finite amount of time, the socalled active scent time, and if a lattice site contains active scent, no other animal may move there. As a result, the terrain naturally subdivides into territories demarcated by the absence of foreign scent. Territories each have a boundary and if two boundaries are juxtaposed, a border is formed. These borders never settle to a stable state. Instead, they continually ebb and flow, albeit at a much slower rate than the animals move. Specifically, the border movement is subdiffusive (i.e. the variance of the border position's probability distribution increases sublinearly with time) since the territories are undergoing an exclusion process [7,8], whereas the animals move diffusively.
Here, we study a modified version of the territorial random walk model where animals are random walkers with an attraction towards a central place, such as a den or nest site where the animals return occasionally [9], or a core area where animals tend to spend most of their time [10]. Similar to the original territorial random walk model, territories emerge whose borders are continually fluctuating. However with central place attraction, the mean square displacement (MSD), i.e. the variance of the border position's probability distribution, tends towards a finite value, as confirmed by stochastic simulations. This causes stable home range patterns to emerge from the territorial dynamics.
To understand better the precise nature of the emergent home range patterns, we compare stochastic simulations of the manybodied non-Markovian central place attraction model with an analytic approximation, following [11,12]. This exploits the time-scale disparity between the rate of animal movement and the slower, subdiffusive territorial borders, to construct an adiabatic approximation for the joint probability distribution of the animal and territory border positions. The model is solved exactly in both 1D and 2D and the resulting marginal distribution for an animal's position allows the macroscopic properties of home range size and overlap to be related to the microscopic details of the animals' movement and interaction processes. In particular, our analytic expressions can be used to infer the longevity of olfactory messages purely by examining data on animal space use. Furthermore, since various properties of space use are predicted to scale allometrically [13], our theory can also be used to give insights into the mechanisms behind these scaling laws. Our results are compared with previous approaches to modelling conspecific avoidance with reaction-diffusion formalisms [9].

Agent-based simulations of territorial central place foragers
Monte Carlo simulations of the territorial random walk system were performed in both 1D and 2D where each animal has a bias of moving towards a central place (CP) (see Methods for details). The MSD of the territory border eventually reached a saturation value that depended on both the strength of attraction towards the CP and the dimensionless quantity representing the time it takes for an animal to move around its territory, D is the animal diffusion constant and r the animal population density. The parameter a~vL=D was used to measure the dimensionless strength of CP attraction, where v is the drift velocity towards the CP and L the distance between CPs of two adjacent territories.
For a fixed a, the amount of border movement arises from the ratio of the active scent time to the diffusive time, which is Z 1 in 1D or Z 2 in 2D (figure 1). Increasing a has the effect of reducing the animal's tendency to move into interstitial regions and claim extra territory. This causes the borders to move less on average as each animal keeps to a small core area well within its territory most of the time. Consequently, when plotting the MSD saturation value against Z 1 or Z 2 , we see that the curves for higher values of a lie below those for lower values (figure 1).
Dynamics of a central place forager within its territory: a reduced analytic model By taking into account the fact that the border movement is much slower than that of the animal, we employed an adiabatic approximation to calculate the probability distribution of an animal inside its fluctuating territory borders (see Methods). The simulated animals are identical, so it is sufficient just to model one animal. Since the MSD of each territority border saturates at long times, the animal probability distribution reaches a steady state.
Movement in 1D. By fixing the CP at the origin for simplicity, we calculated the steady state 1D dimensionless joint probability density function P P 1D ( x x, L L 1 , L L 2 ) of the left (right) border being at dimensionless positions L L 1~L1 =L ( L L 2~L2 =L) and the animal being at position x x~x=L at long times, where L 1 , L 2 and x are dimensional parameters and L is the distance between CPs of adjacent territories (see figure 2 for an illustration and table 1 for details of notation). This is (see Methods for derivations, here and elsewhere) where H(z) is the Heaviside step function (H(z)~0 if zv0, H(z)~1 if z §0), g g 1 ( L L 1 ) (resp. g g 2 ( L L 2 )) is the probability of the dimensionless territory border position x b (x b ) on the dimensionless parameters a~vL=D and Z 1~TAS =T D2 (Z 2~TAS =T D2 ) from stochastic simulation output. The notation S . . . T denotes an ensemble average over stochastic simulations. The border movement is non-dimensionalised by dividing by L, the average distance between central places of adjacent territories. Panel (a) shows output from 1D simulations and panel (b) from 2D simulations. The best-fit lines for the 2D plots are  Each animal moves diffusively with a constant drift towards the CP and constrained to move between the two territory borders to its immediate right and left. The position of the animal studied in the main text is denoted by x. The animals at x A and x C are drawn purely for illustrative purposes. In the Results section, B is assumed to be at 0 and B{A~C{B~L. doi:10.1371/journal.pone.0034033.g002 Table 1. Notation glossary.

Symbol Model Dimension Explanation
T AS S1,S2 T Active scent time: time for which a scent mark is avoided by conspecifics.
Positions of the left and right borders.
K S1,S2,A1,A2 Rate at which territory sizes tend return to the mean size.
x A1 S Position of the animal in 1D.
(r,h) A2 (S, none ) Position of the animal in 2D polar coordinates.
s A2 S Radius of the territory.
F S1,S2 T {1 Rate of jumping to the nearest neighbour.
p S1,S2 none Probability of an animal moving towards its CP next jump.
a S1,S2,A1,A2 none Normalised drift velocity a~vL=D. Glossary of the various symbols used throught the text. The second column details whether the symbol is used in the 1D simulation model (S1), the 2D simulation model (S2), the 1D analytic model (A1) or the 2D analytic model (A2 distribution of the left (right) border and h h 1D ( x xD L L 1 , L L 2 ) is the probability distribution of an animal being at position x x, given that the borders are at L L 1 and L L 2 . The border probability distributions are given by the following expressions g g 1 ( where k~K=(cL 2 ), K is the border generalised diffusion constant, representing the amount the borders tend to move (see methods), and c is the rate at which the territory size tends to return to the expected average value. To visualise these distributions, notice that when k is relatively large, the nw1 terms are negligible, so that g g 1 ( The probability distribution of an animal being at position x x, given that the borders are at L L 1 and L L 2 , is the following normalised version of a Laplacian distribution with average displacement 1=a Movement in 2D. In 2D we assumed that the territory is circular, the CP is at the centre of the circle and the border movement is modelled by fluctuations in the territory radius. The steady state dimensionless joint probability density function P P 2D ( r r,h, s s) for the territory and the animal at long times is P P 2D ( r r,h, s s)&½H( s s){H( s s{1)½H( r r){H( r r{ s s) g g( s s) h h 2D ( r r,hD s s), ð5Þ where s s~s=L is the dimensionless radius, ( r r,h)~(r=L,h) are the dimensionless polar coordinates of the animal, s is the radius and r is the radial component of the animal's coordinates. Here, g g( s s)~1z2 is the probability distribution of the territory radius and is the probability distribution of the animal being at position ( r r,h) inside a territory of radius s s, where E n (z) is the special function defined by E n (z)~Ð ? 1 ds½ exp({zs)=s n . The limit as z?0 of E n (z) is infinite for n~1 and finite for nw1. For large z, E n (z)*e {z =z so the limit as z?? is 0 for every n.
The marginal distribution of the animal Equations (4) and (7) enabled us to calculate the marginal probability distribution of the animal in both 1D and 2D scenarios, where the territory can be anywhere else, by integrating over all possible positions for the territory border. In 1D the dimensionless marginal distribution of the walker at long times is and in 2D, this is The effects that the two parameters a and k have on the marginal distribution (figure 3) can be characterised by observing that a tends to govern the shape of the density function towards the centre of the territory, whereas k governs the degree to which the distribution tails off sharply (high k) or with a shallow gradient (low k).
Expressions (8) and (9) are directly compared with those measured from territorial central place forager simulations. It turns out that the 1D case gives an excellent agreement for all parameter values we tested (figures 3(a-d)). In 2D, a qualitatively close fit is attained only when k and a are sufficiently low. For higher k or a the borders are moving too fast for the adiabatic approximation to be accurate (e.g. figure 3h). However for lower k and a, the terrain contains very little interstitial area at any point in time, so the territories are forced to tesselate the plane. Therefore they each form more of a hexagonal than a circular shape (e.g. figure 3e).

Obtaining active scent time from animal position data
To make use of the present theory, data must be gathered over a sufficiently long period for the animal MSD to saturate. For certain species, the saturation value fails to be reached during the maximal biologically relevant time-window. Male red foxes (Vulpes vulpes), for example, may spend parts of the autumn and winter moving outside their territories to cuckold or disperse [14], so territorial dynamics can only be measured reliably from the animal positions during spring and summer when the males tend to stay within their territories. During those two seasons, the tendency to return to the CP is so weak that the animal MSD continues to increase slowly, never settling [6]. In such cases, it is necessary to use methods developed in [6] to analyse the animal territorial system.
However, if the animal MSD does saturate then the marginal distribution (9) can be fitted to the non-dimensionalised distribution of position locations from the data in order to obtain the parameters a and k. From the theory, the saturation MSD D s s 2 (k)~S( s s{S s sT) 2 T of the territory radius can then be derived from the equation which allows the MSD of the territory radius s s to be computed from k. The MSD of s s is the analogue, in the analytic model, of the dimensionless territory border MSD Dx 2 b~S (x b {Sx b T) 2 T from the simulation model, so we equate Dx 2 b and D s s 2 (k). By using the appropriate curve from the simulation output (figure 1b) related to the value of a calculated from the data, a value for Z 2~TAS Dr is obtained, from which T AS can be derived.
In summary, the active scent time may be obtained from data on animal locations by using the following programme.
1. Fit equation (9) to the data in order to obtain values of a and k. 2. Use this value of k to find the theoretically expected saturation value of the MSD D s s 2 (k) via equation (10). 3. Note that D s s 2 (k) from the analytic model is equal to Dx 2 b from the simulation model. 4. Identify the best-fit line from figure 1b for the value of a found in step 1. 5. Use this line, together with the value of Dx 2 b from step 3, to determine the Z 2 -value from figure 1b for the data being studied. 6. Assuming the user also has values for D and r from the data, T AS can then be derived from Z 2~TAS Dr.

Home range patterns and relations to allometry
Since the animal probability distribution reaches a steady state, it is possible to calculate both the size of the resulting home ranges and the degree to which they overlap. By using the 95% MCP method [15], the dimensionless radius of the home range, after dividing by the mean distance between CPs, is given by R 95% , implicitly defined by the following equation This allowed us to plot R 95% as various functions of k, one for each a (figure 4a). Each of these can be approximated by a sigmoidal function of K~Log 10 (k). Specifically, where Q&0:495{0:010a, n&0:278{0:048a, f&3:18 and g&3:139{0:025a (figure 4a).
For certain values of k and a, the value of R 95% is less than 0:5, meaning that gaps arise between adjacent territories. These socalled buffer zones have been observed between wolf (Canis lupus) territories [4] as a safe place for wolf prey, such as white-tailed deer (Odocoileus virginianus), to occupy.
The allometric predictions of [13] show that the fraction of exclusively used area E is approximately proportional to M {j where j&1=4 and M is the mass of a single animal. In our model E~(1{R 95% ) 2 =R 2 95% so allometric studies predict (1{R 95% )=R 95% !M {j=2 . By using the values of j fitted from the large data sets in [13], the value of R 95% can be estimated for an animal of given mass. Using the trend lines from the simulation plots in figure 1b and equation (10) allows k and a to be related to Z 2 , thus estimating how Z 2 scales with M, as shown in figure 4b.
In [13] the tendency for larger animals to have a lower proportion of exclusive area in their home ranges was explained intuitively, by noticing that they are less efficient than smaller animals in patrolling their territory to deter conspecifics. That is, the time it takes for a larger animal to get around its territory is greater than that of a smaller animal. In our model, this means the diffusive time, T D2 , increases with mass. Our results show that this ability to deter conspecifics is also driven by an additional factor: the active scent time. The ability to maintain exclusive area in fact arises from the ratio of T AS to T D2 . Figure 4b shows that a smaller animal's ability to maintain a higher proportion of exclusive space use arises from maintaining a higher ratio of T AS to T D2 , not just a lower diffusive time.

Comparison with previous approaches
Territoriality in animals with central place attraction has been studied previously in [4] using a reaction-diffusion formalism, which was developed further in [9]. Although both that model and the one presented here use conspecific avoidance mediated by scent marking as the mechanism of territory formation, the present Panels (e-g) compare the two distributions for the 2D system. The black contours show the deciles (i.e. 10%, 20%, 30% etc.) of the height of the probability distribution for the simulation system. The red contours show the same quantities for the analytic approximation. The values used were (e) k~0:011, a~0:80, (f) k~0:0014, a~4:0, (g) k~0:022, a~0:80 and (h) k~0:016, a~4:0. As we increase k or a, the effect of the adiabatic approximation becomes more apparent, since each red contour is further away from the respective black contour. This is due to the fluctuations of the territory border being more pronounced for higher k or a. doi:10.1371/journal.pone.0034033.g003 model is built from the individual-level interaction processes, whereas the reaction-diffusion model relies on a mean-field approximation for the scent mark response. Despite the very different natures of their construction and the resulting expressions, we compare the two models by examining the conditions under which they are numerically similar.
In the reaction-diffusion model, u( x x,t) and w( x x,t) are the dimensionless probability density functions for the left and right animals respectively. In addition, p( x x,t) and q( x x,t) denote the dimensionless densities of the scent of the left and right animals respectively. The dimensionless diffusion constant of each animal is given by d and the dimensionless advection coefficient controlling the strength of motion away from conspecific scent and towards the CP is c. The model also contains a parameter controlling the over-marking response rate: that is, the tendency for animals to scent-mark more having encountered foreign scent. However, since the animals in the model described in the present paper are counter-markers rather than over-markers [16], that is they mark next to conspecific scent but they do not increase marking rate as a response to scent, this parameter is set to 0. With these conditions, the reaction-diffusion system described in [9] has the following dimensionless steady state solution where 0ƒ x xƒ1, together with the probability conservation conditions Equation (12) is equation (6.11) in [9]. The dimensionless parameter b is a function of 5 dimensional parameters, b~c'l=(mLD), where L and D are the same values as used elsewhere in the present study, l is the scent marking rate for the individual or pack, m is the rate of scent-mark decay and c' is the strength of attraction towards the CP. The parameter c' is not the same as the drift velocity v from our model since it has units of space 3 =time rather than space=time. Indeed, the drift velocity at any point x x in the reaction-diffusion model is proportional to the strength of foreign scent at x x (see equations (4.5) and (4.6) in [9]), whereas in the model studied in the present paper the magnitude of the drift velocity is constant throughout space.
The way the rate of scent deposition is modelled also differs between the two approaches. In the reaction-diffusion model, the rate is independent of the magnitude of the animal's diffusion constant. The biological implication being that as the animal's speed increases, consecutive scent marks will be deposited further apart. In our model, the scent marks are deposited every time the animal has moved a distance a (the lattice spacing), regardless of its speed. The reason for our choice is that it is advantageous for animals to ensure that they deposit territorial messages at regularly spaced intervals so that they leave no gaps in the territory boundaries, which might allow conspecifics to intrude.
Scent decay is also modelled in different ways in the two models. In the reaction-diffusion model the scent decays exponentially, whereas we assume scent is ignored after a fixed period of time (T AS ). Whilst exponential decay of scent makes sense regarding the decay of the chemicals that produce the odour, a conspecific may ignore a scent mark it can still smell, if the odour suggests that the mark is old and the territory is no longer being defended. For example, such behaviour has been reported for brown hyaenas (Hyaena brunnea), whose scent marks may still be detectable by conspecifics over a month later, but who tend to ignore scent that is more than about four days old [17].
Making numerical comparisons of our model with the reactiondiffusion model required a further reduction of our 1D analytic model, since the 1D reaction-diffusion model only represents animal movement in the right-hand (left-hand) half of the left-hand (right-hand) territory. Focussing on the left-hand territory, this . Home ranges and allometry. Panel (a) shows how the radius R 95% of the normalised (by dividing by the mean distance between CPs) 95% minimum convex polygon home range depends on k and a in the 2D analytic model. The various shapes (circles, squares, crosses etc.) show the exact values and the solid lines show the least-squares best-fit sigmoidal curves. Notice that whenever R 95% v0:5, a buffer zone appears between adjacent territories. The proportion of exclusive area E~(1{R 95% ) 2 =R 2 95% scales with mass [13] so this value is plotted in panel (b) against the dimensionless parameter Z 2~TAS Dr for various a. Again, solid lines are derived from the best-fit sigmoidal curves whilst the points denoted by various shapes show exact values. doi:10.1371/journal.pone.0034033.g004 required us to simplify our model by fixing g g 1 ( L L 1 )~d( L L 1 ) where d(z) is the Dirac delta function. The resulting marginal distribution for the position of the animal in dimensionless coordinates is This expression is compared with the distribution u( x x) from the reaction-diffusion model, whereas w( x x) is compared with M R (1{ x x). To find the best fit, the square of the difference between the curves of ln½u( x x) and ln½M R (1{ x x) is minimised ( figure 5).
Though the two models are qualitatively very different, if a and k are both very small, it is possible to find a value of b that fits closely (figure 5a). However, if either a or k are increased, even the best fit value of b gives a qualitatively different curve. Conversely, for lower values of b, the best fit curve to the model studied here becomes increasingly different to the curve from the reactiondiffusion model (figure 5b).
To explain the similarities in these parameter regimes, the limit case where the scent marks never decay is examined, so that T AS ?? and k?0. If in addition a?0, the marginal distribution x xw1=2. The analogous limit in the reaction-diffusion model is m?0 so that b??. In this limit case, u( x x) and w( x x) are step functions. By taking the limit numerically as b??, one observes that u( x x)~2 if x xƒ1=2 and u( x x)~0 if x xw1=2 so that u( x x) and M R ( x x) coincide. Similarly, w( x x) and M R (1{ x x) coincide in this limit.
Whilst our model has two parameters, as opposed to one in the reaction-diffusion model, it is possible to collapse our model to one parameter by formally taking the limit a?0 in equation (14), giving the following expression where Ci (z)~{ Ð ? 0 dt cos(t)=t is the cosine integral. This is precisely the limit where the reaction-diffusion model tends to agree best with ours. Plots of equation (15) can be found in the insets (solid lines) of figure (5) for those cases where a~0.

Discussion
A central place foraging model with scent-mediated conspecific avoidance has been constructed where the mechanisms of both the animal movement and the interactions are defined at the level of the individual. Territories naturally arise with slowly fluctuating borders whose probability distribution tends towards a steady state. Stable home range patterns emerge, easily enabling us to quantify the home range size and overlap as a function of the underlying individual-level movement and interaction mechanisms. Whilst this is not the first mathematical model of territoriality in central place foragers, nor is it the first where the movements are built mechanistically from individual-level processes, [9] it is the first where the conspecific avoidance mechanism is built from interactions between individual agents. Though certain predictive inferences have been made using previous approaches, for example regarding what happens when a territory dissolves [18], ours is the first where predictive inferences can be made about the mechanisms of territorial interaction events, in particular the active scent time, from the patterns of animal space-use. Although deterministic reaction-diffusion equations are in general viable approximations to represent spatio-temporal stochastic processes, they are not well suited to model systems in which the individual components are present in low concentrations (e.g. [19,20]). This is precisely the situation of decaying scent marks in our model, which are ignored by conspecifics beyond the time T AS . When that happens, the probability density associated with that scent location is identically zero. In other words, in a scent-mediated interaction process the extinction probability for the scent is non-negligible. The system is intrinsically stochastic, and deterministic approximations, where the occupation probability is coupled to a scent mark profile, may not cope with the discrete nature of the interaction events. A reaction-diffusion formalism may thus provide results that are in complete disagreement with the stochastic description (see e.g. [21,22] in the spatial ecology literature). The particular reaction-diffusion model studied in [9] has been shown here to give very different results to our model, away from the limiting case where scent marks never decay. The similarity in this limit does not come as a surprise, since in this case the scent is never present in low concentrations.
Away from this limit, the choice of model that is most appropriate for a particular data set would depend on both the species involved and the questions to be answered. If one is interested in quantifying both scent marking mechanisms and animal movement processes, a drawback of the reaction-diffusion model is that the dimensionless parameter b governing the animal space use distribution is a product of 5 (dimensional) parameters, including both the strength of central place attraction and details of the scent marking process. This makes it very difficult, if not impossible, to quantify the scent marking mechanism purely by fitting data to the animal probability density function. On the other hand, the present study gives a clear programme for inferring both the strength of the central place attraction and the active scent time by fitting data on animal space use.
This programme for inferring T AS from animal location data was not developed in previous agent-based studies, since the probability distribution of the animal positions never reaches a steady state [6]. In such systems, it is necessary to pick a biologically meaningful time-window over which to measure the extent of home range overlap and thus infer the nature of the border movement and, in turn, the active scent time. This procedure is required for analysing certain animal populations, such as urban red foxes, whose territories, in certain circumstances, may not reach a steady state. However, if a steady state is reached, as shown by a saturating animal MSD, then some aspect of the underlying movement process must be keeping the animal from continually spreading out across the terrain. Such stable home ranges have been reported in a number of species (see e.g. [23]) from wolves (Canis lupus) and coyotes (Canis latrans) [9] to hispid cotton rats (Sigmodon hispidus) [24], cane mice (Zygodontomys brevicauda) [26] and Baird's tapirs (Tapirus bairdii) [25]. One possible mechanism for ensuring this stability is central place attraction, studied here. It may also be possible that some form of memory mechanism keeps the animal in familiar environments and thus causes the probability distribution to saturate [27,28].
Our study also gives insights into the mechanisms behind the allometric scaling of exclusive space use. Previous studies had interpreted the observed scaling laws as a consequence of a greater ability for smaller animals to cover their territory regularly, compared to larger animals. However, by quantifying how the scaling arises from the ratio between the active scent time and the territory coverage time, we have shown that the longevity of territorial messages is also a fundamental quantity. Future studies on allometric scaling of space use should also take into account this mechanism of interaction.

The stochastic simulation model
The 1D simulations consisted of 2 animals on a finite lattice with periodic boundary conditions. The central places (CPs) for each animal were uniformly distributed at a distance L~aN apart, where a is the lattice spacing and N a positive integer. In 2D, 30 animals in a rectangular terrain with periodic boundary conditions were simulated. The CPs were placed at the centroids of a hexagonal lattice, modelling the fact that animal territories tend to be roughly hexagonal in shape [29]. Adjacent CPs were separated by a distance of L. The simulated animals deposited scent at every lattice site they visit, which remained for a time T AS , the active scent time. Animals were unable to visit sites that contained scent of another animal. Besides that constraint, at each step an animal moved to an adjacent site at random but its movement was biased towards the CP. In 1D, this meant that there was a probability of pw1=2 of moving towards the CP and 1{p of moving away. In 2D, the movement probabilities were as follows where (m,n) is the position of the animal and (m c ,n c ) the position of the CP. These probabilities were chosen so that in the continuum limit, they reduce to the form that gives the correct localising tendency in the Holgate-Okubo model (see the section 'Reduced analytic model in 2D'). In particular, they are independent of the distance the animal is away from the den site. This can be shown by replacing m{m c by C(m{m c ) and n{n c by C(n{n c ) in equations (16), for some non-zero constant C, and noticing that all the C-values cancel.
Simulations were run until the MSD of the border had reached a saturation value. Each 1D simulation result was an average of 1,000 simulation runs. In 2D, it was only necessary to average over 100 runs owing to the fact that 15 times as many animals were simulated per run. The simulations were coded in C and compiled on Windows XP OS. To obtain a single saturation MSD value for the 2D simulations took an average of 4 hours 40 minutes CPU time using a 3.0 GHz processor in a 2.96 GB RAM desktop computer.

The reduced analytical model in 1D
To understand the nature of the animal's movement within its territory borders, we considered a simplified analytic model that uses an adiabatic approximation similar to [11] because the animal moves at a much faster rate than the borders. This meant that the joint probability distribution of the animal and the borders could be decomposed as P 1D (x,L 1 ,L 2 ,t)&Q(L 1 ,L 2 ,t)W (x,tDL 1 ,L 2 ) where Q(L 1 ,L 2 ,t) is the probability distribution of the borders to be at positions L 1 and L 2 at time t, and W (x,tDL 1 ,L 2 ) is the probability distribution of an animal to be at position x at time t when constrained to move between the borders at L 1 and L 2 .
Following [11], the borders were modelled using a Fokker-Planck formalism, with time-dependent diffusion constant modelling the subdiffusive nature of the border movement, and quadratic potentials modelling the spring forces (figure 2). Since the CP at B separates L 1 from L 2 , we write Q(L 1 ,L 2 ,t)~Q 1 (L 1 ,t)Q 2 (L 2 ,t) where Q 1 (L 1 ,t) and Q 2 (L 2 ,t) are the probability distributions of L 1 and L 2 respectively. These are governed by the following equations where A is the position of the CP to the left of L 1 , B is the position of the CP between L 1 and L 2 , C is the position of the CP to the right of L 2 , KQ(t) is the time-dependent diffusion constant and cq(t)½L 1 {(AzB)=2 2 =4 (resp. cQ(t)½L 2 {(BzC)=2 2 =4) is the quadratic potential for each spring connected to L 1 (L 2 ). It ensures that the border L 1 (L 2 ) fluctuates around an average position of ½AzB=2 (½BzC=2). Notice that there are two springs connected to L 1 (L 2 ), so that the total resulting potential is cQ(t)½L 1 {(AzB)=2 2 =2 (resp. cQ(t)½L 2 {(BzC)=2 2 =2). As usual for Fokker-Planck equations (see e.g. [31]), this potential appears in equation (17) (resp. 18) after having been differentiated with respect to ). The boundary conditions in equations (17) and (18) ensure that the borders cannot cross over the CPs, since each CP must remain in its territory.
K and c can be measured directly from the simulation model (see e.g. [11]). However, in the steady state solutions (equations 2 and 3), K and c collapse to a single parameter k~K=(cL 2 ). The k parameter can be derived by first measuring the boundary's saturation MSD from the simulations, and then using equation (10).
Equations (17) and (18) can be solved using the method of characteristics [30]. The general solution to (17) is where is defined so that G'(t)~cQ(t) and L 1,0 is the initial value for L 1 at time t~0. By using the method of images [32] to take account of the boundary condition and assuming, for simplicity, that L 1,0~( AzB)=2, we arrive at the following solution Similarly, By making use of the Poisson summation formula [32], equations (20) and (21) can be re-written as follows Since the territories move as tagged objects in a single file diffusion process [6], we have Ð t 0 dsQ(s)*t 1=2 in the 1D system [8].
Therefore the limit as t?? of b(t) is b(t~?)~4K=c. Taking this limit in equations (20) and (21) gives steady state solutions. Furthermore, by setting B~0, B{A~C{A~L, and using dimensionless variables L L k~Lk =L, x x~x=L, Q Q k ( L L k )L Q k (L k ,t~?), g g k ( L L k )~Lg k (L k ,t~?) for k~1,2 and k~K=(cL 2 ), we obtained expressions (2) and (3), displayed earlier in the Results section.
To calculate the animal probability distribution W (x,tDL 1 ,L 2 ), we began by finding the continuous-space limit of the simulation model in the case where the animals and their CPs are infinitely far apart so that they never interact. This corresponds to W (x,tDL 1~? ,L 2~? ), written as W ? (x,t) to ease notation.
The master equation for an animal in this limiting case is where U(n,t) is the probability of the animal being at position n at time t, n c is the position of the CP, F is the jump rate between adjacent lattice sites, and sgn(z)~1 (sgn(z)~0, sgn(z)~{1) if zw0 (z~0, zv0). This can be written as The continuum limit of (25) can be found by taking the limits as a?0, F ??, n??, n c ?? and p? 1 2 such that D~a 2 F , x~an, x c~a n c and v~2aF(2p{1) [33,34]. Physically, D is the diffusion constant of the animal, v the drift velocity towards the CP, x the position of the animal and x c the position of the CP. This procedure leads to the 1D Holgate-Okubo localising tendency model [35,36] LW ? (x,t) Lt~D wherex x~{1, 0 or 1 if xwx c , x~x c or xvx c respectively. This has a non-trivial steady state solution [9], proportional to exp({vDx{x c D=D). Since the animal is constrained to move between the borders at L 1 and L 2 , the probability distribution must be zero for xvL 1 and xwL 2 . As the solution is a steady state, the flux across L 1 and L 2 is automatically zero so it suffices to ensure that the integral of the probability distribution between L 1 and L 2 is equal to 1.
Using dimensionless variables a~vL=D, t~?) and setting x c~0 for simplicity, we obtain equation (4) from the results section.

The reduced analytical model in 2D
In 2D we modelled each territory as a circle with fluctuating radius and the CP at the centre of the circle, assumed to be at the origin for simplicity. As in the 1D scenario, we used an adiabatic approximation so that P 2D (r,h,s,t)&Q(s,t)W (r,h,tDs), where P 2D (r,h,s,t) is the joint probability distribution of the animal to be at position (r,h) in polar coordinates at time t and the territory radius to be s. Q(s,t) is the probability of the territory radius to be s at time t and W (r,h,tDs) is the probability of the animal to be at position (r,h) at time t in a territory of fixed radius s.
Similar to the 1D scenario, Q(s,t) was modelled using a Fokker-Planck formalism with the radius fluctuating around an average value of L=2, where L is the distance between adjacent CPs. As such, it can be calculated using the methods of the previous subsection to be Q(s,t)~½H(s){H(s{L)g(s,t) As the territories are tagged particles in a 2D exclusion process [6], we have Ð t 0 dsQ(s)*t= ln (t) [7]. Taking the limit t?? in equation (28) gives a steady state solution. This gives rise to the expression (6) from the main section by using dimensionless variables s s~s=L, Q Q( s s)~LQ(s,t~?), g g( s s)~Lg(s,t~?). Following our methods in 1D, we calculated W (x,tDL 1 ,L 2 ) by first taking the continuum limit of the master equation governing the movement of an animal unconstrained by other territories (i.e. L 1~L2~? ). This master equation is where U m,n (t) is the probability of the animal being at position (m,n) at time t and (m c ,n c ) is the position of the CP. To find the continuum limit, this is re-written as follows wherex x is the unit vector pointing from the animal at (x,y) towards the CP at (x c ,y c ), or the zero vector if (x,y)~(x c ,y c ), and W ? (r,h,t) is the probability distribution W (r,h,tDL 1 ,L 2 ) in the limit as L 1 ,L 2 ?? where there is no interaction with other animals. As in 1D, (31) has a non-trivial steady state solution [9], which is proportional to E 1 (vr=D). The boundary condition ensuring that rƒs, so that the animal is within its territory, is imposed by normalising the steady state solution so that the integral over the circle, of radius r centred at 0, is equal to 1.
By using dimensionless variables a~vL=D, r r~r=L, h h( r r,hD s s)~L 2 h(r,h,t~?Ds), we obtained equation (7) from the results section.