A Mathematical Approach to Modelling Termite Mound Organization and Population Dynamics

Keywords: Termites, Mound morphogenesis, Tunneling, Food transportation, Population modelling, Stigmergy

Abstract

Termites are tiny but truly fascinating organisms that live based on social structures and models that are beyond our imagination. They are remarkably efficient and adaptable through subtle signals and complex interactions; they are not just small entities but work together as a delicate, fine-tuned machine. This paper explores the foraging behavior and efficiency in food transportation of termites in a simulated environment. Next, the distribution and estimation of termite population using spatial statistics-based and diffusion-based methods is introduced. The paper also examines how the architectural structure of termite mounds is shaped and regulated by pheromones, material density, and termite worker density using positive feedback-loop mechanisms. Moreover, it describes how internal mound structure growth can be modelled as a simple function of wall concavity. The last part of the paper discusses how the outer mound structure is constructed and maintained because of feedback-loop mechanisms and termite behaviour, which itself is influenced by various physical factors like temperature, heat and airflow. This paper provides an in-depth understanding of termite organization, population, and ecological design in natural and engineered environments.

Introduction

Studying termite populations, behaviors, and colony structures is essential for understanding their species and their role and impact in nature. However, their subterranean habitats, concealed nests, and intricate social organization present significant challenges to traditional observation and measurement methods. Confronted by these challenges, some researchers turned to mathematical and computational methods to study termites' behavior and colony dynamics. These methods use data from natural and experimental studies, exploring food transportation, foraging efficiency, population estimation, and mound construction in termites. Termite food transportation strategies vary in efficiency with methods such as bucket brigades and direct deposition in tunnels. Termites’ foraging behavior is studied through adaptive tunneling, optimizing routes to food sources. Population estimation methods, including spatial statistics and diffusion-based models, use metrics like tunnel complexity and termite dispersion patterns to provide non-invasive and accurate assessments of colony size. Models of internal mound geometry reveal how termites respond to pheromone gradients to construct ramps, floors, and tunnels, while accretive growth models explain the characteristic saddle-shaped inner nest structures observed in some species. Additionally, simulations of airflow and heat transfer show how physical and chemical dynamics influence overall mound shape and morphology.

Food Transportation in Termites

Insect societies are composed of individual organisms having their own tasks. However, their ecological success would not be possible if they only worked alone. They often work together, and resource transportation is one of their communal goals. Social insects have unique behaviors in collecting resources and transporting them. The transport method can be classified as follows: individual transportation, bucket brigades, and multistage partitioned tasks with indirect transfer. Individual transportation is a form of transportation where an individual worker directly deposits a piece of resource from the site of the resource to the nest. Bucket brigades, on the other hand, involve an individual transferring a piece of resource to another individual when they come across each other on their way back to the nest. Lastly, in transport using multistage partitioned tasks with indirect transfer, individual workers carry food and drop it in the middle of the pathway and another worker picks it up (Fig. 1) (Anderson et al., 2002).

Fig. 1. Three different transport method in termites. a) individual transportation, b) bucket brigades, and c) multistage partitioned tasks with indirect transfer (Anderson et al., 2002).

Food Searching Shortcuts

Subterranean termites search for food resources through underground tunnels, which can be as far as 100m away from the nest. Other social insects transport their food back to the nest using a more optimized path on the ground. This is more difficult in subterranean termites since they need to excavate a new path, but this shortcut may benefit the colony in the long term (Michael et al., 2023). Michael et al. studied the tunneling behavior of termites when three different paths were given: straight (S), detour (D), and detour + twisting (DT), as Fig. 2 shows. A colony of Coptotermes formosanus Shiraki was used in this study. Arenas (18 x 18 x 0.2) with three different pre-formed paths were used and pictures of the arenas were taken every 6 hours.

Fig. 2. Experimental setup for the behavioral assay where a moistened cellulose media pad (2 × 2 cm) was placed on the top-left side of each arena as a food source (left). An example of tunneling activities of C. formosanus in each arena according to the time (0 h, 12 h and 24 h) (Michael et al., 2023).

Termites in the S arena continuously used the existing tunnels and tunnel distance did not change over time (ANOVA, F = 0.037, df = 4, 20, P = 0.997). The P value greater than 0.05 shows that the change in distance is not statistically significant in this case. However, termites in two other arenas actively excavated new branching tunnels. The tunnel distances from the entrance in the arena to food in both arenas were significantly reduced over time (ANOVA, D: F = 5.30, df = 4, 20, P < 0.01; DT: F = 35.56, df = 4, 20, P < 0.01) (Michael et al., 2023). In both cases, the P values smaller than 0.05 demonstrate that the change in distance is significant (Fig. 3). Interestingly, in the D arena, travel distance began to decrease significantly at 18 hours, while in the DT arena, it started in earlier time, from 6 hours, and continued to decrease constantly. This study suggests that termites tend to find shortcuts to increase food transport efficiency.

Fig. 3. The shortest distance of tunnels to food from the entrance in the arenas used by C. formosanus over time. Red, green, and blue indicate different arena types, S, D, and DT. Black dots inside the box represent the average and the black line indicates connection of the mean line over time (Michael et al., 2023).

Food Transportation Efficiency in Termites

In studies about food transportation efficiency in termites by Lee et al., the researchers used simulated termite to study food transportation behaviour of termites. The tunnel has a narrow width, about 4-5 mm, that only 4 individuals can pass through at the same time. In a space composed of 1,000 grid cells, they built the tunnel using the function below (Eq. 1). Other factors that are likely to delay transposition time are not only the narrow tunnels, but also irregular walls due to heterogeneous soil particles (Fig 4.C), and pore formed between soil particles (Fig 4.D) (Lee et al., 2022).

Fig. 4. (A) Illustration of a simulated termite traveling through a sinusoidal tunnel consisting of a food site and a nest at both ends, (B) a termite passing through a sinusoidal tunnel, (C) a hesitating termite in a simplified irregularity, and (D) a termite delaying time within a simplified soil pore (Lee et al., 2022).

In the transport process, food transfer only occurs between individuals facing each other. When a simulated termite carrying food encounters another termite without food, the food transfer is triggered according to the probability of food transfer, P_T (= 0.0, 0.1, ..., 1.0). There might be a food loss, also according to the probability, P_L (= 0.0, 0.03, ..., 0.3). If an individual meet two or more approaching individuals, one of the individuals is randomly selected, and food transfer is performed, and individuals reverse the direction of their movement (Fig. 5.A). If an individual carries food, transfer does not happen, and they pass each other without interaction (Fig. 5.B). When a food particle is partially lost n times when arrives at the nest, the amount of food is calculated as 1- (n × P_L).

Fig. 5. (A) A simulated termite carrying food encounters another termite without food. The transfer occurs through trophallaxis, a mouth-to-mouth exchange. (B) Two simulated termites without food meet in a tunnel, passing each other without stopping or altering their paths (Lee & Park, 2024).

Some termites would occasionally stop in their path. It causes either other termites to stop as well or, they turn vertically to pass through the gap between the individuals and the tunnel. This is referred as a “traffic jam”. In the simulation, this traffic jam is defined as when four simulated termites meet at the same site and the sum of their direction vector becomes zero. Once a traffic jam occurred, four individuals stayed at that site for 100-time steps.

The food transportation efficiency, E, is calculated according to Eq. 3, where “food(t)” represents the number of food particles delivered to the nest at time t, and N₀ is the number of termites participating in transport. s is a scaling factor with a value of 10,000, which they used to prevent E from becoming too small. P_C, P_T, and P_L value was calculated in a simpler way, using the relation P_C = (k₁-1) × 0.05, P_T = (k₂-1) × 0.05, and P_L = (k3-1) × 0.03, where k₁, k₂, and k₃ (=1, 2, ..., 11) are the variable values (Lee et al., 2022). The variables are presented in Table 1 with their definition for better understating.

For the simulation, 14,641 (=11 × 11 × 11 × 11) combinations of the four variables were performed, each replicated 15 times. The study tested two scenarios: the bucket brigade, where P_T > 0 and the direct deposition, where P_T = 0. When there was no food loss, E increased linearly regardless of P_T in the region t ≤ 1500. However, when the food transport process was longer (t > 1500, P_L = 0.15), E decreased significantly with increasing P_T. Researchers also found various effects of variables on E. First, as P_L increased, E decreased regardless of P_C and P_T. When P_L = 0 and N₀ is low, there was no clear effect on the E pattern. However, when P_C increased, E decreased because termites often stopped in high-curvature zones.  Furthermore, when P_L was higher than 0, E decreased rapidly as P_T increased, and smaller decreases in E when P_C increased. This shows that P_T plays an important role when P_L is greater than 0. As P_L increased more, P_T became more significant as well. Therefore, it can be concluded that the preferable method of transportation depends on the tunnel conditions. In narrow and short-distance tunnels, direct disposition is preferable due to frequent traffic jams. On the other hand, in a wider and longer tunnel, the bucket brigade is better since only a few traffic jams occur. It also improves the stability of food transport, even though there is a slight loss in efficiency (Lee et al., 2022).

Table 1. Definition of some variables used in the food transport efficiency simulation (Lee et al., 2022).

Fig. 6. Food transport efficiency, E (s = 1, T = 1), for food transfer probability, P_T (= 0.1, 0.2, …, 1.0) and food loss values, P_L (= 0.0 and 0.15), occurring during the transfer, where the tunnel curvature effect, P_C, and the number of simulated termites, N₀, was 0.25 and 50, respectively (Lee et al., 2022).

Population Modelling

Estimating populations of termites is challenging because they are hidden, highly social insects with complex nesting structures. Traditional methods like direct sampling and mark-release-recapture (MRR) involve physically capturing termites, marking them, and using the ratio of marked to unmarked individuals to estimate population sizes. The multiple-capture-recapture (MCR) method extends this approach by capturing and releasing termites multiple times for more accurate estimations. While non-destructive, MRR and MCR rely on assumptions of uniform termite distribution and equal recapture probability, which often do not hold due to termites' clustering behavior and complex social organization limiting method effectiveness. In response, newer methods, including diffusion-based models and spatial statistics, offer indirect, minimally invasive approaches for more accurate population estimates without disturbing colonies.

Spatial Statistics Model

In 2022, Sang-Hee Lee and Seung Woo Sim (Lee et al., 2008) developed a spatial statistics-based model to estimate termite populations by analyzing tunnel networks using metrics such as Fractal Dimension (FD), Local Density (LD), and Join Count Statistic (JCS), each of which reflects population size. To find these metrics, high-resolution images of the tunnel network are processed digitally to highlight the tunnel structures. This processed image is then overlaid with a grid of square cells, each representing a specific area. If a cell contains any part of the tunnel, it is marked as filled, if not, it is marked as unfilled. The cell size is gradually reduced, and the marking process is repeated for each trial. This cell-based grid representation of the tunnel network allows for standardized calculations like FD, LD, and JCS, providing a detailed way to estimate termite population size based on the complexity and density of their tunnel systems.

Fractal Dimension (FD) is found using the box-counting method, where grids of progressively smaller sizes cover an image of the tunnel pattern. Researchers count the number of grid boxes that intersect the tunnels at each grid size and plot this count against grid size on a log-log scale. The slope of the resulting line provides the FD, indicating network complexity (Fig. 7). Local Density (LD) measures tunnel concentration within specific areas of the network by counting the number of tunnel cells within a defined space and averaging this count across multiple areas. High local density suggests that tunnels are densely clustered, indicating more termite activity in those areas. Join Count Statistic (JCS) assesses the number of intersections or "joins" between tunnel cells, measuring the connectivity within the network. A higher JCS value implies that the tunnel system is highly interconnected, a feature often associated with larger populations. These metrics collectively depict tunnel structure, linking network characteristics to termite population estimates.

Fig. 7. Calculating Fractal Dimension (FD) of a termite tunnel pattern using the box-counting method. s (size of box) decreases and n (filled boxes increases). For each s and n a log-log plot is created to find the slope which represents FD, complexity (Sim et al., 2022).

To predict the population of a target colony, researchers use the k-nearest neighbors (k-NN) algorithm, which predicts the population by comparing the tunnel metrics of an unknown colony to those of colonies with known populations. In this method, each colony with a known population and set of FD, LD, and JCS values serves as a 'neighbor.' For a new colony, the k-NN model identifies the number of neighbors (k) most like the target colony (based on the metrics) and estimates population size by averaging the populations of these neighbors.

Diffusion Based Model

In 2008, Sang-Hee Lee and Nan-Yao Su (Lee et al., 2008) developed a model used to estimate the population of a termite colony based on the diffusion of the termites in a controlled environment to correlate probability density to population. 11 arenas, each connected by 5 meters of tubing, were set up in a row to simulate a termite colony's foraging environment (Fig. 8). This design mimicked the spatial structure of termite activity, allowing the researchers to study how termites naturally disperse across interconnected spaces.

Fig. 8. Set up of the diffusion-based model. Pictures were taken once a day and in the meantime the arena was kept in dark. (Lee et al., 2008).

10000 worker termites and 1500 soldier termites were roughly divided into 11 arenas. They were allowed to disperse for a week when the researcher added 240 marked worker termites into the central arena. Over the following days, researchers took daily photographs of each arena to count the number of marked and unmarked termites present in each location. This served as observational data on termite dispersion patterns, which was used to correlate diffusion behavior with population estimates. This model accounts for termites’ natural dispersion patterns and non-uniform distribution, avoiding the reliance on assumptions of uniform mixing and equal recapture probabilities. To find the mean number of termites in each area the researchers used Equation 4.

Here, N is the mean number of termites in the given area, α is the total number of arenas or observations included in the calculation, i is the arena index, and M is the total number of marked termites released. Then, n_i is the count of unmarked termites in the i-th arena, and m_i is the count of marked termites in the i-th arena. ω is a diffusion constant that represents the rate of termite dispersion, where the fitting function P(x)~exp(-ωx²) is used to find ω. x_i is the distance of the i-th arena from the central release point, and L is the length of the tubing connecting the arenas, which ensures scaling for termite movement. √(ω/π)*exp(-ωx_i²) represents the probability of termite presence at a given distance x_i from the release point. The variation of N is calculated with the same variables and is given by Equation 5:

Based on the calculations above, the estimated population was compared to the known population using a Z-test. When Z > 1.96, the N was significantly different from the known population.

Modelling Internal Mound Geometry

Termite mounds remain one of nature’s most impressive examples of collective behaviour in insects— a result of highly sophisticated colony organization. Mounds provide termites with an essential microniche; their internal geometry allows for autonomous regulatory processes, such as ventilation and heat transfer, and bulk termite movement (Heyde et al., 2021; Calovi et al., 2019). This striking complexity could hardly be coordinated by one termite alone, yet mound construction still shows uniformity and coherence at the global scale. Research suggests that mechanisms exist to drive local building activity, indicating that termites rely on their local environmental conditions to produce global order (Heyde et al., 2021). Proposed organizational factors include pheromonal gradients, produced either by workers or queens, wall curvature and worker aggregations (Heyde et al., 2021; Facchini et al., 2020; Nakanishi et al., 2017). This section aims to provide several comprehensive models to describe the construction and development of internal mound structures, such as ramps, floors, tunnels and chambers.

Pheromone-induced Construction of Internal Mound Features

Heyde et al. model the emergence of internal features in mounds of the termite genus Apicotermes, which construct small, oval-shaped nests located primarily underground. They observed distinct structural motifs, including both linear and helical ramps connecting parallel floors arranged at regular vertical levels. Continuity between floor thickness and spacing was shown between different nests of Apicotermes lamani, indicating a universal building regime (Fig. 9).

Fig. 9. Digital reconstruction of A. lamani nests via computed tomography (A and B) Cross sections of two different nests show similar internal geometry with evenly spaced floors of similar thickness. (C) Cross section showing presence of (i) linear ramps and (ii) helical ramps (Heyde et al., 2021).

Heyde et al., propose a feedback loop between mound material density, termite worker density and secreted pheromones, which govern building behaviour. These are described by the following spatiotemporal fields, which depend on position vector x and time t: nest material density u(x,t) scaled from no material (u=0) to fully compact (u=1), worker termite density n(x,t) and pheromonal concentration p(x,t). Then, following a reaction-advection-diffusion framework, the three fields can be described by the following differential equations.

Equation 6 models the accumulation and erosion of nest material: ∂_tu , which is the change in nest material density with respect to time. The first term, ∂_zu * [g∂_zu] describes the vertical distribution of material over time, where ∂_zu is the change in material density with respect to the vertical height z , and g is the poroelastic diffusivity of dirt, which describes the capacity of dirt to settle under the influence of gravity. The f₊ and f_- terms describe the rate of nest material addition and removal, respectively, either by worker termites or external environmental influence.

Equation 7 models the spatial distribution of dirt-carrying (n₊) and non-dirt-carrying (n_{_}) termites over time under the influence of diffusion and pheromonal chemotaxis. The first term represents the flux of termites and is dependent on a diffusivity term (1 - u)(D∇n_±) and a chemotactic term Xn_±∇u. The diffusivity term comprises the gradient of termite density ∇n_± and is proportional to D, the effective diffusive constant, and the amount of open nest space (1 - u), which indicates that when nest material density u is high, termite movement decreases. The chemotactic term drives termites toward open, low-density nest regions, with a chemotactic coefficient X representing the tendency to move in response to a chemical stimulus. Finally, the ±(f - f₊) term represents the switch between n₊ and n_- (i.e. termites depositing and picking up dirt). This behaviour is related to the dirt addition and removal rates as well as the average dirt pellet size k, which affects how frequently termites will switch between dirt-carrying and non-dirt-carrying.

Finally, Equation 8 models pheromonal spread as a function of time, i.e. the rate of change of pheromone concentration ∂_tp. Cement pheromones are deposited in dirt material and will gradually diffuse and decay throughout the mound space. Here, the first term ∇*[δ∇p] describes pheromonal diffusion, dependent on a diffusion coefficient δ and the gradient of pheromone concentration ∇p. The Hf₊ term represents the pheromone secreted into deposited soil at a concentration H per dirt pellet. The final -γp term is the natural decay of pheromone over time with an evaporation rate constant γ. The relationship between all three equations is show in Figure 10.

Fig. 10. (A) Schematic of proposed model feedback loop between nest material density μ, termite worker density n, and secreted building pheromones p. (B) Local region of a nest showing the processes described by the model. Termite workers migrate preferentially to low-density regions and cannot travel through very high-density regions. Workers remove dirt arbitrarily in the nest but are most likely to deposit dirt near pheromones, which they release during deposition [Adapted from Heyde et al., 2021]

This model provides the basis for accurate nest simulation, which is numerically solved via a finite difference solver. By initializing the simulation with a uniform nest density and pheromone concentration, and a randomized dispersion of termite workers, nests are generated with regularly spaced floors of thickness (δ / γ)^0.5. Taking literature consistent values of δ = 100 μm²s^-1  and γ = 1.0 hr^-1, this gives an average thickness of 1.74 mm, comparable with observed mound measurements of 1.7 mm and 1.8 mm. The comparison of simulated and measured nest material power spectrums, which provide flood spacing and thickness information, are shown in Figure 11D. Remarkably, the model periodically generates both linear and helical ramps, emerging from defects in the construction process. When the regularly spaced floors become misaligned, edge and screw dislocations lead to the creation of linear and helical ramps, respectively (Fig. 11A and B). These construction misalignments are the direct result of pheromonal influence; Heyde et al. identify two parameters that control ramp development: the scaled pheromone potency Hr₊ / r_{_}, where r₊ and r_- are rate constants relating to the building and removal rates f₊ and f_-, and the scaled pheromone evaporation rate γ / r_- (Fig. 11C). These findings emphasize the role of pheromonal influence in effective mound construction.

Fig. 11. (A) Edge and screw dislocations in floor patterning, resulting from floor misalignment. Edge dislocations can give rise to linear ramps, while screw dislocations can lead to helicoidal ramps that pivot about a slip plane. (B) 3D nest simulation containing one linear ramp (i, blue) and two helicoidal ramps (ii and iii, red). (C) Heat plot for the frequency of helicoidal ramps in a simulated nest as a function of the two model parameters governing pheromone dynamics (potency and evaporation rate). (D) Power spectrum of nest density peaks sharply at a one-floor period for both simulated (gray) and natural (purple) nests, indicating regularly spaced floor structures and thickness [Adapted from Heyde et al., 2021].

Accretive Mound Growth Model based on Wall Curvature

Queen and cement pheromones are not the only proposed driving force of internal mound geometry. Among others are worker termite aggregations, which seem to attract other workers to building sites, pheromonal and gas gradients, and wall concavity (Calovi et al., 2019) (Facchini et al., 2020). Primarily, the extent of the curvature of mound walls seems to act as a cue for excavation and deposition, indicating that termite construction is strongly associated with surface concavity. In arboreal termites of the genus Nasutitermes, construction via material deposition occurs only on the edges of existing wall concavity, leading to tunnel bending and branching throughout the mound and resulting in saddle-like internal structures (Fig. 12) (Facchini et al., 2020). Facchini et al. propose an accretive growth model, where new material is only added at the surface of growing structures, much like crystals or corals. They summarize the results of complex construction dynamics (pheromones, aggregations, gas gradients, etc.) by looking at mound geometries from an isolated perspective and modelling structure growth as a function of its shape and curvature.

Fig. 12. Fragments of nests built by arboreal termite species Nasutitermes walkeri (a) and Nasutitermes ephratae (b). Notice the saddle-shaped curvature throughout the mound (Facchini et al., 2020).

The proposed model is based on the concept of phase fields and treats the two-phase material/ air mound composition as one continuous phase, a scalar field f across the whole mound geometry. The material/air interface can be retrieved as an iso-contour of f: a curve or line in the spatial domain where the field has the same value. Then, surface geometry is retrieved via differentiation. f is defined on the interval [0, 1] and the iso-surface value is given by f = f₀ . Conventionally, f(x) < f₀ represents a point inside the wall material, while represents a point in empty space. The scalar field f is given by Equation 9:

This represents the change in the phase field function over time, or the mound growth. The term -f^α(1 - f)^β is introduced to provide non-linearity to the model; it vanishes at f = 0, 1 and has a maximum at f = f₀ = α / (α + β) , damping any variation of f far from the iso-surface, which is the only space for additions of mound material by worker termites. In other words, there cannot be any growth in the empty mound space far from pre-existing structures nor any growth deep inside a pre-existing wall. The term d(∇*n) represents the divergence of the normal vector n, which describes the mean curvature of the iso-surface multiplied by a growth parameter d. This term mimics the construction behaviour of worker termites with respect to local curvature. When wall curvature is positive, f increases (deposition of material); when wall curvature is negative, f decreases (excavation of material). Finally, the last term represents the diffusion of mean curvature, compensating the previous growth term when local curvature is too high, and introducing a smoothing effect on the growing structures, as observed in natural mounds.

Facchini et al. further simplify this model by setting α = β =1, which assumes structural symmetry and gives f = f₀ = 0.5. They also assume ∇*n ∝ Δf, i.e. that the mean curvature is proportionate to the spatial diffusion of f. Equation 10 gives a finalized function for mound growth:

Note that the Δ²f term represents the higher order diffusion of f, which is involved in interfacial smoothing, as before. The model is then simulated on a cubic grid over time, showing frequent branching and merging of internal structures (Fig. 13). When compared to natural nests, a large part of the structural complexity is maintained (Fig. 14), indicating that complicated mound-building dynamics can be modelled as a response to curvature alone.

Fig. 13. Early evolution of the f = 0.5 iso-surface over time. Red circles indicate early branching and merging events, simulating tunnel formation (Facchini et al., 2020).

Fig. 14. Effect of f = 0 boundary on nest growth. (a) and (b) is a 2D vertical cross-section and a 3D rendering of the simulation at t=128. As the nest simulation approaches the f = 0 boundary, it closes in on itself to form a structure that resembles the nests of Nasutitermes (Facchini et al., 2020).

Mathematical Modelling of Outer Mound Morphology

Termite nests are a remarkable example of functional self-organization, demonstrating how structure and function emerge on multiple length and time scales in ecophysiology (Heyde et al., 2021). In 2019, Ocko et al. presented a model for termite mound morphogenesis that results from the collective behavior of termites in response to local information, producing distinctive morphologies shaped by the interaction of physical and behavioral factors. Their model explores the potential range of mature mounds using a limited numbers of parameters to determine relationships for mound shape and size.

Feedback Loop and Stigmergy

The shape of termite mounds results from the collective behavior of termites as they simultaneously respond to changes in the local environment. This process is characterized by a feedback loop involving both physical conditions and termite behaviors (Fig. 15). Initially, as the mound structure begins to form, it influences heat and air flow, which are also affected by changes in external temperature. Inside the mound, air current transport odors (pheromones or metabolic gases) that serve as signals, guiding termites’ construction behavior. Termite workers adjust their construction behavior in response to local odor concentration and then modify the mound structure, which completes the feedback loop. This system of actions and responses drives the mound’s development until a mature size and shape are achieved (Ocko et al., 2019).

Fig. 15. Schematic of model procedures over a single iteration. (A) The previous mound state is composed of disjoint mound regions. (B) As the ambient temperature oscillates, the heat profile and airflow field are calculated in the mound interior. (C) The odor profile is computed, and the mound wall shape is updated accordingly. Green or red cells indicate sites with an average odor level sufficiently above or below (respectively) the odor threshold, resulting in termites extending the mound outward or allowing the wall to collapse inward. (D) The mound shape is updated and used for the next iteration. This process is repeated until convergence to a steady-state morphology (Ocko et al., 2019).

Moreover, stigmergy is a paradoxical phenomenon of individual insects behaving in a decentralized way but building structures as if being centrally organized. Termites feedback loop is a representation of the positive and negative feedback processes that lead to a decentralized optimal construction of a functioning mound. Positive feedback encourages construction activity in an area where a high concentration of pheromones or odors attracts termite workers. Negative feedback regulates and prevents overbuilding, helping maintain balance. These processes guide termites’ behavior and lead to efficient mound construction (Oberst et al., 2020).

Model of Termite Mound

Ocko et al. model defines four subdomains, which are the mound interior, the ground, the surrounding air, and the wall, which is the boundary between the interior and the air at a fixed thickness. From these subdomains, they modelled the evolution of temperature, airflow, and odor fields as fluids in a porous medium.

The temperature in the mound T(r,t) varies due to thermal conduction and fluid convection. Heat diffuses from warmer to cooler areas within the mound’s interior and walls at a constant rate of speed, which is called thermal diffusivity D_T. This describes how quickly heat can move through the material of the mound due to temperature differences. Heat is also carried by an airflow field u(r,t) in the mound, which is called advection and depends on an advection weight ν. These processes follow the advection-diffusion equation (Eq. 11), which models how heat disperses due to both natural diffusion and airflow.

The airflow field u is modelled as being driven by heat gradients, following Darcy’s law (Eq. 12) as a conserved fluid. Darcy’s law is a principle that describes the flow of a fluid through a porous medium and it depends on the pressure difference (which is driven by heat gradients) and the properties of the medium.

where P(r,t) is fluid pressure, κ(r) is permeability, η is viscosity, and αρg is the buoyancy parameter, with α the thermal expansion coefficient.

The odor concentration field Ø(r,t) in the mound is modelled as the consequence of small particle diffusion and advection by airflows, with odor production. In other words, in addition to small particle diffusion, air currents, which are driven by heat gradients, carry odor particles through the mound. Odor is continuously produced at the nest from metabolic processes or pheromones released by termites, and its rate of production is represented by J. Equation 13 illustrates the odor dynamics for a central nest at ground height.

where ∈ is the time scale of the dynamics, δ(r) represents a point source of odor at the nest, and D(r) is the odor diffusivity. As larger colonies will tend to produce a greater quantity of odor, this high rate will result in a large mound radius.

Henceforth, the behavioral rules are incorporated into the model, along with physical factors and allow the model to update T(r,t), u(r,t), and Ø(r,t) over a diurnal period (daily cycle). The model takes into consideration four length scales, which are the wall thickness, the thermal diffusion length, the Peclet length of odor and a characteristic mound radius. Three dimensionless parameters, which define the morphogenetic processes defining steady-state mound shapes, are obtained by comparing the characteristic mound radius with other length scales. The Peclet number (Pe) is the ratio of the characteristic mound radius over the Peclet length of odor, the relative thickness (Th) is the ratio of wall thickness over the characteristic mound radius, and the Biot number (Bi) is the ratio of the thermal diffusion length over the characteristic mound radius (Ocko et al., 2019).

Simulations of Mound Morphogenesis

They conducted 500 independent simulations of mound morphogenesis based on numerical simulations using Equations 11, 12, 13. Each simulation started with the same conditions, which were a spherical mound wall of a specific, fixed thickness. The model could then evolve over time until it reaches a stable steady-shape morphology. At the same time, the changing size and shape of the mound were tracked (Fig. 16).

Fig. 16. Morphogenesis dynamics of simulated mounds. (A) Various steady-state mound morphologies are plotted in a 2D morphospace of mound sphericity and volume, relative to the critical mound radius. Points are colored and sized by the Bi and Th values used in the simulation. The two black trajectories show two possible paths through this morphospace, from an initially small and spherical mound to the mature shape. (B and C) Progression of mound geometries during morphogenesis. Each row corresponds to a single mound as it progresses along the corresponding trajectory in A (Ocko et al., 2019).

Furthermore, the simulations enabled to determine that mounds with low relative thickness and low Biot number tend to grow to larger volumes. They also found that mound sphericity was strongly linked to low Peclet number, while unusual mound shapes corresponded to high Peclet number. In that case, the movement of particles due to airflow is much stronger than conduction (the diffusion of particles randomly). It means that convection is the primary mechanism for transporting odor particles at high Peclet numbers, transporting the odor in the direction of growth and changing the mound shape (Fig. 17).

Fig. 17. Mapping model parameters to mature mound morphologies. (A) An array of nine simulated mound shapes for varying relative thickness Th and Peclet number Pe, for fixed Biot number. (B) Three simulated mound shapes for varying Bi for fixed Th and Pe. (C) Dimensionless parameter space of our morphogenesis model. Each dot corresponds to a single simulation, with dot position giving the values of Pe, Th, and Bi from the simulation, dot size proportional to the mound volume at steady state, and dot color corresponding to mound sphericity at steady state according to the color bar. (Ocko et al., 2019).

Conclusion

Like other social insects, termites use various food transportation methods, from individual transportation to more cooperative systems like bucket brigades. Studies on subterranean termites show they can find the shortest way to carry resources, which can be energy-intensive work in the short term but has the long-term benefits of shortened foraging paths. Simulations on food transport efficiency reveal how factors like tunnel curvature, food transfer probabilities, and transportation methods affect termites' efficiency in delivering food. These findings help us better understand their complex behavior as individuals and superorganisms.

The spatial statistics-based and diffusion-based methods offer innovative, non-invasive approaches to estimating termite populations, addressing limitations of traditional techniques like mark-release-recapture. Spatial statistics analyze tunnel networks using metrics like Fractal Dimension (FD), Local Density (LD), and Join Count Statistic (JCS), while diffusion-based methods rely on termite dispersion patterns and tools like the Z-test. Despite limitations such as computing costs and resolution needs for spatial models, and assumptions of homogeneity, age and motivation for foraging in diffusion experiments, these methods highlight the value of mathematical modeling for accurate, non-disruptive ecological studies.

The organization of internal mound structures relies on complicated building dynamics. Pheromonal influence on termite workers leads to the emergence of ramps and floor spacings and can be modelled using a reaction-advection-diffusion framework describing the change in density of pheromones, termites and mound material in a positive feedback loop. The results of pheromonal and other external driving forces can be uniquely summarized in an accretive growth model, which describes structure formation as a function of increasing and decreasing wall concavity.

The outer morphology of termite mounds is a result of a feedback loop involving physical conditions, like internal airflows, and behavioral dynamics. A mathematical modelling of mound has been done by Ocko et al. through the evolution of temperature, airflow, and odor fields. Their simulation showed that the Peclet number, the relative thickness and the Biot number all depends on the characteristic mound radius and have impacts on the morphology of the mound. This mathematical modelling is based on the concept of stigmergy, the collective behavior of independent termites results in highly organized and functional structures. These models have their own limitations but each of them gives us insights on the complex dynamics of termite mounds because studying them collectively provides a more comprehensive understanding.