A mathematical model to simulate movement
For ‘attracting features’, such as nesting or roosting sites, we employ potential terms that are logarithmic in distance. Logarithmic potentials have been employed in diffusion models7 such as those involving long-range interactions8. The forces due to these are inversely proportional to distance from the features. Given a choice between locations, an animal would invariably drift towards ones that are closer. Additionally, they also command some influence for longer distances. We did consider alternatives such as a potential that corresponds to an inverse squared force but it diminishes much faster as the distance to the source increases. The ‘repulsive features’ such as human dominated areas are incorporated using Gaussian type potentials that would have an influence only when the animal is close to them. Such forces fall off exponentially fast as one goes away from the source location.
The corresponding Langevin equations can be written as:
$$begin{aligned} frac{dx}{dt}= & {} -gamma sum _{i}frac{ 2alpha times (x – x_{i})}{(x – x_{i})^2 + (y – y_{i})^2} nonumber &+,gamma sum _{j}Big ((x-x_{j}) e^{-((x-x_{j})^2 + (y-y_{j})^2)}Big ) + root 2 of {2D}xi _x(t) end{aligned}$$
(1)
$$begin{aligned} frac{dy}{dt}= & {} -gamma sum _{i}frac{ 2alpha times (y – y_{i})}{(x – x_{i})^2 + (y – y_{i})^2} nonumber &+,gamma sum _{j}Big ((y-y_{j}) e^{-((x-x_{j})^2 + (y-y_{j})^2)}Big ) + root 2 of {2D}xi _y(t) end{aligned}$$
(2)
where x and y denote the coordinates of an animal’s location. ((x_{i}), (y_{i})) and ((x_{j}), (y_{j})) denote locations of i attracting and j repelling features respectively. We only choose nests as points of attraction for breeding hornbills since their diurnal movements are strongly centred around the nests. The white noise terms (xi _x) and (xi _y) are Gaussian in nature and delta correlated—which means that no correlations exist between the noise values at different instances of time. (gamma ) and D denote the drift and diffusion coefficients respectively. The drift coefficient (gamma ) represents the directedness of motion, which could be interpreted as strength of bias towards/against certain features in the landscape. In contrast, D quantifies the strength of random undirected motion. The force term with coefficient (-gamma ) results from negative gradient of the logarithmic potential, whose choice we explained earlier:
$$begin{aligned} U = gamma sum _{i} log left{ (x – x_{i})^2 + (y – y_{i})^2 right} ^{alpha } ,,,. end{aligned}$$
(3)
The value of (alpha ) is determined from calculation of first passage times of the birds (discussed in the following section) and comparison of the values so obtained with observational (telemetry) data. We find that (alpha ) = 8 gives biologically sensible first passage times for hornbills (see “Calculating First Passage Times” in Methods section, Table 3 and Supplementary Tables 1, 2). If one observes an animal’s movement for a very long time, the probability of finding the animal would decrease more drastically away from a central feature for lower values of (alpha ). Such variations are captured by the steady-state probability distributions of space-use that we describe in the following section.
Fokker–Planck methods
Although the Langevin equations can generate trajectories of movement, the corresponding simulations need to be run for very long times to infer reliable information about spatial use. The time steps are further much smaller than the frequency of data recorded by the GPS. The step-lengths thus generated from simulated trajectories do not lend themselves to comparison against those from the recorded data. A convenient alternative is to solve a Fokker–Planck equation which has a direct correspondence with the Langevin equations. For our model, this takes the form:
$$begin{aligned} frac{ partial P(x,y,t)}{partial t}&= frac{partial }{partial x} left{ F_x + D frac{partial }{partial x} right} P(x,y,t) nonumber &quad +, frac{partial }{partial y} left{ F_y + D frac{partial }{partial y}.right} P(x,y,t) end{aligned}$$
(4)
where
$$begin{aligned} F_x&= -gamma sum _{i} frac{ 2 alpha times (x – x_{i})}{(x – x_{i})^2 + (y – y_{i})^2} nonumber &quad+, gamma sum _{j} (x – x_{j}) times e^{-( (x – x_{j})^2 + (y – y_{j})^2)} nonumber F_y&= -gamma sum _{i} frac{ 2 alpha times (y – y_{i})}{(x – x_{i})^2 + (y – y_{i})^2} nonumber &quad+ gamma sum _{j} (y – y_{j}) times e^{-( (x – x_{j})^2 + (y – y_{j})^2)} end{aligned}$$
(5)
The Fokker–Planck equation describes the evolution of the probabilities of occurrence over a given region. The probability distribution eventually reaches a ‘steady state’ which captures the long-term occurrence probabilities for a given bird, and it does not change beyond this point in time. This steady-state probability distribution can be computed by setting the time derivative term to zero in Eq. (4). The numerical solution of the Fokker–Planck equation involves discretizing the spatial derivatives involved. The steady state probability distribution is consequently obtained on a spatial domain of discretized grids.
Interestingly, Giuggioli et al.9 considered logarithmic potentials in their work on home range estimation, where an exponent of 8 was found to have a very similar steady state distribution to that from a harmonic potential. Harmonic potential has been utilized in analyzing home ranges of Peromyscus maniculatus10.
Using the steady-state solution of the Fokker–Planck equation, we compute the mean square displacement averaged over different possible starting locations using the steady state distribution. A discrete version of the mean-square displacement (MSD) can be defined as:
$$begin{aligned} MSD = sum _i^N langle (x – x_i)^2 + (y – y_i)^2 rangle P_{0}(x_i,y_i) end{aligned}$$
(6)
where (P_0(x_i,y_i)) is the distribution of starting locations (x_i) and (y_i) from where displacements are calculated. The inner angular brackets represent a similar weighted average of the mid-points of all grids over the steady-state probability distribution (P_{text {st}}(x,y)). Many of the grids that we define to perform simulations lie outside the known home range of the birds. The probability of choosing a starting location is defined using a Gaussian distribution centred around the nest or the most visited roost site.
The square root of the MSD defines a characteristic length scale. This could be interpreted as home range length when the steady state distribution is computed over an infinite extent9. A logarithmic potential does not lend itself to such computations since it decays much more slowly such that the characteristic length continues to grow with the size of the area considered. We evaluate the characteristic length scale (L) on a domain that is not much larger in size compared to the observed home range.
We also calculate L from empirical data by using the probability of occurrence over space inferred from two-dimensional histograms of location data. The MSD in this case is evaluated in the same vein as above but now the displacements from initial locations are weighted over the probabilities of occurrence derived from the histograms. Since these probabilities are only available for each grid, we choose only the mid-points of grids as possible locations to find the result. The starting locations are chosen from a uniform distribution over the mid-points of the grids. This is definitely a crude way of evaluating L but it does give us some way of comparing our numerical solutions against data. Finding a joint-probability distribution over the two dimensions would have been ideal but it is complicated by the fact that the distribution over space is multi-modal owing to multiple roosts for some hornbills. When inferring MSD from the location coordinates directly, it increases before saturating as the sampling frequency is decreased. For very high sampling frequency (or very small time intervals), diffusion effects dominate which leads to an almost linear increase in MSD. The effects of drift are more prominent compared to diffusion for lower sampling frequencies which marks the saturation of the MSD values10.
A first-passage time model for heterogeneous environments
The temporal information about an animal’s whereabouts is highly scrambled in the data. An important quantity of interest that could be extracted from movement data is the search time to reach a given target. A very useful measure of search times is the ‘first passage time’. Very generally, first passage time is the time taken for a given state variable to reach a particular value. In the case of animal movement, it can be interpreted as the time taken to reach a particular target location. McKenzie et al.11 derived an interesting first passage time model which had a direct correspondence with a Fokker–Planck equation. We use the prescription of Moorcroft et al.12,13 to estimate the drift and diffusion coefficients. This assumes a movement kernel that is a product of exponential distribution of step lengths and von Mises distribution for the turning angles. (This may be seen in the “goodness of fit tests” section in Methods where we assess fit of our data to claimed distributions.) It can be expressed as:
$$begin{aligned} K({mathbf{X}} ,{mathbf{X}}’ ,tau )=, & {} frac{1}{rho } f_tau (rho ) k_tau (phi ) end{aligned}$$
(7)
$$begin{aligned} {rm{where}},,,,,,,,,,,, ,f_tau (rho )=, & {} lambda e^{-lambda rho }end{aligned}$$
(8)
$$begin{aligned} k_tau (phi )=, & {} frac{1}{2 pi I_0(kappa _tau )} exp [kappa _tau cos (phi )] end{aligned}$$
(9)
Here, ({mathbf{X}} ), ({mathbf{X}}’ ) denote the current and previous locations respectively, f is the exponential distribution of step lengths (rho ) with rate parameter (lambda ) and mean (bar{rho }_{tau } = 1/lambda ), and (k_{tau }) is the von Mises distribution of turning angles (phi ). (tau ) refers to the time taken to complete a given step. The turning angles are computed with respect to the nest/roost sites. (kappa _tau ) is the concentration parameter of the von Mises distribution which signifies the departure from a uniform distribution of movement directions. The normalizing factor (I_0(kappa _tau )) is a modified Bessel function of the first kind and of zeroth order. The drift and diffusion coefficients can be reliably estimated as:
$$begin{aligned} gamma= & {} lim _{tau rightarrow 0} frac{bar{rho }_{tau } kappa _tau }{2tau } end{aligned}$$
(10)
$$begin{aligned} D= & {} lim _{tau rightarrow 0} frac{bar{{rho _{tau }}^2}}{4tau } end{aligned}$$
(11)
Employing the formalism in McKenzie11 to derive the equation for the first passage time T, we obtain the following equation:
$$begin{aligned}&gamma sum _{i} left{ frac{ 2alpha times ({mathbf{X}} – {mathbf{X}} _{i})}{(x – x_{i})^2 + (y – y_{i})^2} right} cdot nabla T nonumber &quad -, gamma sum _{j} left{ ({mathbf{X}} – {mathbf{X}} _{j}) e^{-( (x – x_{j})^2 + (y – y_{j})^2)} right} cdot nabla T nonumber &quad +, D nabla ^2 T + 1 = 0 end{aligned}$$
(12)
The terms in dot product with (nabla T) are simply the drift coefficients with spatial dependence.
McKenzie et al.11 had a simpler version of the first passage time equation that only accounted for bias towards the home range centre. The authors mention that the task of solving the first passage time equation is computationally harder with terms that account for more complex types of heterogeneities. We transform the partial differential equation in (12) into polar coordinates which simplifies the process of solving it (see First Passage Time calculation in Methods). The first passage times obtained from this solution also help us fix the value of (alpha ) in the equation above and subsequently in the logarithmic potential in (3), and in Eqs. (1) and (2). On performing this analysis for different hornbills, we see that (alpha ) = 8 works very well for them irrespective of the species and distribution of heterogeneities around them (see First Passage Time calculation in Methods). First passage times are calculated from the roosting/nesting site that lies closest to the home range centre. In case of GHNBr2, we calculate the first passage times from the approximate home range centre where no roosts exist. This ensures that most points considered for computations lie within the actual extent of the bird’s recorded locations. We used the Minimum Convex Polygon method to estimate the approximate home range centre14. This helped in identifying a location for each bird—which was a roost/nest in most cases—from where first passage times were subsequently computed. The method used for home range estimation is not relevant in the context of our proposed model and results presented, and therefore we do not consider other alternatives.
Source: Ecology - nature.com
