Biomimicry at the Landscape Scale: Agent-Based Model Simulating Beaver-Inspired Construction

1. Introduction

Landscape morphology is the emergent result of not only geological processes but also continuous, reciprocal interactions between biological agents and their physical environment. Thus, landscapes are indeed adaptive systems, shaped by environmental change, as well as by the site’s inhabitants, human and non-human alike [1,2,3]. This work extends biomimicry from tools to biomimicry at the landscape scale, emulating both spatial and temporal aspects of design through multi-agent simulations. While Agent-Based Modeling (ABM) has been successfully utilized to simulate human-environment dynamics in urban settings [4,5,6], its application to understanding non-human geomorphological change remains an underexplored frontier.

In this context, ABM provides a structured framework for understanding how specific animal behaviors, modeled as Multi-Agent Systems (MAS) [7], influence the development of adaptive landscapes [8,9]. Designing for adaptive landscapes requires tools and methods that are process- and not merely form-based.

Emergence, a condition in which large-scale patterns arise from simple local interactions, offers a powerful paradigm for environmental design [10]. Rather than relying solely on static, prescribed spatial layouts, an emergent design approach allows researchers and designers to interface directly with ongoing ecological dynamics. By understanding landscapes as dynamic, process-driven entities, we can build in conversation with natural forces—such as seasonal flooding, sedimentation, and vegetation growth—to promote biodiversity and carbon sequestration. In this context, Agent-Based Modeling (ABM) serves as a critical computational tool. Because ABM explicitly translates localized, individual behavioral rules into macro-scale geomorphological patterns, it provides the precise quantitative framework needed to simulate, test, and leverage these emergent ecological processes.

As a clear example, the logic of emergence characterizes the ways North American beavers (Castor canadensis) change landscapes. These remarkable ecosystem engineers [11] construct various structures (dams, lodges, trails, canals), terraforming grounds and water systems as a living example of a process-driven design approach [12,13,14,15]. Collectively, beavers engineer a living, adaptive infrastructure. How might we model and simulate beavers’ transformations? This paper presents a bio-inspired ABM framework to simulate the dynamic, emergent behavior of beavers in landscapes. Indeed, reading and modeling landscapes in a bio-inspired way should capture their adaptive and emergent nature. Complex Adaptive Systems (CAS) [16] are responsive to disturbances, evolving through environmental change, and capable of self-organization. These systems are inherently heterogeneous and continuously updated, capable of absorbing or transforming under dynamic conditions and shifting forces [17].

1.1. Scope of the Work and Methodology

In this work, we propose that the logic of emergence offers a powerful paradigm for understanding, preserving, and designing such adaptive landscapes. We specifically investigate the CAS evolution dynamics through the lens of the non-human agents that continuously and actively shape these environments. By studying the individual behavioral policies that drive collective behaviors, we can gain profound insights into how environmental structures emerge, and how inherent constraints, such as slope, vegetation, and streams, can be efficiently exploited for synergic and sustainable landscape maintenance [18,19,20]. By exploiting the growing availability and accuracy of Geographic Information Systems (GIS) [21], we aim to translate the adaptive strategies of animals [22] into measurable, explainable parametric landscape models.

ABM provides a structured framework for understanding how the behaviors of specific animal groups, modeled as MAS [7], influence adaptive landscapes and geomorphology [8,9,23,24,25,26]. Crucially, recent research highlights a profound convergence between geomorphological transformation and spatial networks [27,28]. This active field utilizes positional graph theory to describe emergent landscapes, shaped mainly by animal agents but also by other external forces, as dynamic topological structures [29]. While much of the existing literature focuses on static descriptions, employing metrics to identify formal properties such as small-world or scale-free characteristics in natural formations [27], growing attention is being directed toward the dynamic evolution of these networks [30].

1.2. Case Study: The Beaver Damming Complex

The activity of the North American beaver represents a primary exemplar of process-driven landscape transformation. Widely recognized as a prolific “ecosystem engineer” [12,13], the beaver constructs a living, adaptive infrastructure that emerges through iteration, feedback, and persistence. While traditional research often focuses on the single-dam framework, our approach reframes beaver building behavior within the broader context of the beaver damming complex (BDC), a landscape-scale ($\approx 1{\mathrm{km}}^2$) network comprising dams, lodges/burrows, trails, canals, scent mounds, and food caches [13,31].

This distributed architectural network fundamentally restructures regional hydrology, vegetation, and biodiversity [32]. Notably, while dams and the resulting ponds are the most visible markers of this engineering, the trail and canal systems that constitute fundamental infrastructure of the BDC remain comparatively understudied [32,33]. These pathways serve critical logistical functions, facilitating the transport of harvested vegetation to central locations for dam material and food storage. The scale and complexity of these decentralized construction processes offer a profound model for adaptive design in human-engineered landscapes.

1.3. Contribution and Long Term Objectives

This paper develops a multi-agent-based simulation framework to investigate how individual beaver behaviors, coupled with environmental feedback [20], drive the emergence of trails and canals [10,34]. The contributions are twofold: (1) we present a simple, open-source simulation suite that integrates Digital Elevation Models (DEMs) and aerial RGB imagery [35,36] with active beaver agent teams, and (2) we introduce biologically plausible behavioral policies that give rise to features qualitatively corresponding to the search trails and persistent ponds observed in the natural systems.

The motivation for this work is inherently interdisciplinary, bridging behavioral ecology, landscape architecture, and wildlife management [10,18]. By releasing the simulation suite code as an open-source testbed [37], we invite the community to extend these behavioral policies and test them against diverse terrains.

While this paper focuses on the simulation framework and emergent morphologies, it serves as a platform for a broader discussion on formal landscape metrics. We envision a shift toward topology-inspired methods, such as positional graph theory or spatial correlation indices, to formally quantify the structure of emergent networks. While a full analysis of these metrics is beyond the primary scope of this work, we provide an initial discussion on how such approaches could eventually validate simulated results quantitatively against real-world BDC configurations.

Furthermore, while the North American beaver serves as our primary case study, the broader ambition of this work is to establish a formal methodology for both targeted quantitative analyses of landscapes shaped by natural and/or artificial agents and the application of animal-inspired interventions to shape landscapes to serve desired ends. We envision a future where these ecologies shaped by animal activity are not only preserved but actively synthesized through the deployment of robotic agents. The latter aim will build on how, by leveraging on-site robotic tools for precise fabrication [38], it becomes possible to modulate grounds with a high degree of sensitivity to matter and form.

The remainder of this paper is organized as follows: Section 2 details the MAS architecture and behavioral policies. Section 3 presents the GIS data acquisition pipeline and the emergent landscape morphologies that result from application of the model. Lastly, Section 4 explores the potential for formal topological metrics and the future of process-based environmental modeling.

2. Materials and Methods

As introduced in Section 1.3, this work develops a multi-agent simulation to investigate the behavioral mechanisms underlying trail formation and canalization in beaver-like agents [39,40]. The simulation is built using the Mesa framework [41] and can be deployed on both DEMs and RGB aerial imagery of real landscapes. The test cases presented in this work utilize data obtained from the National Oceanic and Atmospheric Administration (NOAA) [42] and the United States Geological Survey (USGS) [43]. These datasets are supplemented by experimental data acquisition campaigns conducted across three US-based sites: Acadia National Park (Maine) [44], Rumney Marsh Reservation (Massachusetts) [45], and Rumney Ranch (Montana) [46]. Similar to [9], the framework is structured into two primary modules: (1) the Environment and (2) the Agent [37].

2.1. Dynamic Environment Backend via Empirical Geospatial Data

The BeaversEnvironmentBackend class transforms raw geographic data into an $N\times M$ matrix representing spatial vegetation quality. The choice between the DEM and RGB pipelines depends on the specific structural features required for the simulation:

• DEM and Topographic Map Pathway: Provides topographic robustness and vertical accuracy at a 1.0 m resolution. While the combination of DEM and USGS topographic map data effectively defines the site’s hydraulic constraints, it lacks fine-scale surface features. Consequently, foraging trails and micro-vegetation patterns are supplemented using National Agriculture Imagery Program (NAIP) datasets [43].
• RGB Pathway: Leverages high-resolution aerial imagery (1.4–1.7 cm/pixel) to resolve fine-scale features, such as foraging trails and localized vegetation patches. However, this pathway is sensitive to illumination and weather noise, requiring more intensive preprocessing to isolate water bodies and normalize vegetation quality metrics.

2.1.1. DEM-Based Processing

The DEM pipeline ingests georeferenced CSV elevation points and resamples them onto a regular grid via interpolation. To recover the hydrological network in datasets lacking water-class annotations, a percentile-based thresholding step is applied. Let Z be the elevation field and $P_q\left(Z\right)$ its q-th percentile. A water mask $\mathcal{W}(x,y)$ is generated where:

$$\mathcal{W}(x,y)=\left\{\begin{array}{cc}1,\hfill & Z(x,y)\le P_q\left(Z\right)\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(1)

The parameter q is specific to each dataset and must be calibrated to ensure channel continuity while avoiding over-detection in low-lying areas.

This process is illustrated here using the Rumney Marsh and Acadia sites. For Rumney Marsh, water pixels are already labeled as such in the CSV file; the recovery of the channel network is shown in Figure 1, where the initial Area of Interest (AOI) selection is compared against the final normalized matrix. The sensitivity of the threshold q is demonstrated using the Acadia dataset in Figure 2, showing the progression from conservative to aggressive stream masks. The final matrix maps land to $[0,1]$ and water bodies to $[-1,0)$.

2.1.2. RGB-Based Processing

The RGB pipeline utilizes a multi-temporal aggregation strategy to mitigate the higher noise levels. A baseline is constructed by aggregating multiple images capturing the same geographic area across different time points and lighting conditions. These images are merged into a reference baseline using either pixel-wise median or mean filtering (user-defined) to attenuate transient artifacts such as moving shadows or clouds.

To delineate land from water, the framework supports a range of indices, including the Excess Green index (ExG), the Normalized Difference Vegetation Index (NDVI) visible approximation, and the Color Water Index (CWI). The user can visualize these indices and their distribution histogram, in order to select the thresholded index that provides the best water-land separation for a given site. Target images are shadow-corrected and histogram-matched to the baseline to ensure index consistency. This methodology is applied here to a specific area of the Rumney Ranch site (referred to as A2), as shown in Figure 3. Specifically, we consider a set of 10 images of the A2 region spanning June to August 2018, using ExG to identify water and generate the final vegetation-quality matrix. Due to space constraints, the full image set for the A2 area is not shown here, but is available on GitHub.

A comparison between the shadow-corrected baseline (Figure 3(a)) and the original target image (Figure 3(b)) shows the necessity of the multi-temporal approach; the single target acquisition suffers from non-uniform illumination, particularly in the bottom-right corner. Furthermore, the normalized ExG matrices validate the pipeline’s ability to resolve biological features; in the original target image (Figure 3(c)), faint trail tracks are visible in the bottom-left quadrant. Crucially, they remain distinguishable also in the bottom-left quadrant of the final simulation-ready grid (Figure 3(d)).

Figure 3. Rumney Ranch A2 RGB preprocessing. (a) Aggregated baseline for noise and shadow attenuation; (b) target acquisition displaying illumination gradients; (c) baseline vegetation-quality reference; (d) final normalized matrix highlighting high-frequency trail features in the bottom-left.

Figure 3. Rumney Ranch A2 RGB preprocessing. (a) Aggregated baseline for noise and shadow attenuation; (b) target acquisition displaying illumination gradients; (c) baseline vegetation-quality reference; (d) final normalized matrix highlighting high-frequency trail features in the bottom-left.

2.2. Environmental Representation and Dynamics

The environment is represented as a 2D matrix $\mathcal{E}\in {\mathbb{R}}^{N\times M}$, where each cell $(x,y)$ is associated with a scalar state $v_{x,y}$ representing the local vegetation quality. The backend linearly rescales the vegetation values obtained from the initial DEM/RGB pre-processing to a user-defined interval $[V_{min},V_{max}]$, with a default of $[-10,10]$. The resulting state’s values are interpreted as follows:

• $v_{x,y}<0$: Water-dominant cells, where lower values correlate with increased depth;
• $v_{x,y}\in [0,3)$: Bare terrain or soil, representing low-resource regions;
• $v_{x,y}\in [3,7)$: Grasslands or low-density vegetation;
• $v_{x,y}\in [7,10]$: High-density vegetation (e.g., bushes or trees).

$\mathcal{E}$ acts as the primary sensory input for the agents’ navigation policies and, conversely, is modified by their constructive and destructive behaviors.

2.2.1. Vegetation Dynamics and Seasonal Mean-Reversion

Vegetation regeneration is modeled as a spatially-distributed, seasonally-modulated stochastic process governed by a closed-loop restorative feedback mechanism. At each simulation step, the state $v_{x,y}$ of a cell $(x,y)\in \mathcal{E}$ is updated according to:

$$v_{x,y}(t+\Delta t)=\mathrm{clip}\left(v_{x,y}\left(t\right)+\Delta v_{x,y}\left(t\right),v_{min},v_{max}\right),$$

(2)

where the absolute vegetation state is strictly bounded to the interval $[v_{min},v_{max}]$ to prevent collapse or infinite overgrowth. The increment $\Delta v_{x,y}\left(t\right)$ is sampled for each cell within the valid growth interval $[a,b]$ via the indicator function ${\mathbb{I}}_{[a,b]}$:

$$\Delta v_{x,y}\left(t\right)\sim \mathcal{N}\left({\mu }_{\mathrm{eff}}\left(t\right),{\sigma }_v^2(x,y)\right)·{\mathbb{I}}_{[a,b]}\left(v_{x,y}\left(t\right)\right).$$

(3)

To capture the natural variety of the landscape, the plant growth variance (${\sigma }_v^2(x,y)$) varies with location. We first define a baseline variability ${\sigma }_{\mathrm{base}}={\sigma }_0\left|\gamma \right|$, where ${\sigma }_0$ is a configurable noise parameter and $\gamma$ is the baseline growth rate. This baseline is then amplified locally: growth becomes more variable in rough patches where a cell’s vegetation contrasts heavily with its surroundings. Computationally, we balance out these random changes across the map at each time step (zero-centering) so that the total amount of vegetation does not accidentally drift over time. This ensures that the macroscopic health of the landscape is controlled solely by the effective mean growth rate, ${\mu }_{\mathrm{eff}}\left(t\right)$:

$${\mu }_{\mathrm{eff}}\left(t\right)=\mathrm{clip}\left({\mu }_{\mathrm{base}}\left(t\right)+k\left(v_{\mathrm{target}}\left(t\right)-\overline{v}\left(t\right)\right),{\mu }_{min},{\mu }_{max}\right).$$

(4)

Here, $\overline{v}\left(t\right)$ is the spatial average of the vegetation across the growth mask. The parameter k acts as the proportional gain, creating a feedback that pulls the landscape’s average toward a moving seasonal target, $v_{\mathrm{target}}\left(t\right)$. The growth bounds are defined relative to the base rate $\gamma$ as ${\mu }_{min}=-\gamma$ and ${\mu }_{max}=2\gamma$, allowing the restorative feedback to temporarily exceed standard vernal growth when recovering from severe depletion. The target state $v_{\mathrm{target}}\left(t\right)$ is a continuous linear interpolation between the landscape’s winter minimum and its summer peak. In this framework, $v_{\mathrm{winter}}$ and $v_{\mathrm{summer}}$ are defined as $\pm 20\%$ of the environment’s initial spatial mean, anchoring the feedback loop to the site’s original carrying capacity. Both ${\mu }_{\mathrm{base}}\left(t\right)$ and $v_{\mathrm{target}}\left(t\right)$ are modulated by a normalized seasonal drive $S\left(t\right)$, ensuring that the restorative feedback tracks the vegetation cycle:

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{base}}\left(t\right)& =\gamma \left(2S\left(t\right)-1\right),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill v_{\mathrm{target}}\left(t\right)& =v_{\mathrm{winter}}+(v_{\mathrm{summer}}-v_{\mathrm{winter}})S\left(t\right).\hfill \end{array}$$

(6)

To capture the asymmetry of rapid vernal growth and gradual autumnal decay, $S\left(t\right)$ is defined using the probability density function (PDF) of a skew-normal distribution. Let $\tau \in [0,364]$ represent the day of the year. The seasonal drive is normalized to its peak as:

$$S\left(t\right)={\displaystyle \frac{f_{\mathrm{SN}}(\tau ;\xi ,\omega ,\alpha )}{f_{\mathrm{SN}}(\xi ;\xi ,\omega ,\alpha )}},$$

(7)

where $f_{\mathrm{SN}}$ is the skew-normal PDF:

$$f_{\mathrm{SN}}(\tau ;\xi ,\omega ,\alpha )={\displaystyle \frac{2}{\omega }}\varphi \left({\displaystyle \frac{\tau -\xi }{\omega }}\right)\Phi \left(\alpha {\displaystyle \frac{\tau -\xi }{\omega }}\right).$$

(8)

In this formulation, $\varphi$ and $\Phi$ are the probability density and cumulative distribution functions of the standard normal distribution, respectively. The parameters define the site’s biological rhythm: $\xi$ is the peak growth day, $\omega$ controls the seasonal window width, and $\alpha >0$ defines the right-skewness necessary to model slow autumnal decay. This model ensures that while individual cells exhibit realistic, contrast-driven local stochasticity, the macroscopic landscape preserves a stable, seasonally oscillating topological memory.

2.2.2. Acadia National Park Case Study: Baseline Vegetation Dynamics Validation

Environmental dynamics were validated on the Acadia dataset. The environment was initialized with the DEM-derived matrix from Figure 2(e), scaled to $[-10,10]$.

The baseline vegetation parameters were calibrated to reflect a seasonal cycle characterized by high local stochasticity and a stable global carrying capacity. Table 1 summarizes the mapping between the parameters implemented in the simulation backend and the formal coefficients from Section 2.2.1. The skew-normal distribution’s parameters are hard-coded to $\xi =90$ (peak growth at day 90), $\omega =50$ (seasonal window of 50 days), and $\alpha =5$ (strong right skew).

For this specific validation, we set a grass growth rate $\gamma =5\times {10}^{-4}{\mathrm{h}}^{-1}$, restricted to the growth interval $[0,8]$. To simulate localized micro-climatic variability, the stochastic scaling factor was set to ${\sigma }_0=0.5$. The global stabilizing proportional gain (mean-reversion strength) was set to $k=5\times {10}^{-4}$. The simulation was run for 1440 days ($\sim 4$ years) without beaver agents to observe the unforced system dynamics. The temporal evolution of the landscape state is shown in Figure 4. The system exhibits a bounded seasonal oscillation. In the bottom panel, the growth rate (blue dashed line) is driven by the skew-normal seasonal PDF, exhibiting rapid vernal acceleration, a saturated plateau, and gradual autumnal decay. The simulation shows a phase shift: the maximum growth rate plateau initiates near day 100, while the macroscopic vegetation (top panel) integrates this growth to peak near day 160. This cyclic behavior repeats reliably year after year. Furthermore, despite the injection of high local stochasticity (visualized here via the 10th–90th percentile band), the mean-reversion feedback prevents state divergence. Notably, the vegetation mean (oscillating between $\sim 1.0$ and $\sim 2.5$) consistently tracks higher than the median, which periodically returns to near zero during winter. This positive skew indicates a highly heterogeneous spatial distribution, confirming that the landscape is dominated by expansive barren with highly dense, localized vegetation peaks.

While these macroscopic temporal dynamics remain consistent across different landscapes, the spatial manifestation of the vegetation state $v_{x,y}$ depends on the quality of the dataset. Figure 5 shows the unforced spatial evolution of the environmental matrix $\mathcal{E}$ over a one-year simulation cycle for both the Acadia (DEM-based) and Rumney Ranch A2 (RGB-based) sites. The snapshots are taken every 30 days.

In the Acadia site (Figure 5(a)), the baseline vegetation proxy is derived from a rasterized elevation model. Thus, the resulting spatial field exhibits low-frequency, smooth gradients. Conversely, the Rumney Ranch A2 site (Figure 5(b)) illustrates the distinctly higher spatial granularity typical of RGB-derived environments. The simulation considers the baseline (see Section 2.1.2) as the initial condition, and captures variations in the underlying soil that elevation models inherently miss.

This structural difference validates the proposed dual-pipeline architecture. From a multi-agent perspective, the DEM-based environments present agents with macroscopic navigation gradients. In contrast, the RGB-based environments expose agents to a highly textured landscape, which is better suited for evaluating trailing and foraging behaviors.

Figure 5. Unforced spatial evolution of the environmental matrix $\mathcal{E}$ over one-year simulation. (a) The DEM-based Acadia site (b) The RGB-based Rumney Ranch A2 site.

Figure 5. Unforced spatial evolution of the environmental matrix $\mathcal{E}$ over one-year simulation. (a) The DEM-based Acadia site (b) The RGB-based Rumney Ranch A2 site.

2.2.3. Hydrological Dynamics and Channel Deepening

Complementing the regenerative vegetation cycle, the aquatic regions of the environment undergo a deterministic, passive-dynamic process that simulates continuous hydrological erosion. Within these defined channels, the state of water-cells is updated at each simulation step according to:

$$v_{x,y}(t+\Delta t)=max\left(v_{x,y}\left(t\right)-{\nu }_r·{\mathbb{I}}_{[c,d]}\left(v_{x,y}\left(t\right)\right),c\right),$$

(9)

where ${\nu }_r$ represents a constant channel deepening rate, and the indicator function ${\mathbb{I}}_{[c,d]}$ restricts this erosive process exclusively to the active hydrological interval $[c,d]$. To prevent unbounded deformation of the landscape, the state is clipped at the absolute depth boundary c. Table 2 summarizes these variables and their formal state-space coefficients.

Unlike the highly stochastic vegetation regeneration, this channel deepening establishes a continuous process that gradually steepens the streams over time. As will be detailed in the agent behavioral sections, the current implementation achieves canalization without relying on dam-building. Instead, the agents interact with the streams through stigmergic trailing during foraging. This localized, agent-induced erosion continuously degrades the river banks, facilitating water ingress that progressively fractures the adjacent terrain into auxiliary stream networks and ponds, expanding the colony’s habitable area.

2.2.4. Stigmergic Ledger and Topological Constraints

To support multi-agent decentralized coordination, the environment module maintains a secondary spatial ledger parallel to the primary state matrix $\mathcal{E}$. This visitation matrix, denoted as $\mathcal{V}\in {\mathbb{N}}^{N\times M}$, tracks a slowly decaying traversal history of all agents across the grid. $\mathcal{V}$ serves as the basis for the stigmergic trailing, enabling agents to sense and reinforce established foraging corridors without requiring direct communication.

Furthermore, the agent-induced canalization is governed by a topological boundary condition enforced by the environment. As we will better describe in Section 2.3, when an agent harvests resources from a cell, it subtracts a discrete quantity of vegetation $\Delta v$ from the local state $v_{x,y}$. However, to simulate contiguous water ingress rather than the spontaneous formation of isolated puddles, the environment dictates that an agent can only transition a terrestrial cell into a water cell ($v_{x,y}<0$) if the target cell lies within the immediate Moore neighborhood (D8) of an existing water body. This localized constraint serves as a physical law of the simulation grid, ensuring that agents’ destructive actions naturally lead to the continuous expansion of existing river networks.

2.3. Biomimetic Agent Architecture and Control Policies

The multi-agent system is modeled as a decentralized swarm of hybrid dynamical entities operating within the environmental matrices. To bridge the gap between biological observation and computational modeling, the agent architecture is grounded in two fundamental abstractions regarding energetic costs and spatial awareness.

Assumption A1

(Backpack and Energetic Proxy). Beavers transport collected materials using their mouths and forepaws, effectively limiting the volume they can move in a single trip. In the simulation, this is abstracted as a scalar, $L_{max}$, representing the carrying capacity. This parameter also serves as a proxy for the energetic cost of foraging; a higher $L_{max}$ indicates a more energetic agent capable of longer excursions before returning to a storage site.

Assumption A2

(Localized Mapping). Although natural beavers require time to learn new environments, the simulations herein operate on multi-year time scales. Thus, while agents have access to their global coordinate position, all behavioral and harvesting decisions are strictly constrained to a localized sensory horizon (e.g., a 5-meter radius).

2.3.1. Agent State Representation and Dynamics

Building upon these principles, the comprehensive internal state of the i-th agent (robot_beavers_backend.py) at discrete time step t is defined by the tuple:

$${\mathcal{S}}_i\left(t\right)=\langle {\mathbf{p}}_i\left(t\right),{\mathbf{v}}_i\left(t\right),L_i\left(t\right),{\mathcal{M}}_i^{env}\left(t\right),{\mathcal{M}}_i^{vis}\left(t\right)\rangle$$

(10)

where ${\mathbf{p}}_i\left(t\right)\in {\mathbb{R}}^2$ and ${\mathbf{v}}_i\left(t\right)\in {\mathbb{R}}^2$ represent the continuous 2D position and velocity vectors. The scalar $L_i\left(t\right)\in [0,L_{max}\left(t\right)]$ tracks the current harvested biomass load. The matrices ${\mathcal{M}}_i^{env}\left(t\right)\in {\mathbb{R}}^{N_{env}\times N_{env}}\subseteq \mathcal{E}$ and ${\mathcal{M}}_i^{vis}\left(t\right)\in {\mathbb{R}}^{N_{vis}\times N_{vis}}\subseteq \mathcal{V}$ represent the agent’s localized perceptions of the global environment $\mathcal{E}$ and the stigmergic visitation ledger $\mathcal{V}$, respectively.

On both x- and y-axes, the movement of the agent is governed by a decoupled double-integrator model, incorporating inertial mass m and viscous friction $\mu$. For a single spatial dimension, letting ${\mathbf{z}}_t={[p_t,v_t]}^{\top }$, the discrete-time state-space evolution is formulated as:

$${\mathbf{z}}_{t+\Delta t}=A{\mathbf{z}}_t+B{u}_t,$$

(11)

where $u_t$ is the control force derived from the navigation policy, and the matrices are:

$$A=\left[\begin{array}{cc}1& \Delta t\\ 0& 1-{\displaystyle \frac{\mu }{m}}\end{array}\right],B=\left[\begin{array}{c}0\\ {\displaystyle \frac{\Delta t}{m}}\end{array}\right].$$

(12)

While the inertial mass dampens the transient acceleration phase during movement, the linear dynamics are augmented with a non-linear hybrid stopping condition: if the control input is zero ($u_t=0$), the velocity state is instantaneously truncated to zero ($v_{t+\Delta t}=0$). This hybrid formulation eliminates inertial drift, allowing the agent to execute precise stigmergic actions (such as excavation) without overshooting the intended coordinate.

2.3.2. Hierarchical Cognitive Architecture: Tasks and Actions

The agent’s decision-making is structured as a two-layered hierarchical state machine. This decouples the overarching biological objectives (the Strategic Layer, see Figure 6) from their discrete physical execution within the simulated environment (the Tactical Layer).

Strategic Layer: The Task FSM

The overarching behavioral intent is governed by an atomic FSM evaluated at every temporal increment $\Delta t$. Because the agent relies on the localized spatial patch ${\mathcal{M}}_i^{env}\left(t\right)$ and on the stigmergic visitation ledger ${\mathcal{M}}_i^{vis}\left(t\right)$, it operates reactively without long-horizon path planning. The deterministic state transitions are evaluated sequentially:

• Store: Triggered if the agent’s capacity is saturated ($L_i\left(t\right)\ge \lfloor L_{max}\left(t\right)\rfloor$), or if adaptive parameter decay forces the maximum capacity below $10\%$ of its initial value $L_{max}^{\mathrm{init}}$, compelling a return to base.
• Harvest: Triggered if capacity is available, local resources are valid ($v_{x,y}\in \mathcal{H}\subseteq [V_{min},V_{max}]$), and the agent is not on a lodge site. To introduce natural behavioral variance, this state is gated by a stochastic $ϵ$-greedy condition ($\mathrm{rand}<{ϵ}_{\mathrm{greedy}}$).
• Explore: The default spatial routing state is initialized when immediate resources are depleted, invalid, or stochastically bypassed.

Tactical Layer: The Action FSM

Once a strategic task is assigned, the Tactical Layer dictates its computational execution. The overarching tasks are decomposed into three actions: move, remove_vegetation, or idle. Specifically, an EXPLORE task triggers a move action, whereas a HARVEST task triggers a remove_vegetation action upon reaching the target coordinate.

While the environmental framework is architecturally capable of supporting natural circadian rhythms (e.g., forcing a global idle state during nocturnal phases), the current agent implementation intentionally abstracts this cycle. Because nocturnal resting phases effectively act as a linear temporal scaling factor, their inclusion would introduce redundant computational complexity without fundamentally altering the emergent spatial topology. Thus, the idle action is strictly for transient inactive phases between task evaluations.

2.3.3. Strategic Target Selection and Agent Roles

Once the strategic layer triggers an EXPLORE or HARVEST task, the agent must identify a specific spatial waypoint ${\mathbf{p}}^{*}$ within its localized sensory patch ${\mathcal{M}}_i^{env}\left(t\right)$. To mimic the variability of biological foraging, the framework introduces two distinct agent roles, Explorers and Builders, which alter how an agent evaluates the spatial utility $U(x,y)$ of neighboring cells. The framework also provides a user-configurable parameter (exploration_map) to toggle stigmergic coupling. Let’s recall the global stigmergic ledger $\mathcal{V}\in {\mathbb{R}}^{N\times M}$. To ensure that the magnitude of the stigmergic record remains comparable to the vegetation state during utility calculations, the environment continuously normalizes the entire ledger into a bounded matrix $\tau \in [0,V_{max}]$, as specified in eq. (19) below.

• Mode 1: Resource-Driven (vegetation_quality): Navigation relies exclusively on the primary environmental state. Explorers seek pristine foraging grounds by scaling utility directly with resource density: $U_{\exp }\propto {(\epsilon +v_{x,y})}^2$. Conversely, Builders seek construction sites (e.g., bare terrain or water boundaries) by mathematically inverting the resource gradient: $U_{\mathrm{bld}}\propto {(\epsilon +v_{x,y})}^{-2}$.
• Mode 2: Stigmergically Coupled (vegetation_visits): Navigation integrates the secondary spatial ledger to reinforce established corridors. For Explorers, the resource utility is amplified by historical traffic: $U_{\exp }\propto (\epsilon +{\tau }_{x,y})(\epsilon +v_{x,y})$. For Builders, the inverted gradient is multiplied by the visitation density: $U_{\mathrm{bld}}\propto (\epsilon +{\tau }_{x,y})/(\epsilon +v_{x,y})$.

Remark 1

(Stigmergic Coordination in Target Selection). The stigmergically coupled evaluation mode marks the first active integration of stigmergy into the agent’s decision-making pipeline. In this mode, builders do not merely seek arbitrary bare terrain; the visitation term $(\epsilon +{\tau }_{x,y})$ actively guides their exploratory efforts along highly traversed corridors.

2.3.4. Tactical Execution: Dynamic Control and Obstacle Avoidance

After computing the spatial utility for all valid cells in the immediate neighborhood, the agent calculates the local utility gradient $\Delta U_j$. The probability of selecting a specific neighboring coordinate i as the next waypoint ${\mathbf{p}}^{*}$ is governed by a Softmax distribution:

$$P\left(i\right)={\displaystyle \frac{exp(\Delta U_i)}{{\sum }_{j}exp(\Delta U_j)}}.$$

(13)

When the tactical layer authorizes a move action, the agent must physically navigate to the sampled waypoint ${\mathbf{p}}^{*}$. The control effort $u_t$ driving the double-integrator dynamics (eq. (11)) is generated on each axis by a closed-loop Artificial Potential Field (APF) controller, augmented with a Proportional-Integral-Derivative (PID) goal-seeking law:

$$u_t=\underset{\mathrm{Attractive}\mathrm{Force}}{\underbrace{\left(K_p{\mathbf{e}}_t+K_i\int {\mathbf{e}}_{t}dt+K_d{\displaystyle \frac{d{\mathbf{e}}_t}{dt}}\right)}}-\underset{\mathrm{Repulsive}\mathrm{Force}}{\underbrace{\beta \nabla {\mathcal{M}}_{\mathrm{repulsive}}\left({\mathbf{p}}_t\right)}}+\underset{\mathrm{Water}\mathrm{Drift}}{\underbrace{\gamma \widehat{\mathbf{f}}\left({\mathbf{p}}_t\right)}},$$

(14)

where ${\mathbf{e}}_t={\mathbf{p}}^{*}-{\mathbf{p}}_t$ is the positional error vector, $\beta$ is a scaling weight applied to the environmental penalty gradient, and $\gamma$ represents the strength of stream currents. To accurately reflect terrestrial traversal costs, the repulsive field ${\mathcal{M}}_{\mathrm{repulsive}}$ is dynamically constructed based on the agent’s immediate local neighborhood. Similar to the exploration utility, this penalty supports two configurable modes:

• Resource-Driven Friction (vegetation_quality): Evaluated as ${\mathcal{M}}_{\mathrm{repulsive}}\propto v_{x,y}$. Dense, high-cost patches naturally generate strong localized repulsive gradients, forcing the trajectory to dynamically weave around barriers.
• Stigmergically Mitigated Friction (vegetation_visits): Evaluated as ${\mathcal{M}}_{\mathrm{repulsive}}\propto v_{x,y}/(1+{\tau }_{x,y})$. Historical traversal acts as a physical modifier. Highly visited cells exhibit lowered resistance, simulating trampled vegetation and foraging trails.

Methodological Note: To maintain clarity and isolate emergent behaviors, researchers are advised to apply stigmergic coupling exclusively to either the high-level exploration heuristic (Section 2.3.3) or this low-level APF controller, avoiding redundant systemic feedback.

2.3.5. Hydrological Drift and Riparian Biasing

Distinct from the cognitive obstacle avoidance logic governed by ${\mathcal{M}}_{\mathrm{repulsive}}$, the simulation subjects the agent to direct physical hydrological forces. This is represented by the third term in eq. (14), where $\widehat{\mathbf{f}}$ is an externally defined river flow vector acting upon the aquatic cells of the environment ($v_{x,y}<0$). Because this flow acts as an independent environmental drift, it is not bound by the repulsive weighting factor ($\beta$). Consequently, when the agent navigates aquatic cells, the flow vector is superimposed directly onto the final control effort $u_t$. This current assists downstream and resists upstream movement.

2.3.6. Distance-Gated Local Minima Resolution

A well-documented flaw of standard APF controllers is the tendency for agents to stall in local spatial minima, particularly when a target is adjacent to a strong repulsive obstacle. To resolve this without incurring the computational overhead of global path-planning algorithms (e.g., A* or RRT), the framework uses a distance-gated terminal attraction phase. The repulsive field is strictly evaluated only when the spatial error norm exceeds a critical threshold (e.g., ${∥{\mathbf{e}}_t∥}_2>4.0$ pixels). Once the agent breaches this proximity radius, $\nabla {\mathcal{M}}_{\mathrm{repulsive}}$ is instantaneously forced to zero. The controller drops the obstacle-avoidance constraint and relies solely on the purely attractive PID law to converge to the final target. This effectively allows the agent to "push through" boundary friction to reach its localized objective, mimicking a biological agent finalizing its approach to a specific resource.

2.3.7. Environmental Modification

When a HARVEST task gets to its spatial target, the tactical layer transitions the agent from locomotion to excavation by executing the remove_vegetation action. This process facilitates a direct state transfer between the environmental matrix $\mathcal{E}$ and the agent’s internal resource capacity. At the target coordinate $(x,y)$, the agent extracts a discrete quantum of biomass $\Delta v$, synchronously updating both the environmental state and its internal load:

$$\begin{array}{cc}\hfill v_{x,y}(t+\Delta t)& =v_{x,y}\left(t\right)-\Delta v\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill L_i(t+\Delta t)& =L_i\left(t\right)+\Delta v\hfill \end{array}$$

(16)

As introduced in Section 2.2.4, the environment module enforces a topological constraint during the action: while an agent can freely deplete terrestrial vegetation ($v_{x,y}>0$), it is prohibited from transitioning a terrestrial cell into an aquatic cell ($v_{x,y}<0$) unless the target coordinate lies strictly within the immediate Moore neighborhood (D8) of an existing water body. This constraint allows the agents’ foraging actions to organically fracture the terrain, driving the contiguous expansion of the existing streams into the surroundings.

2.3.8. Writing to the Stigmergic Ledger

The stigmergic trail mechanism is governed by a three-step update cycle at each time step t: (1) agent deposition, (2) environmental decay, and (3) global normalization.

1.

In parallel with physical excavation, each agent navigating the grid deposits a discrete traffic marker $\Delta \tau$ onto the ledger at its current coordinate $(x,y)$:

$${\mathcal{V}}_{x,y}\leftarrow {\mathcal{V}}_{x,y}+\Delta \tau$$

(17)

2.

The ledger values decay, to simulate the fading of trails:

$$\mathcal{V}\leftarrow \xi ·\mathcal{V},\mathrm{with}\xi \in (0,1)$$

(18)

3.

Finally, the entire ledger is normalized to the vegetation scale to yield the local visitation record ${\tau }_{x,y}(t+\Delta t)$:

$${\tau }_{x,y}(t+\Delta t)={\displaystyle \frac{{\mathcal{V}}_{x,y}-min\left(\mathcal{V}\right)}{max\left(\mathcal{V}\right)}}·V_{max}$$

(19)

Once processed through this cycle, the normalized record ${\tau }_{x,y}$ reduces localized friction for subsequent agents (as formalized in Section 2.3.4) and actively shapes the colony’s probabilistic exploration gradients (Section 2.3.3), closing the decentralized feedback loop.

2.3.9. Adaptive Foraging and Parameter Decay

Biological systems rarely operate under static conditions; prolonged failure to acquire resources can induce behavioral shifts. To capture this reality, the architecture implements a time-dependent parameter-decay loop driven by the agent’s foraging efficiency. If the frequency of successful harvest events falls below a critical threshold, the agent registers a stagnation state. To resolve this, the framework applies a multiplicative decay to a pair of configurable behavioral parameters: the harvest interval ($\mathcal{H}$) and the maximum carrying capacity ($L_{max}$). Each decaying parameter induces a specific tactical shift in the agent’s logic and its resulting impact on the simulated landscape:

• Harvest Interval ($\mathcal{H}$): As stagnation persists, the acceptable bounds of vegetation density systematically expand. Behaviorally, the agent lowers its foraging standards out of simulated desperation. Environmentally, this accelerates the depletion of marginal canopy zones that the swarm would normally bypass.
• Carrying Capacity ($L_{max}$): A gradual reduction in $L_{max}$ serves as a biological proxy for metabolic exhaustion. A shrinking capacity lowers the threshold required to trigger the STORE task (Section 2.3.2), terminating excursions and compelling a premature return to the lodge to prevent starvation.

Upon a successful sequence of harvest actions, the agent breaks the stagnation cycle, and all three parameters are instantaneously reset to their initial baseline values. To facilitate reproducibility and establish a parameter space for subsequent experimental validation, the core user-definable variables governing the multi-agent behavioral logic, dynamics, and environmental dynamics are reported in Table 3, Table 4, and Table 5, respectively.

3. Results and Discussion

To experimentally validate the proposed framework, the agent architecture was deployed across two distinct environmental matrices to isolate specific behaviors. First, the Rumney Marsh dataset was utilized to evaluate the low-level kinematic control and hydrological drift mechanics, leveraging its well-defined stream networks. Second, the Rumney Ranch dataset was selected to investigate macroscopic swarm dynamics and emergent stigmergic trail formation. This latter site provides a useful spatial balance; it’s an appropriate size to support long-term simulations without compromising computational tractability, and it also features high-resolution RGB imagery, establishing a baseline for future cross-validation of simulated trajectories against actual biological trails.

3.1. Agent Dynamics and Navigation Policies

This subsection evaluates the agent’s movement capabilities on the Rumney Marsh dataset, focusing on the integration of the double-integrator physics engine with the distance-gated APF controller. First, we establish how the agent’s velocity scales dynamically with both its cognitive spatial horizon and the chosen integration time. We then analyze the necessary tradeoff between computational efficiency and biomimetic accuracy.

3.1.1. Experimental Setup

To isolate the routing dynamics from the multi-agent coordination, a single-agent environment ($n=1$) was considered. Active biomass removal ($\Delta v=0$) and seasonal vegetation growth were disabled, forcing the agent to navigate a static spatial utility field $U(x,y)$. Similarly, adaptive parameter decay was deactivated to ensure the agent’s cognitive thresholds remained constant throughout the simulation. Furthermore, stigmergic coupling was explicitly disabled (exploration_map = vegetation_quality). This setup isolated the agent’s response to primary environmental gradients and prevented "stigmergic self-attraction," namely, an artifact in which an isolated agent becomes trapped in self-reinforcing orbital loops by continuously following its immediate traversal history. Furthermore, the temporal decay of the visitation ledger was set to zero to preserve a cumulative spatial memory of the entire trajectory, as trail fading becomes a relevant control dynamic only during multi-agent stigmergic coordination. Finally, the river flow was set to zero to eliminate the influence of passive hydrological transport.

To govern the double-integrator dynamics, the agent’s physical parameters were established to reflect the empirical biomechanics and environmental resistance experienced by Castor canadensis. Specifically, the inertial mass was set to $m=20.0$ kg, coupled with a viscous friction coefficient of $\mu =3.0$. The continuous control effort driving this physical model was regulated by the closed-loop APF controller (eq. (14)). To ensure numerical stability and realistic decision-making, the system’s state-space integration time step was decoupled to $dt=0.1$ hours (equivalent to 6-minute cognitive intervals).

The controller’s proportional, derivative, and integral gains were dynamically tuned to $K_p=200.0$, $K_d=300.0$, and $K_i={10}^{-5}$, respectively. During traversal, terrain negotiation was dictated by a localized repulsive weight ($\beta =0.2$), calculated over a two-cell spatial neighborhood and based on the underlying vegetation quality matrix. Finally, to counteract the asymptotic slowdown inherent to standard APF algorithms near target destinations, a terminal positional accuracy threshold of $0.5$m was enforced.

3.1.2. Cognitive Spatial Horizons and Velocity Modulation

This section evaluates the impact of the agent’s cognitive spatial horizon on its macroscopic velocity. To isolate the dynamics’ impact on the agent’s decision-making cycle, the action-selection mechanism was constrained using an $ϵ$-greedy policy with $ϵ=0.5$.

This configuration enforces a symmetric probability distribution between exploration and harvesting. Because the harvest state acts as a stationary temporal delay in this specific experimental context ($\Delta v=0$), the agent is forced into a stop-and-go pattern, spending approximately 50% of its time in an idle configuration. The agent’s ability to explore the surrounding area is governed by the radial depth of its spatial horizon, i.e., $N_{vis}$.

To first establish the physics engine’s mechanistic capabilities, a baseline simulation was conducted with a theoretical cognitive spatial horizon of 50 m radially (100 m total diameter). The kinematic impact of this extended horizon is captured in Figure 7(a) and Figure 7(c). The control error norm (Figure 7(a)) shows that destination updates frequently project targets up to 50 meters away, with the double-integrator converging on the 0.50 m accuracy threshold. Notably, minor transient oscillations are observable during this terminal convergence. This behavior is a deliberate consequence of the controller’s tuning; an aggressive PID configuration was explicitly selected to maximize acceleration, ensuring a rigorous assessment of the system’s absolute high-speed traversal behavior. Granted sufficient spatial runway by this expanded horizon, the physical engine reaches peak transit velocities exceeding 150 m/h (Figure 7(c)). Due to the $ϵ$-greedy policy, the agent’s effective macroscopic velocity averages $\approx 11$ m/h.

3.1.3. Kinematic Profiles of Localized Terrestrial Foraging

To transition from theoretical to biomimetic traversal, the spatial horizon was subsequently restricted to 5 m radially (10 m total diameter). This heavily constrained perimeter reflects the localized sensory capabilities of C. canadensis during terrestrial motion. The kinematic bottleneck induced by this restricted decision horizon is demonstrated in Figure 7(b) and Figure 7(d). Because the physical engine is confined to short-distance transits (never exceeding 6 meters, Figure 7(b)), the agent remains within the transient acceleration and deceleration phases of the APF controller. Consequently, the agent never achieves its biomechanical top speed; velocity peaks are capped at approximately 20 to 25 m/h (Figure 7(d)). When combined with the 50% idling from the $ϵ$-greedy policy, this kinematic truncation restricts the agent’s macroscopic displacement rate to ≈ 3 m/h. In nature, free-ranging beavers achieve a higher overall displacement rate (≈ 536 m/h) by leveraging efficient aquatic routes during a compact ≈ 9.7-hour nocturnal window [47].

Because our framework explicitly models a land-only, continuous 24-hour operational cycle without rest, spreading the daily foraging footprint across a full, nonstop diurnal loop dilutes the hourly velocity. While this 3 m/h speed is highly conservative compared to true biological benchmarks, it establishes a safe, numerically feasible lower bound for the simulation while still capturing the cautious terrestrial locomotion and localized patch exploitation typical of foraging C. canadensis.

3.2. Computational Tractability and Temporal Upscaling

While the integration interval $dt=0.1$ h replicates the micro-kinematics of localized foraging, it introduces a computational bottleneck. A single agent requires 7 seconds to simulate a 24-hour period. When extrapolated to multi-agent swarms operating over annual timescales, this overhead becomes prohibitive. To facilitate large-scale simulations, the framework needs temporal upscaling, increasing the integration step to $dt=1.0$ h. However, increasing the integration step alters the stability of the discrete-time engine. Maintaining the baseline velocity ($v\approx 11$ m/h) at $dt=1.0$ h results in discrete positional updates of 11 m per tick. This exceeds the agent’s sensory radius ($r_{max}=5$ m), causing it to overshoot local waypoints. To respect spatial motion boundaries, the allowable velocity must satisfy $v_{max}\le r_{max}/dt$, capping transit speeds at 5 m/h.

Consequently, the physical dynamics and APF control gains were retuned to improve system damping. Mass and friction were adjusted to $m=5$ and $\mu =1$, paired with lower control gains ($K_p=0.5$, $K_d=2.0$, $K_i={10}^{-5}$). As illustrated in Figure 8, this recalibration restricts the peak velocity to around 2.2 m/h, well below the 5 m/h bound. The error norm verifies this stability, showing that spatial projections remain within the 5-meter horizon.

To evaluate the macroscopic impact of this dampening, the simulation duration was extended to 5 days, providing a comparative temporal scale against the baseline model. Because the transit time is extended, and the symmetric action-selection policy ($ϵ=0.5$) forces the agent into a 1-hour idle state during ≈50% of its decision time, the behavioral frequency is bottlenecked. As a consequence, the agent requires 5 days to generate the same number of destination updates achieved in a single day under the baseline ($dt=0.1$ h), transitioning from rapid micro-transit to extended patch exploitation.

Figure 8. Kinematic profile of the temporally upscaled framework ($dt=1.0$ h) over a 5-day simulation on the Rumney Marsh site. (a) The control error norm illustrates a reduction in destination updates, with spatial projections bounded within the 5-meter horizon. Coupled with the $ϵ$-greedy idling penalty ($ϵ=0.5$), it takes 120 hours to accumulate the same number of waypoint visits generated in 24 hours under the baseline (Figure 7(b)), confirming the shift to extended patch exploitation. (b) The velocity profile shows the effect of the dampened physical parameters ($m=5,\mu =1,K_p=0.5,K_d=2.0$), which cap the peak velocity at ∼2.2 m/h. This dampening satisfies the spatial motion boundary condition ($v_{max}\le 5$ m/h), preventing the agent from overshooting the cognitive horizon.

Figure 8. Kinematic profile of the temporally upscaled framework ($dt=1.0$ h) over a 5-day simulation on the Rumney Marsh site. (a) The control error norm illustrates a reduction in destination updates, with spatial projections bounded within the 5-meter horizon. Coupled with the $ϵ$-greedy idling penalty ($ϵ=0.5$), it takes 120 hours to accumulate the same number of waypoint visits generated in 24 hours under the baseline (Figure 7(b)), confirming the shift to extended patch exploitation. (b) The velocity profile shows the effect of the dampened physical parameters ($m=5,\mu =1,K_p=0.5,K_d=2.0$), which cap the peak velocity at ∼2.2 m/h. This dampening satisfies the spatial motion boundary condition ($v_{max}\le 5$ m/h), preventing the agent from overshooting the cognitive horizon.

3.3. Landscape-Scale Navigation and Topographic Routing

Having established the stability of the physical engine ($dt=1.0$ h), the macroscopic spatial routing was evaluated. It is critical to distinguish the scope of the APF from the high-level tactical layer. The APF repulsive gradient ($\nabla {\mathcal{M}}_{\mathrm{repulsive}}$) is intentionally near-sighted. It is designed to induce localized swaying to bypass immediate topographic resistance. It does not act as an exclusionary wall. If the tactical layer designates a waypoint deep within a dense thicket, the APF will minimize the traversal cost while ultimately permitting entry.

To isolate the APF’s pathfinding capabilities, a deterministic control simulation was run. The vegetation removal rate was disabled ($\Delta v=0$) to freeze the cost surface, and the agent was assigned a fixed, long-range destination across the environment. Figure 9 illustrates the emergent trajectories under varying control configurations. Under a standard goal-seeking PID controller (Figure 9(a)), the trajectory follows a biologically unnatural path directly to the destination. Conversely, activating the repulsive field (Figure 9(b)) fundamentally alters the trajectory. Without preventing the agent from reaching its terminal destination, the repulsive gradients deflect the trajectory, resulting in a localized weaving motion that circumnavigates the densest vegetation areas. Lastly, the framework was subjected to physical hydrological biasing to simulate riparian transit.

Utilizing the independent environmental drift vector ($\gamma \widehat{\mathbf{f}}$) defined in eq. (14), the agent was forced to cross a water body under two opposing current conditions. When subjected to a leftward drift vector $[-1,-1]$ (Figure 9(c)), the current displaces the agent off its optimal terrestrial vector, resulting in a pronounced leftward bow during the aquatic crossing. Once the agent exits the aquatic cells, the drift dissipates, and the trajectory resumes its standard terrestrial pathfinding. Reversing the current to induce a rightward drift $[1,1]$ (Figure 9(d)) yields the exact mirrored response, pulling the trajectory into a wide rightward arc. This deterministic test confirms that the framework successfully integrates independent physical environmental forces with cognitive terrain avoidance.

Figure 9. Deterministic spatial routing across a static landscape. Across all panels, the agent’s trajectory is traced in red, originating from the home lodge (green circle) and targeting the terminal destination (magenta circle). (a) A pure attractive PID controller forces a biologically unnatural path through high-cost obstacles. (b) Activation of the near-sighted APF repulsive field induces localized weaving, allowing the agent to circumnavigate vegetation without abandoning its terminal waypoint. (c, d) The introduction of an independent aquatic drift vector ($\gamma \widehat{\mathbf{f}}$) displaces the agent during water crossings, creating pronounced rightward and leftward bows before terrestrial weaving resumes.

Figure 9. Deterministic spatial routing across a static landscape. Across all panels, the agent’s trajectory is traced in red, originating from the home lodge (green circle) and targeting the terminal destination (magenta circle). (a) A pure attractive PID controller forces a biologically unnatural path through high-cost obstacles. (b) Activation of the near-sighted APF repulsive field induces localized weaving, allowing the agent to circumnavigate vegetation without abandoning its terminal waypoint. (c, d) The introduction of an independent aquatic drift vector ($\gamma \widehat{\mathbf{f}}$) displaces the agent during water crossings, creating pronounced rightward and leftward bows before terrestrial weaving resumes.

3.4. Swarm Dynamics and Stigmeric Trail Formation

Having established the APF control’s routing capabilities, this subsection evaluates multi-agent swarm dynamics over extended temporal horizons. To facilitate this macroscopic analysis, the simulation environment is transitioned to the Rumney Ranch dataset. This specific site was selected for its availability of RGB imagery, which provides a baseline for future comparisons of emergent simulated paths with actual biological trails. Specifically, we investigate colony-level spatial organization by contrasting the resource-driven friction map ($\mathtt{vegetation}\_\mathtt{quality}$) against the stigmergic visitation ledger ($\mathtt{vegetation}\_\mathtt{visits}$). We evaluate how physical trampling canalizes swarm locomotion and how this emergent structure is modulated by the maximum payload ($L_{max}$).  

3.4.1. Experimental Setup

To evaluate long-term spatial dynamics, a multi-agent colony ($n=10$) was deployed over a 6-month period. To establish a controlled baseline, seasonal vegetation growth, hydrological drift, and the visitation decay rate were suspended, ensuring that the colony’s spatial footprint was permanently captured and isolating the impact of physical trampling on spatial routing. During locomotion, the visitation ledger was updated with $\Delta \tau =1$, and the vegetation removal rate was set to $\Delta v=1.0$, constrained by a harvesting interval of $\mathcal{H}=[2,6]$, enforcing a moderate patch exploitation, driving agents toward mid-density vegetation while avoiding thickets and barren ground. The maximum load was varied between restricted ($L_{max}=100$) and extended ($L_{max}=1000$) limits, with decay rates of $0.05$ and $0.005$ applied to the maximum load and harvest interval, respectively.

3.4.2. Spatial Diffusion and Lodge-Centered Foraging

Figure 10 contrasts the cumulative traversal footprint of the swarm under two distinct routing policies. The macroscopic divergence in these spatial topologies stems directly from the formulation of the high-level utility matrix, $U(x,y)$, which governs the tactical evaluation of the landscape. Crucially, the manifestation of these emergent structures is depends on the agent’s maximum load capacity ($L_{max}$). Because all foraging sorties originate from and return to the central lodge (see Figure 6), intense radial visitation in the immediate vicinity of the hub is obligate across all models. When $L_{max}$ is strictly constrained (Figure 10(a) and Figure 10(c)), agents must frequently interrupt their foraging to return home, resulting in a persistent trampling around the central lodge.

When the payload capacity is expanded to support long-range foraging expeditions (Figure 10(b) and Figure 10(d)), both configurations exhibit extended exploration and the formation of distinct trajectories leading to displaced harvest patches. Yet, a structural difference emerges in how these resources are exploited. When the framework relies exclusively on the static resource map (Figure 10(b)), the utility is strictly a function of resource density, i.e., $U(x,y)\propto v_{x,y}^2$. Lacking environmental memory, individual agents are blind to their peers’ routing choices. Consequently, while they successfully travel far from the lodge, the expansion front remains highly diffuse. This results in a fragmented spatial footprint, characterized by smaller, scattered harvest patches at the terminal ends of less pronounced trajectories (red boxes), as agents independently fan out across the landscape. Conversely, under stigmergic coordination (Figure 10(d)), the utility matrix integrates the continuous visitation ledger, i.e., $U(x,y)\propto v_{x,y}·{\tau }_{x,y}$. As pioneer agents push further into the landscape, their deposition of $\tau$ markers increases the utility of the cells they traverse. Subsequent agents systematically exploit these gradients, tending to converge on established routes rather than breaking new paths. Rather than fragmenting into scattered individual foraging sites, the harvesting consolidates into larger patches (red boxes). Because agents are drawn to areas where others have navigated, the overall visitation map remains more cohesive.

3.5. Foraging Efficiency and Spatial Stochasticity

While the maximum payload capacity $L_{max}$ governs the spatial range of the swarm, the actual realization of that range is governed by the agent’s cognitive heuristics. Biological entities rarely act as perfect mathematical optimizers; their locomotion incorporates spatial stochasticity driven by territorial patrolling, predator evasion, or sub-optimal decision-making. To evaluate this dynamic, in Section 2.3.2 we introduced an $ϵ$-greedy action-selection policy that determines the probability that an agent harvests local vegetation versus exploring the local surroundings.

To isolate the macroscopic effects of this heuristic, the swarm’s physical range was unconstrained ($L_{max}=1000$), and the colony was simulated under two opposing cognitive regimes over a 6-month integration period. The remainder of the configuration was held constant with respect to Section 3.4. When the swarm operates with a high exploitation probability ($ϵ\text{-}\mathrm{greedy}=0.9$, Figure 11(a)), the agents prioritize immediate energetic returns. Even with the capacity for long-range travel, the deterministic focus forces the swarm to systematically strip local vegetation patches before advancing. This results in a highly concentrated footprint centered around the lodge. Conversely, suppressing the deterministic focus to favor high spatial stochasticity ($ϵ\text{-}\mathrm{greedy}=0.1$, Figure 11(b)) drastically alters the macro-topology of the colony. The agents frequently bypass viable local resources to explore peripheral zones, resulting in a diffuse network of foraging trails.

3.6. Ecosystem Engineering: Foraging vs. Building Phenotypes

Having established the foundational mechanisms of stigmergic canalization and spatial stochasticity, this final experiment evaluates the framework’s capacity to simulate complex, multi-agent ecosystem engineering. To capture this macroscopic division of labor, we exploit the two distinct behavioral phenotypes, Explorers and Builders.

3.6.1. Detailed Experimental Setup

To simulate complex, multi-agent ecosystem engineering over an extended temporal horizon, a heterogeneous swarm ($n=4$, equally divided into 2 Explorers and 2 Builders) was deployed over a continuous 4-year integration period. The temporal resolution was fixed at $dt=1.0$ h, with agents initialized at a central home lodge located at coordinates $(60,60)$ within the Rumney Ranch dataset.

The environment was rendered fully dynamic to capture long-term ecological feedback through coupled vegetation and hydrological models. Vegetation quality v was bounded within $[0,10]$, with regrowth occurring in the interval $[0,8]$ at a base rate of $1\times {10}^{-4}$. This process incorporated stochasticity (${\sigma }_0=0.5$) and a mean-reversion strength of $0.0005$. Crucially, the simulation introduced dynamic hydrology to model water mechanics: water puddles within the interval $[-3,0]$ triggered river expansion with a velocity factor of $1\times {10}^{-4}$. The stream flow was defined by a directional vector $d=[1,1]$ and an intensity factor of $0.2$. To simulate the fading stigmergic trace of traversal markers over the multi-year cycle, a visitation decay factor of $0.99$ was applied. Low-level agent locomotion remained governed by the integrator dynamics detailed in Section 3.2.

The continuous visitation ledger was incremented by $\Delta \tau =1.0$. To simulate progressive environmental depletion, adaptive decay rates of $0.001$ were applied to both the agent’s payload capacity and the bounds on the harvest interval during active sorties. While sharing identical base kinematics and ecosystem impact parameters, the agents were bifurcated into two specialized ecological roles defined by distinct cognitive heuristics:

• The Explorer: Designed to simulate a foraging vanguard, these agents operated with a high payload capacity ($L_{max}=1000$) and a significant vegetation removal rate of $3.0$ per harvest, enabling wide-range spatial diffusion. They used an $ϵ$-greedy policy with $ϵ=0.8$, balancing stochastic search with resource exploitation. Their cognitive map was coupled exclusively to the static resource layer (vegetation_quality) and was constrained to target mature, high-density vegetation within the interval $\mathcal{H}=[6,9]$. This configuration drives uncoordinated exploration based on resource availability.
• The Builder: Functioning as localized environmental engineers, these agents were restricted to a minimal payload capacity ($L_{max}=10$) and a lower vegetation removal rate of $1.0$, enforcing a high-frequency, obligate proximity to the central lodge. They operated with a fully deterministic, exploitative policy ($ϵ=1.0$) and were constrained to target low-density or previously excavated terrain, namely within $\mathcal{H}=[0,3]$ values.
  Their routing was strictly coupled to the dynamic visitation ledger (vegetation_visits), leading them to reinforce heavily trafficked corridors. This design fosters a positive feedback loop where Builders preferentially modify areas of the landscape initially exploited by Explorers and subsequently reinforced by peer Builders.

3.7. Results and Discussion

The simulation results demonstrate the emergence of a stable, self-regulating feedback loop between agent behavior and environmental dynamics. By maintaining conservative vegetation regrowth parameters, we observed the system’s capacity to maintain long-term ecological equilibrium rather than rapid resource exhaustion.

3.7.1. Ecological Stability and Vegetation Dynamics

The analysis of the vegetation layer (Figure 12) confirms a resilient environmental equilibrium over the 4-year cycle. Despite continuous harvesting activity by both agent roles, the mean and median of the vegetation quality oscillate within a narrow band ($v\in [5.5,6.5]$). The "sawtooth" fluctuations in the growth rate (bottom panel, Figure 12) illustrate the regrowth mechanism. This pattern suggests that the system operates under a disturbed equilibrium governed by the mean-reversion dynamics; the environment is regulated by a continuous mean-reversion strength ($0.0005$) that draws the vegetation quality toward a defined annual mean. This mechanism prevents the landscape from reaching maximum saturation ($v=10$) independently of agent activity, while simultaneously ensuring recovery from localized depletion caused by the agents’ disturbance.

3.7.2. Emergent Spatial Patterns and Role Differentiation

The spatial evolution of the landscape (Figure 13) validates the specialized heuristics assigned to the Explorer and Builder roles. The snapshots reveal a distinct concentration of activity in the bottom-left quadrant, driven by proximity to the central lodge and by the specific density intervals ($\mathcal{H}$) targeted by the agents. The differentiation in agent footprints is statistically and visually significant:

• Explorers (Blue): Exhibit a diffuse spatial distribution. This reflects their higher payload capacity and stochastic $ϵ$-greedy policy, allowing them to act as a foraging vanguard that identifies and exploits high-density vegetation patches ($[6,9]$).
• Builders (Red): Demonstrate highly localized, high-density corridors. Due to their restricted payload, these agents are forced to make high-frequency return trips to the lodge. Their coupling to the vegetation_visits map creates a positive feedback loop, reinforcing paths initially opened by peer agents.

Figure 13. Snapshots of the 4-year evolution showing localized depletion and visitation patterns.

Figure 13. Snapshots of the 4-year evolution showing localized depletion and visitation patterns.

This succession confirms that the Builder role functions as an intensifier of existing spatial traces, while the Explorer provides the necessary spatial expansion.

3.7.3. Density-Dependent Geomorphological Change

To investigate the scalability of agent-driven environmental modification, a comparative simulation was conducted with a colony size of 20 agents (Figure 14), in which the increased agent density triggered significant geomorphological shifts.

As illustrated in Figure 14(b), the higher agent density facilitates a more aggressive exploitation of the terrain, crossing a threshold that activates the model’s dynamic hydrology. The increased frequency of excavation within the critical depth interval ($[-3,0]$) results in visible expansion of the river channel into a pond over the four-year cycle.

Figure 14. Comparative spatial evolution across different colony sizes. The bigger colony (bottom) demonstrates the transition from localized harvesting to large-scale geomorphological engineering.

Figure 14. Comparative spatial evolution across different colony sizes. The bigger colony (bottom) demonstrates the transition from localized harvesting to large-scale geomorphological engineering.

4. Conclusion and Future Directions

This paper presented a bio-inspired Multi-Agent System (MAS) framework to simulate emergent landscape morphology. By coupling a dynamic environment with bifurcated agents (Explorers and Builders), we demonstrated that beaver trail and canal networks arise from simple behavioral heuristics. Transitioning the framework into a quantitative predictive tool requires rigorous empirical validation. Accordingly, we now discuss current limitations and outline critical next steps.

4.1. Empirical Validation and Automated Feature Extraction

The primary bottleneck in validating simulated spatial networks against empirical data is the extraction of emergent environmental traces. As seen in the Rumney Ranch case study, BDC trail systems are often optically faint and heavily blended with surrounding vegetation (Figure 15(a)).

While standard preprocessing techniques (e.g., vegetation indices) effectively map general biomass, they struggle to autonomously isolate these narrow, high-frequency features (Figure 15(b)). Extracting a clean topological network currently requires manual annotation to enhance visual contrast (Figure 15(c) and Figure 15(d)). Although this yields the high-fidelity structural data necessary to confirm network presence, it is inherently unscalable and precludes landscape-level analysis.

Consequently, future work must prioritize aggregating a large-scale dataset of empirical BDC topologies to train automated feature-extraction pipelines. By leveraging computer vision techniques, such as semantic segmentation, we aim to automate the translation of raw aerial imagery into discrete spatial graphs. Resolving this challenge is a critical prerequisite for deploying formal topological metrics and enabling direct, quantitative comparisons between our ABM’s simulated traces and the physical ground truth.

Figure 15. The feature extraction bottleneck in empirical validation. While standard preprocessing (b) is insufficient for isolating high-frequency spatial networks from raw imagery (a), manual annotation (c) successfully yields the topological graphs (d) required for analysis. Automating the transition from (a) to (d) via computer vision represents a critical next step for large-scale environmental modeling.

Figure 15. The feature extraction bottleneck in empirical validation. While standard preprocessing (b) is insufficient for isolating high-frequency spatial networks from raw imagery (a), manual annotation (c) successfully yields the topological graphs (d) required for analysis. Automating the transition from (a) to (d) via computer vision represents a critical next step for large-scale environmental modeling.

4.2. Quantitative Comparison via Spatial Network Analysis

An initial qualitative comparison between empirical and simulated spatial traces is shown in Figure 16. The empirical data reflect three months of preprocessed aerial activity, loaded directly into the simulation suite, contrasted against the four-month temporal landmark of a four-agent colony (Figure 13(a)). Despite exhibiting similar macroscopic structures, such as primary foraging corridors and radial expansion, visual correlation remains insufficient for rigorous scientific validation.

To mathematically measure structural similarity and operational throughput, we propose integrating geographic information science (GISc) with network theory. The first methodological step requires extracting discrete spatial graphs $G=\{N,L\}$, where nodes N represent georeferenced junctions and links L represent connecting pathways, from the continuous environmental matrices of both the MAS and empirical pipelines.

Following graph extraction, the second step involves directly comparing these topologies using established network theory measures:

• Degree Distribution ($P\left(k\right)$) and Clustering Coefficient ($C\left(i\right)$): To evaluate the local connectivity and structural clustering of trail intersections.
• Average Shortest Path Length ($\langle l\rangle$): To quantify the global logistical efficiency of biomass transport to the lodge.
• Betweenness Centrality ($g\left(i\right)$): To identify bottleneck corridors and key spatial nodes.

By calculating these properties, future research can mathematically validate whether the simulated Explorer and Builder heuristics successfully replicate the structural complexity and spatial efficiency of natural ecosystem engineering.

Figure 16. Abstraction of emergent environmental traces into complex spatial networks. By converting both physical ground truth (a) and simulated outputs (b) into discrete graph structures, the morphological features of the landscapes can be quantitatively compared using topological metrics such as degree distribution and betweenness centrality.

Figure 16. Abstraction of emergent environmental traces into complex spatial networks. By converting both physical ground truth (a) and simulated outputs (b) into discrete graph structures, the morphological features of the landscapes can be quantitatively compared using topological metrics such as degree distribution and betweenness centrality.

4.3. Open-Source Framework and Community Engagement

Ultimately, the ambition of this work extends beyond simulating a single species. By providing this MAS framework as an open-source, process-based computational platform [37], we aim to establish a robust testbench for the wider interdisciplinary community. We actively invite ecologists, landscape architects, and control engineers to engage with, challenge, and expand upon this framework. Through collaborative iteration, we hope this framework will accelerate research into non-human environmental agency, serving as a catalyst for the synthesis of truly regenerative, adaptive landscapes.