Self-organized biotectonics of termite nests

To answer these questions, we first collected, scanned, and statistically analyzed the structure of nests in the species Apicotermes lamani. Two nests (MeMo14, Fig. 1A, and MeMo13, Fig. 1B) were collected near Libreville, Gabon (SI Appendix, Fig. S1), although one of these (MeMo13) showed too much internal damage for a statistical analysis. A third nest (MeMo80, Fig. 1C) was collected in 2008 near Pointe Noire in the Republic of the Congo (SI Appendix, Fig. S1); in this nest, we measured that workers had length 4 to 5 mm and height 2 to 3 mm, while soldiers had length 8 to 10 mm with a height of 2 to 3 mm (12) (SI Appendix, Table S1). We digitized and imaged the excavated nests with a medical X-ray computed-tomography (CT) scanner (Materials and Methods), allowing for nondestructive exploration of the complete three-dimensional internal structure of each nest (Movie S1).

Both nests showed a similar external appearance and, based on our X-ray CT scans, also shared a similar internal architecture consisting of many floor layers arranged in parallel, with several vertical pillars connecting adjacent floors (Fig. 1 A and B). To investigate these structures, we systematically searched our two digitized examples for open corridors connecting floors (Materials and Methods), and we identified several occurrences of two types of ramps (Fig. 1C): linear ramps (15 in MeMo80 and 23 in MeMo14), in which a path connects floors along a linear incline, and helicoidal ramps (5 in MeMo80), in which the path spirals around a central pillar as it spans floors. The helicoidal ramps had no clearly preferred chirality (three left, two right), although a pair of co-occurring left-handed and right-handed helicoidal ramps was identified (Fig. 1 C, ii).

We then examined the internal architecture of the two A. lamani nests (MeMo80 and MeMo14; SI Appendix, Table S1) for which we had data, to quantify the spacing and arrangement of floors and ramps. For each nest, we assembled a series of 359 and 177 vertical slices, respectively (Fig. 2A), and measured the thickness of each floor and the vertical spacing between floors (Fig. 2B). In total, our dataset included 35 floors (17 in MeMo80 and 18 in MeMo14), sufficient to conduct separate analyses on each of the two nests. The floor thickness was on average 1.7 mm (median = 1.6 mm) in nest MeMo80 and on average 1.8 mm (median = 1.2 mm) in nest MeMo14. The floor spacing was on average 4.6 mm (median = 4.7 mm) in nest MeMo80 and on average 7.2 mm (median = 7.7 mm) in nest MeMo14. For both nests, the variance in floor thickness was low (coefficient of variation [CV] = 0.31 in MeMo80, CV = 0.59 in MeMo14), as was the variance in floor spacing (CV = 0.19 in MeMo80, CV = 0.22 in MeMo14), demonstrating that the flooring of the termite nests is arranged at regular vertical frequencies, with a strong continuity between adjacent floors (SI Appendix, Fig. S1 B and C).

Given this consistent arrangement of floors, we next investigated the positioning of ramps within the nests by measuring the horizontal distance from each ramp to the nearest ramp on the same floor (Fig. 2C). In MeMo80, we found a mean same-floor distance of 27.4 mm (median = 25.4 mm, CV = 0.41), and in MeMo14 this mean distance was 25.7 mm (median = 27.1 mm, CV = 0.33); these both deviated from a null model of randomly scattered ramps (Kolmogorov–Smirnov [KS] test, $P=0.047$ for MeMo80, $P=0.022$ for MeMo14; SI Appendix). We compared these measurements with the horizontal distance between ramps on adjacent floors. Because these ramps were often directly connected, this distance was much smaller for both nests, with a mean of 12.9 mm (median = 9.2 mm, CV = 0.79) in MeMo80 and a mean of 7.1 mm (median = 6.5 mm, CV = 0.47) in MeMo14, deviating in the opposite direction from the null model (KS test, $P=0.046$ for MeMo80, $P=0.040$ for MeMo14). Hence, while ramps tended to be spaced at fairly large, regular intervals across a floor, they tended to be positioned near ramps on adjacent floors, allowing for efficient vertical traversal between nest levels.

To address how these structures are formed, we turn to a series of observations, old and new, to motivate a minimal continuum model of nest construction. With no evidence of design or a designer, termites build in response to local cues such as the nearby mound structure and secreted pheromones (17, 20). The resulting architecture both enables and constrains the movement of pheromones and termites and thereby modifies behavior, so that nest construction can be seen as a result of a feedback loop linking physical and behavioral dynamics (Fig. 3A). The architecture of the nest dictates which spaces are accessible to termite workers, the density of termite workers in turn controls the concentration of secreted pheromone, and the information carried by the pheromone profile in turn serves as a template for the ongoing remodeling of the nest architecture, thereby completing the feedback loop of nest morphogenesis.

We therefore model the spatiotemporal dynamics of three fields that depend on location denoted by the vector $\mathbf{\text{x}}$ and time $t$: the nest material density $u\left(\mathbf{\text{x}},t\right)$, scaled from no nest material ($u=0$) to fully compacted ($u=1$), as well as the termite worker density $n\left(\mathbf{\text{x}},t\right)$ and the pheromone concentration $\rho \left(\mathbf{\text{x}},t\right)$, both in units of mass per volume. These three fields jointly evolve according to a set of conservation laws, expressed as the following differential equations (SI Appendix) for the nest material, termite workers, and pheromone levels:where ${\partial }_t$ and ${\partial }_z$ are the differential operators in time and height, respectively; $\nabla \cdot$ is the three-dimensional divergence operator; and $\nabla$ is the three-dimensional gradient operator. Here Eq. 1 reflects the addition and removal of the nest material, with $f_{+}$ and $f_{-}$ denoting the building and removal rates of dirt (whose functional form is discussed below). Moreover, the poroelastic diffusivity of dirt $g$ captures its capacity to settle under the influence of gravity (21); this term breaks the rotational symmetry of the model equations and establishes a defined vertical orientation for the nest. Eq. 2 reflects the dynamics of termites carrying dirt $n_{+}$ or not carrying dirt $n_{-}$, with a flux that has two components: an effective diffusivity that is proportional to the amount of open space, i.e., $\left(1-u\right)D$, as well as a chemotactic term that drives termite workers into open and low-density nest regions with chemotactic coefficient $\chi$ (22, 23). The final term represents switching between $n_{+}$ and $n_{-}$, with the pellet size $k$ denoting the average amount of dirt transported by a worker. Finally, Eq. 3 reflects the dynamics of pheromones as they diffuse with a diffusion coefficient $\delta$, are secreted into deposited soil at a level $H$ per pellet (19), and degrade through evaporation at a rate $\gamma$. We note that the pheromone is relatively nonvolatile and thus not advected by the ambient fluid, since the subterranean nests of Apicotermes do not typically have noticeable temperature gradients (24).

To complete the formulation of our model, we require functional forms for the building rate $f_{+}$ and removal rate $f_{-}$. Although these rates have not been directly measured in A. lamani termites, the building rate has been investigated for the neotropical termite Procornitermes araujoi (25), and both rates have also been quantified in the ant L. niger (19). In all cases, we find that the local nest density scales directly with the building rate and inversely with the removal rate, resulting in high-density regions becoming increasingly compacted while low-density regions are excavated. Moreover, this positive feedback is mediated by a pheromone added by the workers to the building material (Fig. 3B), so that the termites are more likely to deposit dirt in regions marked by pheromone as being sites of active remodeling and so that building does not take place either where there is no nest material ($u=0$) or where the nest material is already fully compacted ($u=1$). It is reasonable also that the removal rate scales with the density of termite workers not carrying dirt and with the density of dirt $u$ available for removal. Given these basic considerations, the simplest possible functional forms for $f_{+}$ and $f_{-}$ arewhere $r_{\pm }$ are rate constants. With these choices, the dynamics of nest morphogenesis as described by our model can be captured by only a small number of nondimensional parameters: the scaled pheromone potency $H{r}_{+}/r_{-}$; the scaled evaporative flux $\gamma /r_{-}$; the pellet size $k$; and the scaled diffusivities $\delta /D$, $\chi /D$, and $g/D$ (SI Appendix). For simplicity, and to directly study the self-contained organization of this system, we model the nest as a closed three-dimensional domain, using no-flux boundary conditions, such that when integrated over the nest domain the total number of termites $\left|n\right|$ and the total quantity of dirt $u$ are conserved. With this formulation, our model describes the simultaneous disassembly and reassembly of the nest, which are features of both initial nest construction and its ongoing renovation during its lifespan (19).

We conducted 500 independent simulations of nest construction by initializing our model each time with uniform nest density $u\left(\mathbf{\text{x}},t\right)$ and pheromone concentration $\rho \left(\mathbf{\text{x}},t\right)$, and with a randomly dispersed density of termite workers $n\left(\mathbf{\text{x}},t\right)$, using measured and inferred parameter values from several studies on termite nest structure and behavior (SI Appendix, Table S2) and a numerical implementation of a stable Euler differencing scheme (Materials and Methods) run until a steady-state distribution of floor thickness and separation is achieved. Nest construction time ranged from 25.5 to 94.6 d (95% confidence interval; SI Appendix, Fig. S6). The model generated simulated nests (SI Appendix, Fig. S3 A–C) with regularly spaced floors of roughly planar orientation (orthogonal to the direction of gravity in the model) and a floor thickness that scales as ${\left(\delta /\gamma \right)}^{1/2}$ in some regimes (Materials and Methods and SI Appendix, Fig. S3D). For a pheromone diffusivity $\delta =100\hspace{0.17em}\mathrm{\mu }{\mathrm{\text{m}}}^2$/s and an intermediate degradation rate $\gamma =1.0$/h, this gives an average prediction of a 1.74-mm floor thickness, in line with our observed measurements of 1.7 and 1.8 mm in MeMo80 and MeMo14, respectively. To evaluate model sensitivity, we also allow all parameters to vary across realistic ranges (SI Appendix, Table S2)—including three orders of magnitude for the degradation rate $\gamma$, as this value is the most variable in the literature (19). Across all parameter ranges, we observed floor thicknesses in the relatively constrained range 1.2 to 2.5 mm, with the greatest sensitivity to changes in the diffusivity $\delta$ and degradation rate $\gamma$ as predicted by our scaling analysis (SI Appendix, Fig. S4).

Our simulated nests contained surface edges (topological defects, Fig. 3C) spanning across floors and sometimes forming traversable connections (Fig. 3D) corresponding to linear ramps—in which one floor terminates to provide the vertical space necessary for a simple ramp between adjacent floors—and helicoidal ramps, which wind about a vertical pillar that serves as a dislocation line spanning multiple floors (26) (SI Appendix, Fig. S5). By identifying the eventual location of a ramp, we observe the dynamical process by which ramps emerge from these defects during the construction process (Fig. 4A). Ramps first arise from regularly spaced floors in adjacent nest regions that are misaligned at the boundary via edge dislocation (in which an intermediate floor terminates at the boundary), which are connected to form a linear ramp to the floors above or below, or via screw dislocation (rotational misalignment about a slip plane), providing the axis of rotation about which a helicoidal ramp can be built (Fig. 3C). This misalignment of floor spacing is required for ramp formation; in follow-up simulations that initialized the mound with perfect floor regularity, ramps were altogether absent (SI Appendix). Moreover, when varying our model parameters, we find that two parameters control whether helicoidal ramps form at steady state: the scaled pheromone potency $H{r}_{+}/r_{-}$ and the scaled pheromone evaporation rate $\gamma /r_{-}$, yielding a phase space for their emergence (Fig. 4B). This is consistent with previous agent-based simulations of ant nest construction that pointed to the central role of the pheromone evaporation rate in influencing nest structures (19).

A fly-through of our simulated nests (Movie S2) reveals the structural similarity of the simulated nests produced by our model with the scanned A. lamani nests. To quantitatively evaluate this similarity, we gathered an identical set of statistics to describe the simulated nest topologies. We found that floor thickness and spacing between floors were similarly consistent in both the simulated and scanned nests, indicated by distributions with a similarly narrow variance (Fig. 4C). Moreover, the horizontal spacing between neighboring ramps on the same floor was greater than the horizontal spacing between ramps on adjacent floors that often connected directly, as observed in the scanned nests (Fig. 4D). To assess the similarity of the simulated and scanned nests on large scales, we carried out a spectral analysis of the average nest density as a function of height for both the scanned and simulated nests and found that these density profiles are very similar (Fig. 4E). Together, these metrics indicate that the simulated nests resemble natural nests across global and local scales.

Our minimal theoretical framework links three spatiotemporal fields involved in the biotectonics of termite nests: dirt to constitute the nest, termite workers to shape the nest, and secreted pheromone to mark active regions of the nest. They allow us to capture two key geometrical and topological features of termite nests, namely the regular vertical spacing of floors and the horizontal spacing of ramps on the same or adjacent floors, as well as the months-long timing of the nest construction process (12). This feedback loop involving architecture, behavior, and information that drives nest construction is quite general in biological systems; relatedly, it is implicated also in the macromorphogenesis of termite mounds (18) and might provide insight into emergence of other animal architectures, such as ant nests or mounds and nests of other termite subfamilies (1, 19, 27). Understanding the complex ways in which these systems differ, while still adhering to the same fundamental principles, is a task for the future. Nevertheless, our finding that key termite nest structures emerge spontaneously in our model for a range of parameter values points to it being a nonequilibrium steady state that is robust and adaptable—a direct consequence of the ability of a colony to both lay down and respond to cues in its environment.

An important open question is the role of heat gradients and airflow within the nest (5, 6, 15), which are determined by the nest structure and in turn help to propagate pheromones and other signals, such as carbon dioxide levels (28). Although the diffusion rate of pheromones can range widely depending on their molecular weight and chemical properties, the pheromones we consider here have low diffusivities (20), and thus only provide a localized signal useful in nest construction. Given these low diffusivities, heat-driven airflow would result in a negligible change when incorporated into our model (SI Appendix, Fig. S6), although other potential odors with a higher diffusivity, such as those secreted by the queen, would be more susceptible to internal nest flow and hence more important to the larger-scale problem of the global templating and emergence of overall mound shape (18, 20). In the subterranean nests of Apicotermes, we expect that the role of heat gradients and airflow is minor relative to that in above-ground nests which are directly subjected to sunlight and air currents. Furthermore, because termite workers can move with relative ease through each nest, the initial dispersion pattern of workers within the nest would have minimal effect on the overall nest structure (SI Appendix, Fig. S7).

In constructing the functional, palatial structures of their nests, termites demonstrate a remarkable capacity for spontaneous self-organization through stigmergic interactions, coupling their behavior to deposited cues in their local environment—a process that stands in stark contrast to the planned building processes that give rise to human architecture. By studying the natural architectures that termites and other insects build and inhabit, and by understanding their ability to leverage their surroundings for information transfer, we can better understand how multiagent system approaches can be adapted for our own lives, from swarm intelligence to architectural design (29). After all, while we as humans have a mere $\sim$400,000 y of experience with constructing primitive structures (30), mound-building termites of the family Termitidae have spent close to 50 million y (31) developing one of the grandest examples of architecture in the natural world.

We thank A. Robert for providing nest MeMo80 and the natural history museum in Paris for providing nests MeMo13 and MeMo14, both collected by P. P. Grassé. We are indebted to our reviewers for their careful reading and insightful suggestions. This work was supported by NSF Grant DGE-1144152 (to A.H.), Agence Nationale de la Recherche Grant ANR-06-BYOS-0008 (to G.T.), a China Scholarship Council PhD Grant (to L.G.), and the NSF Physics of Living Systems Grant PHY1606895 (to L.M.).