Front. Phys.Frontiers in PhysicsFront. Phys.2296-424XFrontiers Media S.A.10.3389/fphy.2019.00129PhysicsOriginal ResearchFrom Micro-to-Macro: How the Movement Statistics of Individual Walkers Affect the Formation of Segregated Territories in the Territorial Random Walk ModelSarvaharmanSeeralan1Heiblum RoblesAlexandro2GiuggioliLuca13*1Department of Engineering Mathematics, University of Bristol, Bristol, United Kingdom2Foro Consultivo Científico y Tecnológico, Mexico City, Mexico3Bristol Centre for Complexity Sciences, University of Bristol, Bristol, United Kingdom
Edited by: Ralf Metzler, University of Potsdam, Germany
Reviewed by: Diego R. Amancio, University of São Paulo, Brazil; Muktish Acharyya, Presidency University, India
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Animal territoriality is a widespread phenomena in many vertebrate species. In mammals it is often associated with territorial marking with which individuals make their presence conspicuous to others by leaving trace of their passage, often in the form of deposited scent. A simple interaction mechanism consisting of retreating upon the encounter of a foreign scent is sufficient to observe the emergence of territorial patterns at the population level. With the introduction of the so-called territorial random walk model this local avoidance mechanism coupled with a simple diffusive movement of the individuals has been shown to generate long-lasting patterns of segregation at much larger spatial scales. To shed further light on the micro-to-macro connection of this collective movement model we study how the movement statistics of the individuals affect the formation of the segregated scented territories. We represent individual animals as correlated random walkers and we analyse the spatial ordering of the population as a function of the length of time a scent mark remains active after deposition and as a function of the degree of correlation of the movement steps. For low and intermediate correlation strength we find that territories undergo a liquid-hexatic-solid transition as active scent time is increased. Increased spatial order also appears by increasing the correlation strength but only if well away from the ballistic limit. We ascribe this non-monotonic dependence to the coverage efficiency of the individual walkers mainly controlled by the correlation and the mobility of the territories mainly controlled by the active scent time.
territorial random walkcorrelated random walkKTHNY meltingtopological defectsanimal territorialityEP/I013717/1Engineering and Physical Sciences Research Council10.13039/5011000002661. Introduction and Background
In biology it is rather common to find a system in which the underlying movement of its constituent parts is not diffusive, often owed to the out-of-equilibrium nature of the processes involved. Examples of anomalous transport can be found at all scales: from the active and passive microrheology of various nanoparticles inside the molecularly crowded environment of cells [1] and the two stage diffusion of macromolecules on cell membranes [2, 3] to the superdiffusive displacement of epithelial cells [4], and all the way to the anomalous dynamics of whole organism while foraging randomly [5] or during memory biased search [6, 7].
Among the cited biological examples characterized by anomalous diffusion we have focused for this special issue on the movement of whole organisms and on the collective effects of many such organisms spatially excluding one another. Our interest here lies in understanding how the local statistical features of the movement of individual animals affect the emergence of the collective patterns at the population level. We choose a movement statistics that becomes random at long time scales, but retains a degree of persistence at shorter time scales, the so-called correlated random walker [8], a paradigmatic movement model in animal ecology [9, 10]. Among the forms of spatial exclusion or avoidance we are particularly interested in a very common behavior of vertebrate populations in 2D: the subdivision of the terrain into spatially segregated regions [11]. When these regions are of exclusive ownership of a single individual or a single family unit they are called territories [12, 13].
The purposes of territoriality change from species to species or even throughout the year [14]. These include, roosting, mating, nesting, and harboring resources. Accordingly the interaction mechanisms that animals rely upon to form and maintain these territories are quite rich and depend on the type of signals that animals exchange. Broadly speaking one may distinguish them between direct and indirect, that is whether the time scale for the signals to travel from the emitter to the receiver is short or long relative to the time scale for the emitter to move. Examples of direct interactions, often employed by birds [15], is the use of visual displays and audio calls. In this case the signal of the emitter is detected by the receiver nearly instantaneously. Examples of indirect interaction, used by a large number of mammals [16], is the use of olfactory cues. In this case the scent that an animal deposits is nearly static. For the period over which the scent remains detectable, which could be a long time after deposition, the environment retains the memory of the passage of the emitter. Any individual that comes in close proximity to the deposited scent becomes a receiver and acquires information about the (past) presence of the emitter.
In this study we analyse a specific case of olfactory-based territorial formation, often called conspecific avoidance, whereby animals mark the terrain wherever they go, and other individuals passing by respond to these olfactory cues by retreating from the region or area where foreign scent was encountered. This indirect animal interaction, that occurs through the modification of the environment, is called stigmergy [17]. It is a common form of interaction in eusocial insects [18], but it has been shown to occur also in territorial animals [19].
The mathematical study of scent-marked territorial patterns has a relatively long history dating back to the early '90s when the first reaction-diffusion model representing a pair of animals avoiding each other scent was formulated [20]. This model coupled the occupation probability of two Brownian walkers tethered to their respective den or burrow and their scent profiles. It was later generalized to include the effects of landscape heterogeneity and animal movement responses and applied to movement data on wolves and coyotes [21].
In subsequent modeling studies on the formation of scented territories a different approach was followed [22]. That approach becomes necessary when one aims to account for the sharp spatial dependence of the interaction. In these cases the field nature of the interaction potentials (or forces), which tacitly assumes that they are defined at every point in space, may not be adequate. It is more convenient to account for the interaction dynamics through localized walls or spatial partitions representing the deposited scent [23]. By doing so one remains faithful to the biology of scent-marking species for which individuals react to the encounter of foreign scent only if it is informative, that is only if the deposition occurred before a certain time in the past [see e.g., Alberts [24]].
As this approach requires tracking the movement and interaction of the entire population to determine the time-dependent position of each animal and the age of the scent deposited, it was formulated as an individual based model [22]. This increased complexity, however, allows to consider animals that do not have a den or a burrow as well as to study the emergence of territorial patterns as a collective phenomena rather than a two-body problem between neighboring individuals.
The collective movement model of territory formation, termed the territorial random walk (TRW) model [23], lends itself naturally to questions on the nature of the emerging patterns, e.g., if macroscopic order and/or disorder phases appear and how the microscopy of the movement and interaction rules influence the emergence of the macroscopic patterns. Along these lines a very recent investigation on the presence of order-disorder phase transitions in the TRW model as a function of the scent decay time and the population density has been conducted [25]. In that study it has been shown that the emerging territorial patterns display a solid to liquid melting scenario analogous to the one supported by the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory of melting [26–28], with the appearance of an intermediate partially ordered phase, called hexatic, due to the importance of geometrical arrangements of the first six neighbors.
Here we extend that analysis by considering a more realistic movement model for the individual animals. We modify the movement statistics of the animal by representing them as correlated random walkers and we ask how the “microscopic” movement of the individuals alters the collective dynamics of the system and the order-disorder scenario.
The paper is organized as follows. In section 2 we introduce the model and we present details of the stochastic simulations and the thermalization of the system. The analysis has two parts to it. At the “microscopic” level, that is at the level of the individuals, we study the influence of the persistence of the walk on the variance of the occupation probability, that is the mean square displacement (MSD). This is presented in section 3 together with an analysis of the spatial coverage of each walker in its own territory. As we scale up to the level of the territories, we analyse the effects of the animal movement statistics on the appearance of ordered phases in the system. This is dealt with in section 4. And finally section 5 presents concluding remarks.
2. The Correlated Territorial Random Walker Model
The correlated territorial random walker (cTRW) model, an extension of the TRW model, is a lattice-based collective movement model where individuals move with some degree of persistence, that is the movement direction at each time step depends on the previous step direction. The model can be run in any dimension and in discrete or continuous time. For computational efficiency we have run it in discrete time and we have only analyzed the 2D case. To study the structure of the emerging territories it is important not to introduce spatial frustration effects and we have thus used a triangular lattice with periodic boundary conditions.
In the absence of interactions a walker's steps are correlated. To model the persistent random walk on a triangular lattice we start from a continuous turning angle variable θ ∈ (−π, π], and we create six bins of width π3 for each of the turns a walker can make. The angle θ is then drawn randomly from a wrapped Cauchy distribution, C(θ)=(1-ρ2)[2π(1+ρ2-2ρcos(θ))]-1 where −π ≤ θ < π and ρ is the persistence parameter, or the mean cosine of the distribution. It indicates the tendency of a walker to continue to move in the same direction from where it came from at the previous step. In the limit ρ → 0, the distribution reduces to C(θ)=12π i.e., the turning angles are uniformly distributed and the movement of the walker becomes random. The opposite limit corresponds to the ballistic case C(θ)→δ(θ), where δ(θ) is the Dirac delta function as ρ → 1.
Avoidance between individuals occurs as follows. As animals deposit marks wherever they go, at each time step any individual may encounter one of its own marks or may step on a lattice site with one or multiple foreign marks. When an animal encounters its own marks, no interaction occurs. On the other hand upon encountering a foreign mark within a time period shorter or equal to the so called active scent time TA from when it was deposited, an interaction occurs. In order to avoid intruding further into neighboring territories. An interaction consists of one of three types of retreat at the following time step. These are in order of preference, first toward a neighboring site with no scent marks, second toward a neighboring site with the animal's own scent and lastly in the rare case where neither of the previous options are available the animal moves randomly to one of the six neighboring sites. At the subsequent step the animal reorients itself with a new direction chosen at random between the six potential choices.
As a result of the above interactions, foreign scents serve as spatial barriers effectively making individuals avoid entering further into others' territories. At any given time t a territory is represented by all the sites that contain an active scent of a specific animal. Alternatively it represents the sites an individual has visited within the interval (t−TA, t). Except at the boundaries, where a given territory may overlap with one or more neighboring ones, the avoidance of region recently visited by other individuals create spatially segregated areas of exclusive ownership for each animal.
The spatio-temporal dynamics of the territories can be fast and slow depending on the choice of parameters, which are, respectively, the inverse of the population density or specific volume ν = L2/N where L2 is the domain size, the active scent time TA and the persistence parameter ρ. With small TA and large ν one observes highly mobile and morphing territories as exclusion interactions are not very frequent. The opposite happens with large TA and small ν as interactions occur more frequently.
The dependence on the initial condition of the system is dealt with by running simulations and waiting for the system to thermalize before making any spatial and/or temporal measurement on appropriate quantities of the system. The thermalization is performed as follows. Two different initial conditions are used. The first is initialized with the scent profile tessellating the available space, i.e., the territory of each walker is a perfect hexagon with all boundary lattice sites overlapping with the neighboring territories. The second is initialized with randomly distributed walkers with no scent profile. The two are left to run until the standard deviation of the territory size across the entire population converges to the same value in the two cases. When that happens the system is deemed to be thermalized.
Rather than exploring the dynamics for different population density we have selected an intermediate value of ν, namely ν = 48 for multiple reasons. On one hand the rich (slow and fast) dynamics can be attained within a limited range of TA values so that disordered and ordered phases of the system could be observed and studied as a function of ρ. For smaller population densities, as interactions become more infrequent, it becomes harder to observe spatial ordering and we have thus not explored the regime with larger ν. On the other hand smaller values of ν makes the computation very expensive and makes it harder to reach thermalization except for very small values of TA.
The specific range and resolution of the parameter space ρ − TA with ν = 48 has been as follows. The active scent time TA has been chosen in the range between 89 and 1,424 with a resolution of 89, except for specific cases detailed in the figure captions. For the correlation parameter ρ we have used the values 0, 0.15, 0.35, 0.55, 0.75, and 0.95. To ensure that the territories, with the appropriate choice of TA, can tessellate with equal hexagons the 2D domain, the linear domain size (with periodic boundary conditions) is set to L = 1200 with N = 3 × 104, which corresponds to ν = 48, although we have also used smaller domains e.g., in section 3.
3. Partial “Caging” and Spatial Coverage of the Individual Walkers
The spatial exclusion of the scented territories is a complex collective phenomenon whereby an animal remains confined in certain region of space depending on when and where foreign marks as well as its own marks have been deposited. But this confinement is only partial since a mark has a lifetime and disappears unless an animal remarks that same location within TA steps. Animals may thus be trapped for some time before escaping. It is possible to observe such partial “caging” by plotting the MSD as a function of time. We do so in Figure 1 where we draw the normalized MSD for TA = 712, 1,068, and 1,424 for various values of ρ by averaging over all individuals in one simulation (after the system has thermalized). The MSD plots for TA = 1 and TA = 89 are additionally shown in Figure A1. We term these snapshot averages to distinguish them from ensemble averages when multiple copies of the system are run.
MSD Δ(t) of a walker in a system with ν = 48 and L = 120, where each curve is a snapshot average. The active scent times (TA) for panels (A–C) are 712, 1,068, and 1,424, respectively. At short timescale the walkers tend to exhibit superdiffusion indicated by the Δ(t)∝tα relationship where α > 1. For intermediate timescales the dynamics of the walkers are governed by the mobility of their territories. High values of TA yield slowly moving territories and increase in caging caused by neighbors. This caging saturation is highly sensitive to the correlation of the walkers (see Figure A2) while the time spent at the caging saturation before increasing again is related to the coverage efficiency of the walker.
In moving from panel (a) to (b) and (c) in the figure, that is by increasing TA, one clearly observes a reduced growth of the curves at intermediate times for each value of ρ. Since a long-time saturating MSD is expected whenever a walker roams inside a finite domain [29], the flattening of the MSD at intermediate times is a manifestation of increased caging from the neighbors. As marks remain active for longer they make the territories progressively less mobile and the MSD curves start increasing again at later times.
A comparison of different values of ρ within the same panel shows that this caging effect becomes stronger—curves remain flat for longer time—for larger ρ except when one approaches the ballistic regime. This is clearly visible in panel (c) when TA = 1, 424 where this non-monotonicity as a function of ρ can be seen more clearly. There also appears some non-monotonic dependence of the intermediate-time saturation value of the MSD, but this non-monotonicity is only apparent. As we have verified and shown in Figure A2 in the appendix, the intermediate saturation of the MSD is highly erratic as it is very sensitive to the actual shape of the occupation probability for each animal. We consider the time it takes the MSD curve to start increasing again the more significant feature, and to understand the non-monotonicity as a function of ρ we look in detail at the walker spatial coverage.
Considering initially a fixed finite domain, for a walk to improve spatial coverage it needs to reduce spatial oversampling, which clearly happens by increasing ρ. However, in a confined domain, when the degree of correlation is too large, the walker would retrace back its steps thus increasing again the sampling of lattice sites already visited. It is well known in fact that the mean coverage time of an independent walker in a finite domain can be minimized for intermediate values of walk persistence [30, 31].
To test whether this understanding suffices to explain the non-monotonic ρ dependence of the caging effect, we study the walkers' spatial coverage in the cTRW model. As territories are mobile and change shape, we cannot compare values of coverage time for each walker. However, it is possible to obtain information about the efficiency of spatial coverage by plotting the average number of sites an individual covers in a time TA as a function of ρ. The outcome is shown in Figure 2 where we have plotted Cp, the fraction of sites covered in a time TA, for different values of ρ and TA. From the figure it is evident that the partial coverage of the terrain is maximal for intermediate values of ρ and for TA sufficiently large. As TA affects the size of the confining domain, the correlation parameter that maximizes coverage depends on how long marks remain on the terrain.
The fraction of the territory defined at time t covered within TA time steps by walkers as a function of the correlation parameter ρ, and the active scent time TA. For each pair of ρ and TA a snapshot average is taken on a thermalized system with L = 120 and ν = 48. For very low TA the territories are highly volatile, i.e., within a TA time window the territories are no longer the ones defined at time t which results in low Cp independent of ρ. For higher values of TA the territories are less mobile and we see the maximization of territory coverage for an intermediate values of ρ.
4. Spatial Ordering of the Emerging Territories
The interaction dynamics of the individual walkers have a long lasting effects on the spatial structure that emerges in the population, as was shown in the case of the simpler TRW model [25]. The spatial ordering of the emerging territories is dictated by the “caging” effect whose size corresponds to the size of the territory. To analyse the spatial ordering as a function of ρ and TA we look at the system with a coarser temporal and spatial resolution because no order is present at distances smaller than the size of the territory. We neglect the dynamics of the walkers and look at the territory centroids over a TA time resolution. As territorial centroids are calculated only by considering the mean position of all lattice sites with active scent every TA steps, we obtain a coarse-grained (continuous-space) mesoscopic time-dependent description of the cTRW model.
4.1. The Pair Correlation Function
For a coarse measure of spatial ordering, we plot in Figure 3 the pair correlation function [32], g(r)=ν〈∑i≠0δ(|r→-ri→|)〉, where δ(z) is the Dirac delta function, as a function of the centroid distance r=|r→| for different ρ and TA values. In particle systems the pair correlation function gives the likelihood of finding a neighboring particle as a function of distance relative to the ideal gas case. In an ideal gas, g(r) equals 1 for any r as individual particles do not interact with each other and are equally likely to be anywhere in space. A value of zero indicates instead the impossibility of having another particle at that distance in space, e.g., for hard-core interaction of particles of radius σ, g(r) remains zero for all distances up to r = 2σ. Except for extremely small values of TA, the emerging territories in our model have an effective hard core interaction as can be observed for any of the curves in Figure 3. Those curves remain close to zero up to 2σ where σ=aν/2, with a the lattice spacing, represents the distance between two evenly spaced centroids when hexagonal territories tessellate the entire domain (see explicit calculation for σ in Heiblum Robles and Giuggioli [25]). A rise from zero to the value one with no evident oscillations would be indicative of a gas with purely hard-core particles. On the other hand a rise beyond the value one with subsequent damped oscillations toward one represents the arrangement of other particles into “shells” of neighbors with the decay pointing to the radial distance over which such spatial ordering persists. Two phases may have qualitatively these characteristics: a liquid but also a hexatic phase. While it is not possible to distinguish a liquid from a hexatic phase from a pair correlation plot, it is possible to determine the appearance of a solid phase because the shells of neighbors are not arranged anymore radially. As a solid possesses hexagonally arranged particles, the radial symmetry is not present anymore and the smoothness of peaks and troughs in the pair correlation function is lost.
The pair correlation function g(r) for a system with L = 1, 200, ν = 48. Moving from panel (A–C), the TA values are 712, 1,068, and 1,424, respectively. Within each panel, the g(r) function has been shifted upwards from zero for visualization purposes. For the case ρ = 0 we have used the results from a previous study and as such the TA values are not precisely the same as the one used for the cases when ρ≠0; they are 707, 1,070, and 1,409 for panels (A–C), respectively. These differences do not take away from the qualitative differences presented here.
A feature that can be evinced by looking at Figure 3 is that the solid phase for sufficiently large TA appears for intermediate values of persistence ρ. This is particularly evident in panel (b) where the system has radial symmetry at ρ = 0, but then loses it as it approaches ρ = 0.3, and then regains it beyond ρ = 0.55. In panel (c) although the system is already in a solid phase at ρ = 0, the shape of the various g(r) plots also points to a progressive loss, even though only slightly, of radial symmetry as ρ is increased to 0.55, whereas radial symmetry is present when ρ = 0.75. This interesting dependence of the appearance of a solid phase depending both on the value of ρ and TA matches qualitatively the region in parameter space where the coverage efficiency of an individual is maximized, displayed earlier in Figure 2.
The partial coverage analysis in conjunction with the pair correlation plots indicate that an animal, for a given size of its own territory, that is for a given value of the active scent time TA, may select the most appropriate correlation statistics to be able to remark the majority of the scented terrain in TA steps. When that occurs, it implies that the neighbors are kept outside of the terrain that the animal defends. This in turns makes the territories less mobile, reducing the chance of having territory shape far from hexagonal and with large variability in sizes, and the entire population possess a bigger spatial order.
While Figures 2, 3 brings good evidence on how the microscopy of the movement statistics affect the emerging spatial order of the territories for small and intermediate ρ values, it does not seem a viable explanation for large values of ρ. From Figure 2 one would expect in fact a very high order also for ρ = 0.95 and high TA, but that does not appear to be the case in Figure 3C. However, the system may still possess a great degree of order if it were in an hexatic phase. To determine if that is the case, we first try to use appropriate order parameters that should help us to map out more precisely when the system is in a solid phase as well as to distinguish between a liquid and a hexatic phase as a function of ρ and TA.
4.2. Order Parameters
In Heiblum Robles and Giuggioli [25] some of the present authors have brought support for a continuous solid-hexatic-liquid transition akin to the KTHNY theory of melting as a function of TA in the TRW model. Starting from a very high TA in a perfect crystal configuration, by decreasing TA below a certain value makes the system loose translational order even though it retains orientational order and no phase coexistence has been observed, the transition being continuous. With a further decrease of TA the system looses also the orientational order to become liquid. We have extended these results by constructing the orientational and translational order parameters as a function both of ρ and TA.
We use the so-called local bond-orientational parameter [32] ψ6(ℓ)=1Nℓ∑jNℓexp(i6θℓ,j), where θℓ, j is the angle that particle ℓ, located at r→ℓ, makes with the j-th neighbor relative to a reference axis, and the j-summation is over the Nℓ neighbors, the latter obtained through the Voronoi construction [33]. In a perfect crystal, Nℓ = 6 and θℓ,j=π6, hence ψ6(ℓ) = 1. In Figure 4 we visualize qualitatively the changes in orientational order of the system by plotting the argument of the orientational order parameter, arg(ψ6), for each territory centroid. The orientational ordering of the system, that is when ψ6 → 1, corresponds to when arg(ψ6) → 0. From the various subplots one clearly notices that for a given ρ, if one increases TA, the hexatic phase does not always appear. It is also evident that for TA < 800 the system cannot reach the hexatic phase independently of the correlation parameter. On the other hand for intermediate values of TA, e.g., for a given TA in the range 800–1,100, the system is in a liquid state but then becomes eventually hexatic with sufficiently large ρ, and then liquid again with further increase in ρ. These findings support the pair correlation analysis depicted in Figure 3 and, in addition, they might also give a better idea of whether the system is in an hexatic or a liquid phase, which was not possible by looking only at g(r).
Change in the local bond orientational parameter ψ6(ℓ) for each territory centroid located at r→ℓ as a function of ρ and TA. The color of the centroid is determined by arg(ψ6(ℓ)), while the size of each centroid is proportional to the size of its territory. The system parameters are L = 1, 200 and ν = 48. Each subplot represents subdomain with width Ls = 50, a small fraction of the L2 domain, taken at a random location and at some random time after thermalization.
Within the range of where the bond-orientational parameter appears to indicate a hexatic phase, the system may actually be in a solid phase. This can be determined by analyzing the local translational order parameter ψt(ℓ)=exp(iG→·Δr→ℓ), where G→ is one of the two reciprocal vectors of the simple hexagonal lattice and Δr→ℓ is the displacement of the ℓ-th territorial centroid from its ideal lattice site if it were a perfect crystal (ψt(ℓ) = 1). The result of such an analysis is shown in Figure 5 where we plot the argument of the local translational order parameter arg(ψt(r→ℓ)) for each territory. The subplots point to the appearance of a solid phase for intermediate values of ρ and sufficiently large values of TA, which is when the system is also hexatic when one compares (Figure 4).
Change in the local translational order parameter ψt(ℓ) for each territory centroid located at r→ℓ as a function of ρ and TA for the system parameters utilized in Figure 4. The color of the centroid is determined by arg(ψt(ℓ)), while the size of each centroid is proportional to the size of its territory. Each subplot represents subdomain with width Ls = 50, a small fraction of the L2 domain, taken at a random location and at some random time after thermalization.
Given the limited resolution we have employed for the parameters of the system, ρ in particular, it is hard to identify the precise region in parameter space where phase transitions occur. However, by comparing Figures 4, 5 a small hexatic region may exist in a narrow region of parameter space, potentially for 1, 100 ≲ TA ≲ 1, 300 when ρ = 0.15, for 900 ≲ TA ≲ 1, 200 when ρ = 0.35, and for 900 ≲ TA ≲ 1, 200 when ρ = 0.55. This hexatic region should thus be limited on the left by the liquid-hexatic transition and on the right by an hexatic-solid transition.
To try to confirm these qualitative findings we construct the global bond-orientational order parameter Ψ6=|(1/N)∑ℓψ6(ℓ)| and the global translational order parameter Ψt=|(1/N)∑ℓψt(ℓ)| and calculate their susceptibility, χ6=N〈(ψ6−〈ψ6〉)2〉 and χt=N〈(ψt−〈ψt〉)2〉 [32, 34], respectively, as a function of TA. In Figure 6 we show the outcome of that analysis, which does not, however, help us to identify precisely the transition points. We thus turn to the analysis of topological defects to help us pinpoint the hexatic region in phase space.
The susceptibilities χ6 and χt of the global order parameters Ψ6, Ψt, respectively. A system with parameters L = 1, 200 and ν = 48 was subdivided into square subdomains of widths Lw. For each of these subdomains the global order parameters and the susceptibilities are computed and then a variance of those measurements is taken.
4.3. Topological Defect Analysis
Two types of topological defects accompany the KTHNY melting scenario, namely dislocations and disclinations. Dislocations are translational defects, which destroy long range translational order when are isolated or free. On the other hand, disclinations are orientational defects, which destroy the long range orientational order when they are free. In a perfect 2D hexagonal crystal all atoms have six neighbors. Disclination cores have only 7 or more neighbors whereas antidisclination cores have 5 or less neighbors [35, 36]. For a 2D crystal, an isolated tightly bound 5–7 fold disclination pair is a dislocation [32]. In an imperfect solid, dislocation defects are always found in pairs because they do not destroy long range order as shown in Figure C1 for our cTRW model.
According to the KTHNY theory of defect induced melting, by increasing temperature the system first undergoes a continuous transition from a solid where only pairs of dislocations are possible, to a hexatic phase where the pairs of dislocations “unbind” spawning free dislocations (see in Figure C1 an example of how a free disclination destroys the orientational order in the cTRW model). These free dislocations result in the loss of long range translational order found in crystals. With further increase in temperature the system then undergoes another continuous transition from hexatic to liquid where the dislocations themselves unbind resulting in free disclinations and the loss of long range orientational order [37, 38].
We analyse these topological defects in our cTRW model by considering their fraction λ in the system shown in Figure 7. For each type of defect λ is their number relative to the total number of centroids. Looking sequentially at panels (a) to (d) in Figure 7 one can see that as free disclinations drop significantly there is a sharp increase of dislocations signaling the onset of the liquid-hexatic transition at around, respectively, TA = 870, 890, 712, and 712. The hexatic phase appears once there are no free disclinations in the system. The faster decay of the free dislocations vs. the dislocation pairs that is observed by increasing TA is indicative of the fact that free dislocations are binding together to form pairs of dislocations. Once free dislocations disappear the system has entered the solid phase. By looking at these values of TA when free disclinations disappear (hexatic) and free dislocations disappear (solid) we are able to identify the hexatic region in the TA − ρ parameter space as shown in Table 1.
The fraction, λ of each of the defect types considered as a function of TA for thermalized systems with parameters L = 1, 200 and ν = 48. The absence of free disclinations marks the transition from liquid to hexatic phase, while the absence of free dislocations is indicative of the system having transitioned into the solid phase from the hexatic phase. The panels (A–F) correspond, respectively, to ρ = 0.00, 0.15, 0.35, 0.55, 0.75, and 0.95. Defects with different number of neighbors are at least two orders of magnitude smaller than the ones drawn and we have omitted plotting them.
Approximate phase regions in TA for each value of ρ extracted from the appearance/disappearance of topological defects (see Figure 7).
ρ
Hexatic region (TA)
Solid region (TA)
0.00
1273–1409
>1409
0.15
1157–1335
>1335
0.35
1068–1157
>1157
0.55
979–1157
>1157
Within the hexatic region there are free dislocation defects but no free disclination defects, whereas in the solid region there are pairs of dislocation defects but no free dislocation defects.
Once in the solid phase with further increase of TA the system approaches progressively the perfect crystalline arrangement by shedding paired dislocations. For the case ρ = 0.75 there is the possibility that ordered phases exist but one ought to look at values of TA beyond those considered here. On the other hand, when ρ = 0.95 the defects show no sign of decaying. This results from the breaking down of the centroid coarse-graining analysis, which is further discussed in the following section.
4.4. Territory Fragmentation
While the previous sections have shed light on the links between the movement statistics of the walkers and the phases of the territorial system for ρ ≤ 0.55, we have a less clear picture of what occurs for large ρ. While stronger correlation implies that an animal would cross its own territory more quickly, it does not necessarily mean that the terrain is covered more efficiently. On the contrary the coverage efficiency is expected to decrease beyond a certain value of ρ as demonstrated in the ideal scenario of Figure B1. With neighboring individuals not covering efficiently their own territories, an animal spending more time at the boundaries has more chance to carve away part of the terrain recently occupied by its own neighbors. In so doing the territorial shapes become less and less convex as ρ is increased.
With further increase in ρ at any given time an animal may have its own scented territory separated into multiple islands, which are sets of contiguous lattice sites that contain the scent of one walker only. Those islands where a walker is not present are cut off quickly from its owner and get absorbed by neighboring territories. While islands get created and decay away continuously, we expect the rate at which they form to become larger than the rate at which they dissolve beyond a certain value of ρ. When that happens the centroid of an animal territory could easily be located in between multiple islands, in areas that are also covered by the scent of other individuals, or even outside its own scented region when the territory shape is highly concave. In these cases the centroid analysis, utilized in previous sections to identify the phases of the system, is not appropriate any more and there are no biological meaningful ways to define a territorial centroid.
We can, however, quantify the degree of fragmentation as a function of the system parameters by considering the ratio κ = N/I where N is the number animals in the system and I is the total number of islands. This fragmentation index is plotted in Figure 8. For low TA values, traces of past presence of an individual get lost quickly, and thus it is not so frequent for an animal to encounter foreign scent. While such encounters increase as TA gets larger, potentially reducing I, they are rare and are thus relatively independent of how straight the animal movement trajectories are. Beyond a certain value of TA foreign scent encounter events become very frequent and the rate of increase is affected by the shape of the territories. For ρ ≲ 0.75 the encounters occur mainly at boundary sites and thus further increase in TA does not change the number of islands in the system, which remains at κ = 1 when the number of islands corresponds to the number of territories. This is not the case for ρ > 0.75 where the number of islands becomes larger than N.
The fragmentation index κ as a function of the active scent time TA for thermalized systems with parameters L = 1, 200 and ν = 48. With ρ ≲ 0.75, at low TA values the territories are concave and as a result are more prone to fragmentation resulting in low κ. As the TA is increased, the territories become more convex becoming more resistant to fragmentation resulting in κ = 1 for TA ≳ 800. On the other hand, for ρ ≳ 0.75 fragmentation persists in the system regardless of TA.
The appearance of fragmented territories for ρ ≳ 0.75 as shown in Figure 8 helps us interpret Figures 4, 5 when ρ ≥ 0.5, which indicated the presence of a fluid for any TA. In light of Figure 8 it is clear that we cannot conclude that for larger TA no ordered phase appears because the territory centroids are not an appropriate coarse grained representation of the animal scented territories. Although we believe that an analysis that takes into account the actual shape of the territories might help us develop other ways to identify order in fragmented terrain, we have decided not to perform such investigation here given the heavy computational cost that such a study would entail.
5. Conclusions
Territoriality is a common behavior in natural populations as it allows individual animals or family units to defend resources being food, shelter, or mates. We have studied here the case where the mechanisms through which animals defend their resources is via leaving trace of their passage on the terrain. For our theoretical analysis we have used a collective movement model recently employed to study red fox movement data [39] called the territorial random walk model. We have extended its prediction by determining how the individual movement statistics affect territory formation at the population level. And we have focused in particular on relaxing the diffusive assumption allowing the walkers to have a variable degree of correlation in their movement steps.
By studying the main parameters of the systems, namely the degree of correlation of the walk ρ and the active scent time TA, we have brought evidence pointing to the existence of the KTHNY melting scenario with liquid, hexatic and solid phases. We have identified the coverage efficiency of an individual territory, which is controlled by ρ and TA, as the micro level behavior that explains why different phases of territories are present at the macro level.
As the movement of the individuals becomes more correlated i.e., by increasing ρ from 0, they increase territory coverage. Individuals are able to re-scent more quickly and more readily defend their territory from neighbors. Territories are less mobile and the whole system becomes more ordered.
Beyond an intermediate degree of correlation, that is for high values of ρ, the walkers are not only inefficient at covering their own territory, but also spend more time at the borders. Their increased presence at the borders gives them more opportunities to take over neighboring territories. This results in large variation in the shape and size of territories and ultimately their fragmentation. In this scenario, the coarse grained approximation is no longer suitable to determine the spatial ordering.
Moreover, we have demonstrated how it may be useful to look at the presence of different types of topological defects to pinpoint phase transitions in the correlated territorial random walk model. As thermodynamic limits represent a big abstraction from the realm of biological systems we believe that this latter part of our analysis will be instrumental to determine empirically the spatial order of an actual population of territorial animals. With the recent development in animal tracking [40] it is in fact realistic to have in the near future simultaneous highly resolved movement data of neighboring territorial individuals from which defect densities and movement correlation statistics can be extracted.
Applications of this study extend beyond an ecological context. For example, in multi-robot online area coverage [31] tasks such as surveillance, search and harvesting [41–44] involve coordinating a large number of robots. In such cases the natural approach is to use a decentralized controller [45]. The observations and insights of this study may help refine stigmergic control systems that have been successfully demonstrated in previous studies [46, 47]. As an avenue for future research, it might be interesting to study the difference of the cTRW model on scale-free and small-world networks and compare it the bioinspired machine learning models used for semi-supervised learning [48, 49].
Data Availability
The datasets generated for this study are available on request to the corresponding author. Data used can be found at https://www.dropbox.com/sh/06ts5rwb922yqwk/AADeY66Q_7fqn5axMKcCK-Tra?dl=0.
Author Contributions
LG designed the study and wrote the paper. SS and AH created the computational code. SS analyzed the data.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2019.00129/full#supplementary-material
ReferencesGoychukIKharchenkoVOMetzlerR. How molecular motors work in the crowded environment of living cells: coexistence and efficiency of normal and anomalous transport. PLoS ONE. (2014) 9:e91700. 10.1371/journal.pone.009170024626511KusumiANakadaCRitchieKMuraseKSuzukiKMurakoshiH. Paradigm shift of the plasma membrane concept from the two-dimensional continuum fluid to the partitioned fluid: high-speed single-molecule tracking of membrane molecules. Annu Rev Biophys Biomol Struct. (2005) 34:351–78. 10.1146/annurev.biophys.34.040204.14463715869394KenkreVMGiuggioliLKalayZ. Molecular motion in cell membranes: analytic study of fence-hindered random walks. Phys Rev E. (2008) 77:051907. 10.1103/PhysRevE.77.05190718643102DieterichPKlagesRPreussRSchwabA. Anomalous dynamics of cell migration. Proc Natl Acad Sci USA. (2008) 105:459–63. 10.1073/pnas.070760310518182493ViswanathanGMda LuzMGERaposoEPStanleyHE. The Physics of Foraging: An Introduction to Random Searches and Biological Encounters. Carmbridge: Cambridge Univ. Press (2011).BoyerDCrofootMCWalshPD. Non-random walks in monkeys and humans. J Roy Soc Interface. (2011) 9:842–7. 10.1098/rsif.2011.058222031731Falcón-CortésABoyerDGiuggioliLMajumdarSN. Localization transition induced by learning in random searches. Phys Rev Lett. (2017) 119:140603. 10.1103/PhysRevLett.119.14060329053283GoldsteinS. On diffusion by discontinuous movements, and on the telegraph equation. Q J Mech Appl Math. (1951) 4:129–56. 10.1093/qjmam/4.2.129CodlingEAPlanckMJBenhamouS. Random walk models in biology. J Roy Soc Interface. (2008) 95:813–34. 10.1098/rsif.2008.0014BenhamouS. Of scales and stationarity in animal movements. Ecol Lett. (2014) 17:261–72. 10.1111/ele.1222524350897BrownJLOriansGH. Spacing patterns in mobile animals. Ann Rev Ecol Syst. (1970) 1:239–62. 10.1146/annurev.es.01.110170.001323BurtWH. Territoriality and home range concepts as applied to mammals. J Mammal. (1943) 24:346–52.MaherCALottDF. Definitions of territoriality used in the study of variation in vertebrate spacing systems. Anim Behav. (1995) 49:1581–97.AdamsES. Approaches to the study of territory size and shape. Ann Rev Ecol Syst. (2001) 32:277–303. 10.1146/annurev.ecolsys.32.081501.114034DaviesNBHoustonAI. Territory economics. In: KrebsJRDaviesNB editors. Behavioural Ecology: An Evolutionary Approach, 2nd ed. Oxford: Blackwell Sci. (1984). p. 148–69.GoslingLMRobertsSC. Scent-marking by male mammals: cheat-proof signals to competitors and mates. Adv Stud Behav. (2001) 30:169–217. 10.1016/S0065-3454(01)80007-3GrasséPP. La reconstruction du nid et les coordinations interindividuelles chez Bellicositermes natalensis et Cubitermes sp. La théorie de la stigmergie: essai dÕinterprétation du comportement des termites constructeurs. Insec Soc. (1959) 6:41–83.TheraulazGBonabeauE. A brief history of stigmergy. Artificial Life. (1999) 5:97–116. 10633572GiuggioliLPottsJRRubensteinDILevinSA. Stigmergy, collective actions, and animal social spacing. Proc Natl Acad Sci USA. (2013) 110:16904–9. 10.1073/pnas.130707111024082100LewisMAMurrayJD. Modelling territoriality and wolf-deer interactions. Ecology. (1993) 366:738–40.MoorcroftPRLewisMA. Mechanistic Home Range Analysis. Princeton, NJ: Princeton University Press (2006).GiuggioliLPottsJRHarrisS. Animal interactions and the emergence of territoriality. PLoS Comput Biol. (2011) 7:e1002008. 10.1371/journal.pcbi.100200821423708GiuggioliLKenkreVM. Consequences of animal interactions on their dynamics: emergence of home ranges and territoriality. Mov Ecol. (2014) 2:20. 10.1186/s40462-014-0020-725709829AlbertsAC. Constraints on the design of chemical communication systems in terrestrial vertebrates. Am Nat. (1992) 139:S62–89.Heiblum RoblesAGiuggioliL. Phase transitions in stigmergic territorial systems. Phys Rev E. (2018) 98:042115. 10.1103/PhysRevE.98.042115HalperinBINelsonDR. Theory of two-dimensional melting. Phys Rev Lett. (1978) 41:121–4.KosterllitzDJThoulessJM. Ordering, metastability and phase transitions in two-dimensional systems. J Phys C Solid State Phys. (1973) 1181:1181–203.YoungAP. Melting and the vector Coulomb gas in two dimensions. Phys Rev B. (1979) 19:1855–66.GiuggioliLAbramsonGKenkreVMParmenterRRYatesTL. Theory of home range estimation from displacement measurements of animal populations. J Theor Biol. (2006) 240:126–35. 10.1016/j.jtbi.2005.09.00216221478ChupeauMBénichouOVoituriezR. Cover times of random searches. Nat Phys. (2015) 11:844. 10.1038/nphys3413GiuggioliLAryeIHeiblum RoblesAKaminkaGA. From ants to birds: a novel bio-inspired approach to online area coverage. In: GroßRKollingABermanSFrazzoliEMartinoliAMatsunoF., editors. Distributed Autonomous Robotic Systems: The 13th International Symposium. Cham: Springer International Publishing (2018). p. 31–43.BinderKKobW. Glassy Materials and Disordered Solids: An Introduction to Their Statistical Mechanics. Singapore: World scientific (2011).FraserDPZuckermannMJMouritsenOG. Simulation technique for hard-disk models in two dimensions. Phys Rev A. (1990) 42:3186. 9904396JasterA. Computer simulations of the two-dimensional melting transition using hard disks. Phys Rev E. (1999) 59:2594.HarrisWF. Disclinations. Sci Ame. (1977) 237:130–45.DeWitR. Relation between dislocations and disclinations. J Appl Phys. (1971) 42:3304–8.QiWGantaparaAPDijkstraM. Two-stage melting induced by dislocations and grain boundaries in monolayers of hard spheres. Soft Matter. (2014) 10:5449. 10.1039/C4SM00125G24752821QuinnRAGoreeJ. Experimental test of two-dimensional melting through disclination unbinding. Phys Rev E. (2001) 64:051404. 10.1103/PhysRevE.64.05140411735920PottsJRHarrisSGiuggioliL. Quantifying behavioral changes in territorial animals caused by sudden population declines. Am Nat. (2013) 182:E73–82. 10.1086/67126023933730NathanRMGiuggioliL. A milestone for movement ecology research. Move Ecol. (2013) 1:1. 10.1186/2051-3933-1-125709815Schneider-FontánMMatarićMJ. Territorial multi-robot task division. IEEE Trans Robot Automat. (1998) 14:815–22.LarionovaSAlmeidaNMarquesLDe AlmeidaAT. Olfactory coordinated area coverage. Autonom Robots. (2006) 20:251–60. 10.1007/s10514-006-7099-7PaullLSaeediSSetoMLiH. Sensor-driven online coverage planning for autonomous underwater vehicles. IEEE/ASME Trans Mechat. (2013) 18:1827–38. 10.1109/TMECH.2012.2213607FazliPDavoodiAMacKworthAK. Multi-robot repeated area coverage. Autonom Robot. (2013) 34:251–76. 10.1007/s10514-012-9319-7MarjoviANunesJGMarquesLDe AlmeidaA. Multi-robot exploration and fire searching. In: 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS 2009. St. Louis, MO: IEEE. (2009). p. 1929–34. 10.1109/IROS.2009.5354598ScheidtDStipesJ. Cooperating unmanned vehicles. In: Proceedings 2005 IEEE Networking, Sensing and Control, 2005. Tucson, AZ: IEEE (2005). p. 326–31.StipesJHawthorneRScheidtDPacificoD. Cooperative localization and mapping. In: 2006 IEEE International Conference on Networking, Sensing and Control. Fort Lauderdale, FL: IEEE (2006). p. 596–601.de ArrudaHFSilvaFNCostaLdFAmancioDR. Knowledge acquisition: a complex networks approach. Inform Sci. (2017) 421:154–66. 10.1016/j.ins.2017.08.091LimaTSde ArrudaHFSilvaFNCominCHAmancioDRCostaLdF. The dynamics of knowledge acquisition via self-learning in complex networks. Chaos Interdiscipl J Nonlinear Sci. (2018) 28:083106. 10.1063/1.502700730180654
Funding. LG and SS acknowledge funding from, respectively, the Engineering and Physical Research Council Grant nos. EP/I013717/1 and S108151-111.