Nonlocal Mechanistic Models in Ecology: Numerical Methods and Parameter Inferences

: Animals utilize their surroundings to make decisions on how to navigate and establish their territories. Some species gather information about competing groups by observing them from a distance, detecting scent markings, or relying on memories of encounters with rival populations. Gathering such information involves a nonlocal process, prompting the development of mechanistic models that incorporate nonlocal terms to explore species movement. These models, however, pose analytical and computational challenges. In this study, we focus on a multi-species advection– diffusion model that incorporates nonlocal advection. To efﬁciently compute solutions for this system involving a large number of interacting species, we introduce a numerical scheme using spectral methods. Additionally, we examine the inﬂuence of various parameters and interaction potentials on population densities. Our investigation aims to provide a method to identify the primary factors driving species movements, and we validate our approach using synthetic data.


Introduction
The ability to learn from the local environment plays a crucial role in the survival of species.Various processes, such as interactions within social groups, competition with members of the same species, foraging behavior, and predator avoidance, profoundly influence how a species navigates and behaves in its surroundings [1][2][3].Understanding and predicting the territories of different species depend on comprehending these processes.A compelling example of this phenomenon can be observed in wolves, where their movements are guided by vital information, including the scent of their pack, the scents of fellow pack members, pack density, and the location of wolf pups [4].Spatial ecology has a keen focus on unraveling the mechanisms that underlie species movement and how these movements contribute to territory decisions.In light of the effects of climate change on the environment [5,6] and its impact on species' territories [7], these issues have gained increasing importance in scientific research.
Mathematical models leverage the mechanisms underlying territory decisions to predict how territories adapt in response to environmental changes, such as alterations in conspecific territories or resource density [8].Specifically, mechanistic models offer a framework to investigate territory formation, building upon theories developed by field scientists for species movement.These models provide some advantages over statistical models, as they facilitate hypothesis testing for the driving forces behind animal movement, exhibit strong predictive capabilities, and can be validated using field measurements [9].Mechanistic models have proven useful in describing various behaviors, including foraging [10,11], aggregation [12], and home ranges [13,14].They encompass movement processes, conspecific influences, environmental pressures, and space utilization [15][16][17].Earlier models, such as those proposed by [18,19], incorporated diffusion and attraction to a central location, like a den site, to describe territories across different species.Over time, these models have been enriched with additional processes to better replicate territory patterns based on empirical data, including aspects like scent marking [9] and habitat selection [20].Local advection-diffusion (AD) models have proven successful in predicting territories for species like wolves [4], coyotes [9], and meerkats [17].For instance, in the case of meerkats, location data and environmental information are integrated into a local mechanistic model, enabling parameter estimation to describe territory formation among social groups of meerkats.One crucial feature of these models, which in some cases ensures their well-posedness, is an attraction to fixed den sites.In a way, den-site attraction represents the use of very specific nonlocal information, making these models well-suited for species known to maintain fixed den sites.
Numerous mechanistic models proposed to understand animal territory formation have traditionally assumed that animals base their decisions solely on local information.However, the process of gathering information about the environment inherently involves nonlocal aspects.Animals rely on observations, sounds, and smells to gain insights into distant surroundings [21][22][23].Additionally, studies have highlighted the significance of nonlocal advection in maintaining swarm coherence [24] and ensuring well-posed solutions in certain models [25,26].Furthermore, field observations of meerkats have provided evidence that they actively avoid locations where they have encountered members of rival groups previously, further emphasizing the importance of nonlocal interactions in certain species [17,27].As a result, mechanistic models incorporating nonlocal advection have garnered recent interest [25,28,29], along with relevant references.This study focuses on a nonlocal advection-diffusion model that elucidates the territory formation of social species with subgroups.However, nonlocal models pose both mathematical and computational challenges, as classical theories may not be applicable, and numerical methods can be computationally demanding [25,30,31].Overcoming these challenges is crucial to investigating the behavior of two-dimensional solutions involving multiple interacting species and effectively connecting the model with empirical data.
In this study, we present an efficient spectral numerical scheme designed to solve a nonlocal advection-diffusion equation, specifically aimed at handling a large number of interacting species.Spectral methods have been shown to be efficient in a variety of contexts, including fluid mechanics [32], higher order reaction-diffusion equations [33], and integro-differential equations [34].We dive into investigating the impact of various factors, such as the interaction potential, environmental conditions, diffusion strength, and the number of groups, on the behavior of equilibrium solutions.Our findings reveal that a change in the variance of an interaction potential has a pronounced effect on a population's territory.Moreover, we observe that the delicate balance between diffusion and aggregation strengths significantly influences the computational efficiency of the numerical solutions.Furthermore, our approach incorporates environmental data, which plays a crucial role in influencing species movement through diverse mechanisms.To test the validity of our method, we employ synthetic data and implement a data fitting technique that employs maximum likelihood estimation and stochastic gradient descent.Through this process, we aim to discern the primary driving forces behind the movement of social groups, shedding light on the intricate dynamics governing species' behaviors.
Outline: In Section 2, we provide the motivation for our study and introduce the mechanistic model under investigation.Section 3 is dedicated to describing the numerical scheme designed for a general nonlocal advection-diffusion equation in two dimensions.Here, we explore the influence of model parameters on the behavior of solutions and discuss the computational implications of adding more groups to the system.Moving on to Section 4, we integrate synthetic data into the model and employ maximum likelihood estimation to analyze the primary driving forces behind species movement.Finally, our conclusions are presented in Section 6, summarizing the key findings and implications of our study.

Background and Model
In their work [17], Bateman and colleagues explored various local advection-diffusion models to study meerkat movement.They linked these models to spatial-temporal field data, encompassing the location of meerkats from 13 different groups.The core model they considered for N groups is as follows: Here, x ∈ Ω ⊆ R d , t > 0, and u i (x, t) = u 0 (x).The density of group i is represented by u i , D denotes the spatial diffusion rate, and C i is the velocity field.Each group is advected towards its home center by the unit vector vi .The model has two versions: one where groups move away from direct interactions with conspecifics, and another where they move away from scent markings by conspecifics.The authors incorporated data from three distinct time periods into various versions of Equation ( 1) and assessed which models better fit the data for each period.The specific conclusions of their study are not mentioned.However, one drawback of models like (1) is the reliance on artificial home centers to maintain group coherence, which is not observed in meerkat populations.To address this, the current study aims to develop a nonlocal model that eliminates the need for such artificial centers.An efficient numerical scheme is introduced to explore the qualitative behavior of the system and facilitate the data-fitting process.The proposed nonlocal mechanistic model is given by: where The model considers information gathering through local and nonlocal mechanisms and incorporates.The first term appearing on the right-hand-side of Equation ( 2) is an intra-group dispersal modulated by the parameter η.Note that the diffusion is not passive but rather a power-law, modeling an overcrowding effect, which results in higher group dispersal rates with larger population densities.The convolution terms capture long-range intra-group aggregation (governed by potential K 1 ) and inter-group repulsion (governed by potential K 2 ).These terms are responsible for maintaining group coherence and promoting segregation between groups.Note that we assume that the potentials K 1 and K 2 are symmetric and satisfy conditions that guarantee that aggregation.The term U describes the environment's favorability.As an example, for meerkat territory development, it could represent the density of sour grass.The study in [17] takes into account various environmental traits, such as sand type, interfaces between sand types, and elevation, as detailed in Section 5.1.
Table 1 offers a concise overview of the descriptions associated with the terms involved in Equation ( 2).This equation represents a competition involving various factors: the diffusion term facilitates population spread, intra-group long-range aggregation concentrates the population (assuming K 1 serves as an aggregation potential), inter-group interaction disperses territories, and the environment assigns significance to additional environmental elements.The qualitative behavior of solutions to this system will consequently hinge on which factors dominate the dynamics.We adopt the notation used in [35] and set K 1 = γ 1 K and K 2 = γ 2 K, where γ 1 and γ 2 are positive constants, and K is a kernel.Based on this notation, we formulate an associated energy functional for (2) as follows: When γ 1 = γ 2 , the energy functional is non-increasing over time [35].Equilibrium solutions of the system can be obtained by minimizing this energy, and a way to do this numerically was explored in [36].One can employ a scaling argument based on the energy functional (3) to demonstrate that when we have aggregative potentials K 1 and K 2 , the behavior of the minimizers will be influenced by the relative strengths of these terms.
Inter-group non-local segregation A variant of System (2) was applied in the study of animal ecosystems with multiple groups [37].In that context, it was assumed that each population can detect opposing populations within a local neighborhood through direct observations, interactions, communication via scent marking, or memory of past interactions with opposing populations.Another version of System (2) was introduced in [38] to model social segregation.While the case when N = 1 has received significant attention [31,[39][40][41], as well as versions for N = 2 [42,43], the model has been utilized to study diverse phenomena such as animal territories [26], predator/prey dynamics [44], and human gangs [45].Its behavior is further investigated in [36], where the authors explore methods to find equilibrium solutions for multiple groups in two dimensions, a critical aspect in connecting the model with data.Additionally, in [35], the authors investigate the efficient solution of a multi-species model using spectral methods in one dimension and two interacting species.

Numerical Method and Solutions
In this section, we introduce the numerical methods to solve System (2) and investigate how the parameters, potentials, and terms in the model affect the behavior of equilibrium solutions.This method is similar to the method used to study the solutions to two interacting species in one-dimension in [35].

Discretization Scheme
We outline the process to numerically find equilibrium solutions of System (2).We take advantage of working with a bounded, periodic grid and use the Fast Fourier Transform to move into frequency space, compute derivatives and convolutions, and use the Inverse Fast Fourier Transform to move back to physical space.This process gives us a discrete representation of the right hand side of our system, which allows us to discretize in time and numerically solve the system.
Consider the Fourier Transform in two dimensions, Because we are solving a nonlocal system of PDEs, we must address computing convolutions and derivatives efficiently.Therefore, there we employ two convenient properties of the Fourier Transform in our computations: As we are working on a bounded, periodic, M x M grid, we can discretize in x y,and t such that x m = m∆x, y j = j∆y, where m, j ∈ {0, ..., M − 1}, and t n = n∆t, n ∈ N. Thus, u(x m , y j , t n ) = u n mj .Below, we drop the dependency on n for notational convenience.For a fixed time t n , we compute derivatives and convolutions in frequency space using the Discrete Fourier Transform, which allows us to pass into frequency space.The Inverse Discrete Fourier Transform, allows us to return to physical space.We utilize Equation ( 6) to move into the Fourier domain, we use properties ( 4) and ( 5) to compute derivatives and convolutions in Fourier space, and use (7) to move back to physical space.In fact, we make use of the Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT) to move in and out of frequency space efficiently, reducing number of operations from M 2 to M log M.
We first compute the right-hand side of System (2) on the bounded, periodic grid for a fixed t = n∆t.Let u i,n mj be the value of u i at time n∆t for group i ∈ {1, ..., N} at grid point (m, j) for m, j ∈ {0, 1, ..., M − 1}, and let ûi,n kv represent the DFT for k, v ∈ {0, 1, ..., M − 1}.We begin by using the FFT to move into Fourier space to determine ûi,n kv , K, and Û.We compute the nonlocal terms in System (2) by multiplying ûi,n kj and K for each group i, utilizing the property in Equation ( 5).We remain in Fourier space to compute the gradient of the nonlocal terms and U, using the property in Equation ( 4 We use the IFFT to return the physical domain and multiply these terms by u i,n mj .We move back into frequency space in order to compute ∇ , and ∇ • (u∇(U)), using property (4).For example, Finally, we use the IFFT to move back to physical space.This gives us a discrete representation of the right-hand side for group i at time t = n∆t, where f is the process of computing the right-hand side using the FFT and IFFT, as described above.We can solve for equilibrium solutions directly by setting ∂ t u i,n = 0 in (8) for i = 1...N.To investigate the dynamics of the system, we can time step until we reach equilibrium.One can use Equation ( 8) and an ODE solver to find u i,n+1 .In particular, we use the Runge-Kutta method in MatLab.
We determine error at time n∆t by computing the L 2 -norm of f (u i,n ) for i = 1...N.For an equilibrium solution, this value should be zero.We continue time-stepping until this value is below a chosen tolerance.

Computation Times
To integrate location data into the system defined by Equation (2), we must determine equilibrium solutions for a large number of groups and a wide range of parameter values.This motivates us to explore equilibrium solutions for the system across one to seven groups and varying parameter values.In the context of one group, there is only intra-species short range diffusion and long-range aggregation with additional environmental influence.It is only when we have two or more groups that the segregation term comes into play.In this analysis, we focus on observing how the number of groups and the different parameter regimes impact the computation times.To streamline our discussion moving forward, we will work with a special form of Equation (2): Note that we are assuming that In the simulations, we manipulate the parameters η, a, and b.Recall that η represents the strength of diffusivity.The two new parameters a and b are constants multiplying U and K, respectively, measuring the environment and interaction potential strength.Regarding the environmental potentials employed in these trials, we use a Gaussian function defined as U(x, y) = Ae − 1 10 (x 2 +y 2 ) , where A is a normalizing parameter ensuring that U integrates to one.The interaction potential is given by K(x, y) = Be − (x 2 +y 2 )

2
, where B is determined such that K integrates to one over the specified domain.The simulations share a common initial condition for each trial, denoted as u i,0 = e −(x−x i 0 ) 2 −(y−y i 0 ) 2 .For N = 1, the initial position is set as (x 0 , y 0 ) = (0, 0).For N > 1, the initial positions (x i 0 , y i 0 ) are symmetrically distributed about a circle centered at the origin with a radius of five.The simulations were executed until the error for each equation in System (2), as described in Section 3.1, dropped below a predefined tolerance value of tol = .0005,chosen after running trials and calculating the difference between consecutive solutions.
Table 2 displays the computation times required to compute equilibrium solutions for various combinations of η, a, and b when considering different numbers of groups, specifically N = 1, 3, 5 and 7.As expected, as the number of groups increases, the computational expense also rises.However, it is noteworthy that the computational cost is significantly lower compared to the findings in [36], where it took hours to compute equilibrium solutions in two dimensions with N = 1.Furthermore, certain combinations of η and b lead to more efficient computations for higher values of N. When the parameter b is larger, indicating stronger aggregation within groups, a higher value of η helps balance this effect.An observation from Table 2 is that lower η values take more time to compute when b = 0.25, compared to when b = 0.2.In cases where η and b are well-balanced, the expense of adding groups remains relatively constant.However, in situations where η is small relative to b, adding groups becomes more computationally expensive.A factor influencing the computational duration is the proximity of the initial data to the equilibrium solution.Additionally, we observe that when the value of b was excessively large compared to η, or if the value of a was too large, our numerical techniques exhibited instability and failed to reach convergence.The balance between η and b not only impacts the computational cost of solving the system, but also influences the behavior of solutions.For instance, Figure 1 displays equilibrium solutions for the scenario when b = a = 0.25, corresponding to the bottom-left quadrant of Table 2.In this context, the first column (representing more aggregated groups) takes the longest time to compute compared to other simulations with the same number of groups, while the second and third columns are the fastest, but similar to the remaining four.Note that the increase in the parameter η leads to more overlap between territories of competing groups., where B is determined such that K integrates to one over the specified domain, with varied values of N and η.

The Effect of Parameters on Equilibrium Solutions
In this section, we explore the influence of the overcrowding effect, interaction potentials, and environmental force on the qualitative properties of the territories of each group.The equilibrium solutions are obtained using the techniques described in Section 3.1, with various parameters being adjusted.All simulations utilize the Gaussian potential, given by K(x, y) = Be − (x 2 +y 2 ) 2σ 2 , where we vary the variance σ 2 .The parameter B is chosen such that the integral of K over the domain equals one.Additionally, we also assume that the environment is given by a Gaussian potential.Similar to previous simulations, we have parameters a and b, which are multiplied by U and K, respectively, to control the strength of the environment and the interaction potential.

Single Group Dynamics
We start with a single group to isolate the impact of changing potentials and parameters on a group before introducing inter-group interactions.As previously discussed, the parameter η influences the diffusivity strength within the group, which competes with the force of the nonlocal term, K * u, responsible for group aggregation.The interplay between these two terms determines the group's territory density and size.The top row of Figure 2 presents an illustrative example of how equilibrium solutions change when the diffusion strength, η, is varied while keeping a and b fixed.This means that both the environmental and interaction potential strengths remain constant.As anticipated, the territory size increases with η; for instance, in Figure 2c, when η = 9, the territory is more spread out compared to Figure 2a,b, where η = 5 and η = 7, respectively.The second row of Figure 2 shows a variation in the strength of the interaction potential.In Figure 2d-f, the parameter b is adjusted, while all other parameters remain fixed.It is important to note that the population densities are displayed relative to the environment, and as b increases, a clear aggregation of the population can be observed.The parameters of the potential also play a crucial role in influencing the aggregation of population densities.Figure 3 provides a visual representation of the solution for a single group using Gaussians with varying variances.It is evident that larger variances lead to a more spread-out population territory.
The territory arrangement of the groups can also be influenced by environmental factors.Figure 4 illustrates the equilibrium population density of a group existing in a fragmented environment.This group is subject to an interacting potential denoted as K, following a Gaussian distribution with σ = 1.The environmental potential, U(x), is composed of four Gaussians specified by U(x) = A[e − 1 2 (x 2 +(y−5) 2 ) + e − 1 2 ((x−5) 2 +y 2 ) + e − 1 2 (x 2 +(y+5) 2 ) + e − 1 2 ((x+5) 2 +y 2 ) ], and it is normalized to integrate to one.The upper row offers an bird's eye view of the territories, while the lower row depicts a three-dimensional view showcasing both the environmental potential and the population density.All simulations showcased in Figure 4 utilize the parameters η = 3 and b = 0.5.The two columns present distinct strengths of the environmental potential.Figure 4a,c illustrate different viewpoints of the same simulation, with a value of a = 0.25.Remarkably, in this instance, the environment is less attractive, leading to territory aggregation.As the environmental potency intensifies, as demonstrated in Figure 4b,d, the dominance of the environmental potential surpasses the aggregation strength within the group, ultimately leading to territorial segregation.

Multi-Group Dynamics
Next, we look at the scenario involving multiple groups.Figure 5 illustrates the case of three groups with U(x) = Ae −0.05(x 2 +y 2 ) and K a Gaussian with σ = 1, η = 5 and a = 0.25.Figure 5 provides two views of the solution when the interaction potential strength is at its weakest and b = 0.25.We observe that as b increases the territories become more segregated.

Incorporation of Synthetic Data
We are able to use synthetic data obtained from equilibrium solutions to System (2) to explore parameter inference via maximum likelihood estimation.In this section, we use the same environment U and potential K as was used Figure 1.We generate synthetic data by rejection sampling.Specifically, we choose the coordinate pair (x i , y i ) at random from the domain and z i at random from a uniform distribution on the interval (0, 1).If f (x i , y i ) > z i , where f is the probability density function found from the equilibrium solution to (2), we add the data value (x i , y i ) to the set.We continue until we have the desired number of data points.
Given the set of synthetic data with M data points for each of N groups, {y m i |i = 1, ...N, m = 1, 2, ..., M}, we minimize the negative log-likelihood function over the set of parameters, θ: The probability density functions, f (x|θ), are found using the methods described in Section 3, with initial condition determined by the kernel density estimation of the set of synthetic data.We use stochastic gradient descent (SGD) to minimize the negative log-likelihood function.SGD uses one or a few data points to update the parameter at each iteration.This leads to more variance in the update, thus allowing a possibility to move away from a local minimizer.We update the parameter as follows: where m is a randomly chosen from the set {1, 2, ..., M}, and α is the learning parameter, which typically decreases with time.To keep within a stable parameter range, we have the bounds 0.125 ≤ b ≤ 0.44 and 3 ≤ η ≤ 8.

Data Fitting with One Parameter
In this section, we explore parameter estimation, with an emphasis on a single parameter.To assess potential variations in the estimation of different parameters, we have opted to separately investigate the estimation of the parameters η and b.Given that territory formation hinges on the interplay between aggregation and diffusion, and our chosen territory primarily serves the purpose of aggregating the population, we opt to exclude the parameter a.

N = 1
Figure 6 illustrates an instance of fitting parameter b to synthetic data for a particular group through the aforementioned methodologies.The process of Stochastic Gradient Descent (SGD) is depicted in Figure 6a, showcasing the progression of negative log-likelihood values and the corresponding estimated parameter values throughout each iteration.The equilibrium solution employed in generating the data are represented in Figure 6b, while the equilibrium solution derived from the parameters estimated via SGD is displayed in Figure 6c.The inferred equilibrium solution closely aligns with the genuine equilibrium solution's territory, although the latter exhibits greater symmetry and a slightly heightened peak value.To assess convergence frequency and the impact of dataset size, we conducted 60 maximum likelihood estimation trials to fit b using data generated from b = 0.25.We explored three dataset sizes: M = 100, M = 200, and M = 300, presenting the outcomes in Figure 7.A significant majority of trials clustered around the true parameter value of 0.25, with fewer deviating.This trend intensified as dataset size increased, supported by decreasing variances: 0.0085, 0.0072, and 0.0060, respectively.In a similar vein, Figure 8 fits the parameter η using data generated from η = 5.The iterations of SGD are depicted in Figure 8a, the equilibrium solution employed in generating the data are represented in Figure 8b, while the equilibrium solution derived from the parameters estimated via SGD is displayed in Figure 8c.As we saw previously, the estimated equilibrium solution closely aligns with the true equilibrium solution's territory with less symmetry.In Figure 10, we estimate b from data generated with b = 0.25, η = 5, a = 1 in the case of two interacting species.In this case, the algorithm does not converge to the parameter the data were generated with, but converges to a parameter value that has a lower log-likelihood value.The parameter b is slightly overestimated, which leads to a more aggregated estimated territory in Figure 10c than the genuine territory depicted in Figure 10b.Additionally, the estimated territory has a higher maximum value and is less symmetric than the true territory.In Figure 11, we similarly estimate η from data generated with η = 5, b = 0.25, a = 1, and N = 2.The algorithm underestimates η and converges to a lower log-likelihood value than that of the true parameter value.This underestimation of η is in line with the pattern we saw in Figure 9.It is again worth comparing the territory used to generate the data with the predicted territory.Figure 11c shows the equilibrium solution found using the parameters predicted from SGD, and Figure 11b shows the equilibrium solution used to generate the data.The territory boundaries are similar, with the true territory being more symmetric, more spread out, and having a slightly lower maximum value.

Data Fitting with Two Parameters
In this section, we fit two parameters to the two-species model.As in previous sections, we fit parameters to synthetic data generated for N = 1 and N = 2, but in this case, we fit both b and η to the data.

N = 1
When we simultaneously adjust both b and η to match the model, our focus lies in the ratio of these parameters.This emphasis arises from the fact that the specific environment used induces territory aggregation as b increases, whereas territory dispersion occurs as η increases.This behavior can be seen mathematically through a scaling argument applied to the energy function (3).As a consequence, overestimating or underestimating both parameters tends to offset each other and results in a comparable territory outcome.In Figure 12, we fit the model with data generated using η = 5 and b = 0.25.Figure 12a illustrates that the algorithm converges to approximately η = 4.25 and b = 0.19, with both parameters being underestimated.However, while a lower b leads to less aggregated territory, a lower η also results in a less spread-out territory.These effects can counterbalance each other, resulting in a territory resembling the one used to generate the data.The equilibrium solutions for System (2) with the parameters used for data generation and the estimated parameter set are depicted in Figure 12b,c.These parameter sets yield similar territory boundaries, highlighting this phenomenon.

N = 2
In Figure 13, we demonstrate fitting b and η to synthetic data when N = 2. Figure 13a illustrates that b and η are both overestimated.As we have seen previously, overestimating both parameters can balance out and lead to a comparable territory to the true equilibrium.This is demonstrated when comparing the true territory, Figure 13b, to the estimated territory, Figure 13c.The estimated territory is similar to the true territory, it is more aggregated and has a higher maximum value.As is true of other estimated territories, the genuine equilibrium solution is also more symmetric.

Real World Application
In this section, we present a small case study centered on meerkat population territories within the Kalahari Desert.Meerkats (Suricata suricatta), belonging to the mongoose family, are well-known for their social behavior.They live in small social groups known as clans, typically consisting of an average of 20 members [46,47].Competition among these clans plays a pivotal role in the social structure of meerkats.This competition arises due to the scarcity of resources within the desert environment, with each clan needing to safeguard its territory.These territories encompass burrow systems, foraging areas, and essential water sources.Meerkats employ mobbing calls as a means of intra-clan communication regarding the presence of rival clans [48].Meerkats exhibit a preference for landscapes with sparse vegetation, as it enhances their ability to defend against potential predators.Furthermore, their inclination toward specific sandy environments is a focal point of our exploration here.

Non-Smooth and Non-Periodic Environments and Location Data
In practical applications, environmental data often exhibit irregularities and a lack of periodicity.In this section, we address these challenges and demonstrate our approach using relevant environments for meerkats from the study conducted by Bateman et al. (2015) [17].We refer the interested reader to that study for detailed information about the data acquisition process for the meerkat habitat [17].Notably, the type of sand and the interfaces between different sand types are believed to influence meerkat movement.The environmental data collected from meerkat habitats in South Africa, representing sand types and sand-type interfaces, are depicted in Figure 14a,b, respectively.For ease of analysis, the data has been normalized to a range between zero and one.In the case of sand-type data (Figure 14a), meerkats seem to prefer ferrous sand over clay sand, thus zero corresponds to clay sand, and one corresponds to ferrous sand.In the case of edge data (Figure 14b), zero represents the interface between different sand types, while one represents a distinct sand type.To address the irregularities in the environmental data, we apply a smoothing technique by convolving it with an approximation to the identity, resulting in a regularized representation denoted as U.This regularization process ensures uniformity throughout the data while preserving its influence on animal movement.The regularized environment functions are visualized in Figure 14c,d.Additionally, it is unsurprising that meerkat habitats lack periodicity, which poses a challenge for our numerical scheme that assumes periodic boundary conditions.To mitigate boundary effects, we extend our domain and set U(x) = 0 for locations x in the extended part of the domain.This extension is depicted by the blue border in Figure 14c,d.We incorporate meerkat location data from the Kalahari Meerkat Project for illustrative purposes.The data, as presented in [17], cover three distinct time periods: a steady period characterized by relatively stable territories, a period during which a social group divided into two new groups, and a period when territory locations underwent shifts.The dataset encompasses observations from thirteen meerkat groups, where each data point represents the time and location of a sighting of a group member.Figure 15a visualizes the location data from the steady period.In order to utilize initial conditions based on this dataset, we conduct kernel density estimation for each group.To achieve convergence in numerical simulations with these initial conditions, we mollified the kernel density estimation for each group and added a small constant of 0.001 to avoid densities close to zero.Note that this shifts all densities and has no bearing on the results.The initial conditions with the corresponding data are shown in Figure 15b.Please note that Figure 15b is presented without the extended domain.However, when utilizing the extended domain, the initial values for each group are further adjusted away from the boundaries, thereby reducing the impact of periodic boundary conditions.

Environment and Location Data
Finally, we find equilibrium solutions informed by the data we intend to incorporate into the system.This allows us to understand challenges that come with using more complicated environments and computation times associated with more detailed environments and large N. We use the location data, Figure 15a, take the kernel density estimation of the data, mollify it, and then add a small constant to avoid proximity to zero.Once again, this has no bearing on the results.The result is illustrated in Figure 15b.We use this as an initial condition.Then, for a fixed set of parameters, we find the equilibrium solutions using the two environment types, mollified with extended domains, U = EDGE, shown in Figure 14d, and U = SAND, shown in Figure 14c.
Figure 16 shows equilibrium solutions with two different environments, as well as the location data.The figures were both generated with the same parameter values, but Figure 16a was generated with the environment seen in Figures 14d and 16b was generated with the environment seen in Figure 14c.Although similar, these different environments lead to different territories.The negative log-likelihood value was lower with a value of 1.2232 × 10 4 for the EDGE environment data than the SAND environment, which had a negative log-likelihood value of 1.2284 × 10 4 .

Discussion and Conclusions
Nonlocal mechanistic models provide a more realistic approach for modeling species with consideration of nonlocal information during territory formation.However, these models come with higher computational costs compared to local models.Incorporating data into such models requires multiple iterations across parameter space, especially for complex systems with multiple species.Thus, achieving efficient solutions becomes crucial.In this study, we use a nonlocal mechanistic model to depict territorial behaviors reliant on nonlocal information.Our approach efficiently handles system complexities using spectral methods, known for their efficiency in numerical schemes.We investigate computation times across different parameter settings and species numbers, examining the influence of parameters and system terms on solutions.The results in Table 2 indicate that finding the equilibrium for seven species takes approximately one hour.We find that balanced diffusion and aggregation lead to quicker equilibrium solutions.
Utilizing our numerical method, we generate synthetic data and use maximum likelihood estimation to infer model parameters.We perform parameter estimation for cases involving one or two parameters and one or two species.Multiple trials enhance estimation accuracy, with larger datasets narrowing parameter prediction ranges.Our primary findings emphasize the efficiency of spectral methods for nonlocal systems, coupled with the feasibility of data incorporation through minimizing the negative log-likelihood function using stochastic gradient descent.However, a limitation is the need for periodic boundary conditions, which may lack physical significance.We discuss potential remedies for this limitation.Looking ahead, future research prospects include integrating meerkat location data and environmental information into the model, along with implementing model selection.

Figure 2 .
Figure 2. The equilibrium solutions to System (2) with U(x) = Ae −0.05x 2 −0.05y 2 , where A is chosen so U integrates to 1, and K is a Gaussian with σ = 1.The simulations illustrated in the top row fix b = 0.375 and vary the value of η.The bottom row fixes η = 5, and varies b.In all simulations, the environment U is displayed in black and the population density in teal.

1 Figure 3 .
Figure 3.The equilibrium solution to a single group with U = 0, K a Gaussian, η = 7, b = 0.25, with σ varied.The top row illustrates the territory of the population, and the bottom row illustrates the interaction potential.

Figure 5 .
Figure 5.The equilibrium solutions to System (2) with N = 3, U(x) = Ae −0.05(x 2 +y 2 ) , normalized to integrate to one, K a Gaussian with σ = 1, η = 5 and a = 0.25.The top row illustrates a bird's eye view of the equilibrium solution, while the bottom row provides a three-dimensional view with U shown in black and the population density in teal.

Figure 6 .
Figure 6.Iterations of SGD minimizing the negative log-likelihood function to estimate b using data generated with M = 300, N = 1, b = 0.25, a = 1, and η = 5, and the equilibrium solutions to System (2) using the estimated parameter and the parameter used to generate the data.(a) Iterations of SGD and corresponding estimated b (top) and negative log-likelihood value (bottom).(b) PDF used to generate data.(c) PDF from estimated parameters.

Figure 7 .
Figure 7.The histogram of estimated b values from 60 trials of SGD for various values of M, with data generated from b = 0.25, a = 1, η = 5.The average values from the simulations are 0.2537, 0.2342, 0.2464, respectively.The variances from the simulations are 0.0085, 0.0072, 0.0060, respectively.

Figure 8 .
Figure 8. Iterations of SGD minimizing the negative log-likelihood function to estimate η using data generated with M = 300, N = 1, b = 0.25, a = 1, and η = 5, and the equilibrium solutions to System (2) using the estimated parameter and the parameter used to generate the data.(a) Iterations of SGD and corresponding estimated η (top) and negative log-likelihood value (bottom).(b) PDF used to generate data.(c) PDF from estimated parameters.

Figure 9
Figure 9 analogously assesses convergence frequency and the effect of data set size for η.Via 60 maximum likelihood estimation trials, we fit η to data generated from η = 5 for three different data set sizes.Average values closely aligned with the true parameter although, interestingly, trial averages consistently underestimated the actual parameter value.Furthermore, as dataset size increased, variance diminished: 0.1171, 0.0772, and 0.0462, respectively.

Figure 9 .
Figure 9.The histogram displays estimated η values from 60 SGD trials for different M values, using data generated with b = 0.25, a = 1, and η = 5.The average values across simulations are 4.7295, 4.6357, and 4.6667, with corresponding variances of 0.1171, 0.0772, and 0.0462 for each case.4.1.2.N = 2In Figure10, we estimate b from data generated with b = 0.25, η = 5, a = 1 in the case of two interacting species.In this case, the algorithm does not converge to the parameter the data were generated with, but converges to a parameter value that has a lower log-likelihood value.The parameter b is slightly overestimated, which leads to a more aggregated estimated territory in Figure10cthan the genuine territory depicted in Figure10b.Additionally, the estimated territory has a higher maximum value and is less symmetric than the true territory.

Figure 10 .
Figure 10.Iterations of SGD minimizing the negative log-likelihood function to estimate b using data generated with M = 150, N = 2, b = 0.25, a = 1, and η = 5 and the equilibrium solutions to System (2) using the parameter used to generate the data and the estimated parameter.(a) Iterations of SGD and corresponding estimated b (top) and negative log-likelihood value (bottom).(b) PDF used to generate data.(c) PDF from estimated parameters.

Figure 11 .
Figure 11.Iterations of SGD minimizing the negative log-likelihood function to estimate η using data generated with M = 150, N = 2, b = 0.25, a = 1, and η = 5 and the equilibrium solutions to System (2) using the parameter used to generate the data and the estimated parameter.(a) Iterations of SGD and corresponding estimated η (top) and negative log-likelihood value (bottom).(b) PDF used to generate data.(c) PDF from estimated parameters.

Figure 12 .
Figure 12.Iterations of SGD minimizing the negative log-likelihood function to estimate η and b using data generated with M = 300, N = 1, b = 0.25, a = 1, and η = 5 and the equilibrium solutions to System (2) using the parameters used to generate the data and the estimated parameters.(a) Iterations of SGD and corresponding estimated b and η (top, middle) and negative log-likelihood value (bottom).(b) PDF used to generate data.(c) PDF from estimated parameters.

Figure 13 .
Figure 13.Iterations of SGD minimizing the negative log-likelihood function to estimate η and b using data generated with M = 300, N = 2, b = 0.25, a = 1, and η = 5 and the equilibrium solutions to System (2) using the parameters used to generate the data and the estimated parameters.(a) Iterations of SGD and corresponding estimated b and η (top, middle) and negative log-likelihood value (bottom).(b) PDF used to generate data.(c) PDF from estimated parameters.

Figure 14 .
Figure 14.The meerkat environmental data and their regularized counterpart on an extended domain.

Figure 15 .
Figure 15.(a) The location data of meerkats are presented within the same geographical area as depicted in Figure 14.Each data point corresponds to a sighting of a meerkat from a specific group, and the different colors represent distinct meerkat groups.(b) The distribution approximation of the meerkat location data via mollified kernel density estimation added to a small, positive constant.

Table 1 .
Physical description of the terms in the model given by Equation (2).

Table 2 .
The computation time (in seconds) needed to find the equilibrium solutions to System (2) with varying environment strength (a), potential strength (b), number of groups (N), and diffusivity parameter (η).Note that the number of groups vary down the second column and the value of η varies across the second row.