A growth model driven by curvature reproduces geometric features of arboreal termite nests
Giulio Facchini, Alexandre Lazarescu, Andrea Perna, Stéphan Douady
AA growth model driven by curvature reproducesgeometric features of arboreal termite nests
G. Facchini , A. Lazarescu , A. Perna , and S. Douady Life Sciences Department, University of Roehampton, London, UK Institut de Recherche en Mathématique et Physique, UC Louvain,Louvain-la-Neuve, Belgium Laboratoire Matière et Systèmes Complexe, Université ParisDiderot, Paris, France
Abstract
We present a simple three-dimensional model to describe the autonomousexpansion of a substrate which grows driven by the local mean curvature ofits surface. The model aims to reproduce the nest construction process inarboreal
Nasutitermes termites, whose cooperation may similarly be mediatedby the shape of the structure they are walking on, for example focusing thebuilding activity of termites where local mean curvature is high. We adopt aphase-field model where the nest is described by one continuous scalar fieldand its growth is governed by a single nonlinear equation with one adjustableparameter d . When d is large enough the equation is linearly unstable andfairly reproduces a growth process where the initial walls expand, branch andmerge, while progressively invading all the available space, which is consistentwith the intricate structures of real nests. Interestingly, the linear problemassociated to our growth equation is analogous to the buckling of a thin elasticplate under symmetric in-plane compression which is also known to producerich pattern through non linear and secondary instabilities. We validated ourmodel by collecting nests of two species of arboreal Nasutitermes from thefield and imaging their structure with a micro-CT scanner. We found a strongresemblance between real and simulated nests, characterised by the emergenceof a characteristic length-scale and by the abundance of saddle-shaped surfaceswith zero-mean curvature which validates the choice of the driving mechanismof our growth model. a r X i v : . [ n li n . AO ] F e b Introduction
Self-organised growth phenomena are ubiquitous in physics, societies and bi-ology, with examples ranging from dunes and ripples in the sand, to crystalsand cities, to developing plants and embryos (see [4] for a review). All suchphenomena are associated with autonomous processes in which a substrateincreases its size while simultaneously developing a characteristic shape.In some systems growth is driven directly by volumetric expansion. Thisprocess is common for instance in living tissues, where cell proliferation canhappen both at the surface and inside the volume. In other systems, growthis driven mainly by accretion processes, whereby new material is added orremoved only at the surface of the growing object (crystals and concretions,but also corals and shells are example of this latter type of growth)[36]. Asa consequence, accretive growth is inevitably related to the dynamics of theinterfaces between the growing substrate and the external environment. Im-portantly the shape of the interface can contribute to focus the growth-drivingquantity at specific positions thus causing the onset of positive feedback loopsthat produce fingering and tip-splitting. The growth-driving quantity itselfcan be different in different systems, for example it could be pressure, as inthe seminal work of [34] - where a first fluid invades a more viscous one -, orit could be the concentration of a protein as in the more recent work of [9],aimed at modelling the branching morphogenesis of lungs. At small lengthscales negative feedback (viscous dissipation, capillarity) stabilizes the systemand fixes the typical scale of the final pattern [31]. Growth is then often aself-sustained process where the form of the substrate is the main ingredientas well as the main outcome of the growth process.A particular example of self-organised accretive growth is nest building bysocial insects. At first sight, collective nest building by social insects may lookdifferent from the examples above because the transport and deposition ofbuilding material is mediated by active agents - the insects - that are capableof multiple regulations of activity, and of complex movement patterns. Infact, there is extensive evidence that nest building can be influenced by anumber of factors, including insects responding to environmental gradients[18], crowding [41], the insects using their own body as a template [24], andindividual abilities such as proprioception [3] (see [32] for a review). Yet, whileinsects movement itself can follow complex rules, the decisions that insects takeof adding or removing pellets at a particular place obey simple “stigmergic”principles whereby the local configuration of the environment is the main drivethat determines the probability of deposition and collection of material [16,35]. When insect density is high the effect of such stigmergic regulationsis likely to dominate over the other factors and the nest growth resembles n accretive growth. In this continuum approach the construction processshould be described by a Turing-like system of coupled equations, one for eachof the interacting fields - for example density of agents, pheromones, CO concentration etc. - which are constrained by the shape of the structure underconstruction. In particular we expect the spatial boundaries to affect not onlythe spatial distribution but also the interactions among the different fields atplay. This suggests to consider an alternative and simplified approach whereonly such a spatial constraints are retained and the nest grows by itself as afunction of its own shape, similarly to a crystal or a cell’s tissue.Here, we focus in particular on modelling the nests built by arboreal ter-mites of the genus Nasutitermes . These are lightweight - but resistant - struc-tures that are usually built up on trees, typically around a branch or a treefork. Their structure is often isotropic and very homogeneous, to the pointthat different part of the nest can hardly been distinguished; galleries show atypical length-scale and there are no chambers or other specialised structures(see Fig. 1). As a consequence these nests resemble a continuum substratewith a coherent morphology, which makes them particularly suitable to betreated as the result of an accretive growth. Our task is to identify a relevantand minimal “stigmergic function” that summarizes the results of the interac-tions of termites with the building material (i.e. the local rules of growth andremodelling of the built structure) and which is sufficient to produce structurescomparable to the observed ones.The actual construction process is difficult to observe, as termites con-centrate their building activity in short and unpredictable bursts. The onlyavailable descriptions of termite behaviour during nest expansions come fromthe laboratory experiments of [21] (from the same author see also [20]). Jonesobserved that construction always happens by deposition of material on theedge of existing walls, and that these edges progressively bend, split, and mergein order to form an intricate matrix, where the local surface of the walls oftenresembles to a saddle. These observations suggest that the local curvature ofthe wall surface focuses the growth activity. The idea that curvature mightguide termite aggregation and building behaviour was investigated much morerecently in experiments by [7] who showed the existence of a positive correlationbetween termites activity and surface curvature in
Macrotermes michaelseni termites (an African termite species not too closely related to
Nasutitermes ).The observation that termite activity is driven by curvature does not excludethat cement pheromones or inter-individual interactions [14, 17] are involvedin the nest construction. Simply we suggest that either these quantities aredirectly affected by nest curvature, or they co-vary with it.Consistently with these observations, we develop a model of surface evo-lution that depends directly on surface curvature. The model is reminiscent cm 1 cm Figure 1: Fragments of nests built by
Nasutitermes walkeri (left) and
Nasutitermesephratae (right). of gradient growth processes, where the local curvature of the surface con-centrates gradients and fuels instabilities. In these processes, the instability isstrong and results in separate branches that compete together [25], and cannotfuse [10]. Here the nests surfaces are tightly interconnected so we have to finda new way of writing this growth, allowing both separation and fusion.In section 2 we derive our minimal model and perform the linear analysisof our growth equation. In section 3 we briefly describe the numerical toolsthat allow to solve the full nonlinear equation and in section 4 we present theresults of numerical simulations. In section 5 we compare our simulated nestto the tomographies of real nest fragments collected in the field. Finally insection 6 we discuss the results of the present study and illustrate the futureapplications of our model.
A large set of growth phenomena can be identified with the dynamics at the in-terface that divides the growing substrate from the medium that supports thegrowth. One is then confronted with a two-phases problem where boundaryconditions change in time and take the form of a curve or a surface, dependingon whether the geometry of the problem is highly confined in one direction or ully three dimensional like our nests. Whenever this curve or surface largelydeviates from a planar front, it is quite challenging to treat the problem analyt-ically, and predictions can be done only locally. Implementing a time evolvingboundary condition in three dimensions is also difficult and would certainlyrequire additional computing resources compared to a regular domain. As analternative approach, here we adopt a phase-field model [12], or diffuse inter-face method, which consists in replacing the substrate and the medium with acontinuous phase endowed with a scalar field f which is defined everywhere inthe considered volume and determines the presence of the first or the secondentity according to a threshold value. The frontier between the two phasescan then be easily retrieved as an iso-contour of f in the considered volume.As long as f is continuous, the geometry of the frontier can be also obtainedfrom f via differentiation [15], the growth model is then extremely simple andconsists of only one differential equation for the field f . The growth equationshould be nonlinear because we want the growth process to saturate at somepoint, and we also want to impose that new walls are only built adjacent topre-existing ones, i.e. we want to describe an accretion process. We define ourscalar field f in the interval [0 , f = f , for example f = 0 .
5. We will thensay that a point x is inside the wall whenever f ( x ) > f and outside otherwise.If the nest growth is a function of its own shape, the growth equation musttake the form: ∂f∂t = A ( f ) (1)where the operator A is meant to mimic the main features of the growth processas they are inferred from nest observations, and we construct it as follows: A ( f ) = − f α (1 − f ) β · d ( ∇ · ˆ n ) − ∆( ∇ · ˆ n ) (2)where the constant d is strictly positive and the equation is presented in anon dimensional form. The unit vector ˆ n is the normal to the local iso-surfaceand is given by ˆ n = ∇ f / |∇ f | , for example ˆ n ( x ) = ( ∇ f / |∇ f | ) | x will be thenormal to the iso-surface f ( x ) = f ( x ) in the point x .The first term of the equation represents the actual growth and can bedecomposed in two factors: a differential operator − d ( ∇ · ˆ n ) and a non linearkernel f α (1 − f ) β . The first factor is precisely the mean curvature of the localiso-surface times a constant, and is intended to mimic how local curvatureenhances the building or digging activity: the field f increaseas (decreases)whenever the local curvature is positive (negative). The second factor vanishesat f = 0 , f = f = α/ ( α + β ) and is meant to damp anyvariation of f far from the isosurface f = f which is the only accessible space or our walking agents: there can’t be any growth neither in the empty spacefar from the pre-existent nest, nor deep inside a pre-existent wall.The second term in equation (2) is proportional to the diffusion of the meancurvature, which consequently compensates the growth term whenever thelocal curvature is too high. Technically this is a dissipative term which providesa cutoff at the highest wavelength and is necessary to produce a meaningfulpattern as is highlighted below in the linear analysis. Phenomenologically thereare several reasons to introduce this term. On one side, the growth of a wallhappens (field and experimental observations) by the sequential deposition offaecal pellets and soil/sands particles of finite size, thus we can expect theexistence of a minimum scale in the spatial structure. On the other side onecan expect that, while ruled, the contributions of thousands of agents must sumup with some noise whose average effect likely produces some diffusion. Finallyan explicit justification can be found in the experimental observations of [22]which include smoothing process among the tasks performed by individualworkers.Now that a growth equation is defined we want to make a further ap-proximation. In fact the curvature operator ∇ · ˆ n is highly non-analytic (andnonlinear) while we would like our equation to be as simple as possible forboth analytical and numerical purposes. In the hypothesis that |∇ f | is almostconstant in the nearby of f = 0 .
5, one can approximate ∇ · ˆ n ∝ ∆ f . In ad-dition we take α = β = 1, which fixes the maximum of the nonlinear kerneland the contour of the actual nest at f = 0 .
5. Thus, our simplified growthequation reads: ∂f∂t = − f (1 − f ) d ∆ f − ∆ f, (3)where ∆ is the standard laplacian or diffusion operator and ∆ is the bi-laplacian operator which commonly enters in elasticity problems [39] and wasrecently included in a phenomenological model for vegetation patterns [33].One observes that ∆ f appears here with the opposite sign of a common dif-fusion equation, that is the first term in equation (3) is anti-diffusive, whilethe second term is a diffusive term of higher order or a hyper-diffusion. Notethat the difference in the derivative order of the two operators is the key in-gredient for pattern formation as it will be shown below in the linear analysis.Finally one may remark that equation (3) does not conserve the total mass (orvolume) because of the non linear pre-factor f (1 − f ). We stress that thesetwo elements are necessary to produce an organized structure like our termitenests as already observed by [11]. .1 Linear analysis In its simplified form (3) the growth equation can be easily linearised around f = 0 . f = f +ˆ f exp [ i k · x + σt ] which give the following dispersion relation: σ ( k ) = ˜ f dk − k (4)where for simplicity we wrote | k | and ˜ f = f (1 − f ) = 0 .
25. The disper-sion relation is shown in Fig. S1. The parameter d is strictly positive whichmeans that σ ( k ) changes sign at k = ( ˜ f d ) / and has a local maximum at k max = ( ˜ f d/ / : equation (3) is linearly unstable with the most unstablelength at λ max = 2 π/k max while all the small scale disturbances beyond thecutoff length λ = 2 π/k are damped. Note that k > d , thus the system is formally always unstable. Nonetheless a threshold isfixed by the smallest k which can appear in a given domain of size L which is k = 2 π/L . The instability threshold d c is then fixed by L the size of the do-main, with d c = (2 π/L ) / ˜ f being the marginal condition. Interestingly in twodimensions the same eigenvalues problem arises when studying the bucklingmodes of a thin plate under symmetric in-plane compression [40]. The station-ary two-dimensional version of equation (3) coincides with a particular case ofFöppl–von Kármán equations and were discovered already more than one cen-tury ago [13]. Nonetheless a very similar problem has gathered new attentionin the last decades to explain delamination patterns occurring in multi-layeredmaterial where the coating material does not have the same elastic propertiesof the underlying substrate [30, 5, 1]. Notably a similar mechanism was alsoinvoked to explain the formation of Miura-ori folding patterns that approxi-mately appears in nature, for example in insects wings and blooming leaves[27]. In most cases emerging patterns significantly differ from the planar so-lution (i.e. the solution we chose at the beginning of this section) and likelyarise from secondary instabilities [2] that rather select a superposition of pla-nar waves which reflects the invariance of their orientation at the instabilityonset. We simulate equation (3) on a cubic grid with a finite difference code writ-ten in
Python and automatically parallelised with the
Numba compiler. Thelargest domain we considered is 256 x 256 x 256 where a growth simulationcan be accomplished within about 5000 cpu hours. Both laplacian and bi-laplacian operators are implemented with a first order explicit Euler scheme. f for a simulation starting witha noisy half sphere at the bottom of the simulation domain. Boundary conditions areperiodic at the lateral boundaries, and f = const at the top and bottom boundaries.The time is normalised by the growth rate of the most unstable mode in Fig. S1. Two different boundary conditions were considered: in the first configurationboundary conditions are triply periodic while in the second one lateral facesof the simulation box are still periodic but Dirichlet boundary conditions (i.e. f = const ) are imposed on bottom and top faces. Initial conditions are alsovaried. In the simplest case the whole simulation domain is initiated with whitenoise, mainly to characterize the long term (or quasi-stationary) solutions ofthe equation within the bulk. Alternatively we initiate simulations with partof the volume copied from a previous simulation where a distinct pattern hasalready appeared, while the remaining part of the domain is set to zero. Inthis case we also forbid any evolution for all the points that are initially abovea certain threshold which maintains the initial seed unchanged. This configu-ration is intended to mimic the expansion process of a pre-existing nest, or anew one which starts from an arbitrary shaped support (in the following wewill refer to that as a nest-like seed). As a first result we report that our model reproduces a growth process. InFig. 2 we show the evolution of the field f for a simulation initiated witha nest-like seed (see previous section) with the shape of a half-sphere at thebottom of an empty domain ( f = 0). One observes that the initial seed grows,and progressively invades the empty space at a rate that is comparable withthat of the most unstable mode in Fig. S1. Even more interestingly we noticethat far from the interface between the empty space and the growing nest thefield f is essentially stationary (for example there is little rearrangement ofthe pre-existing nest). This indicates that equation (3) describes an accretionphenomenon happening at the interface between the growing walls and theempty space, even if we did not implement any freezing or aging effect and a) t=38 (b) t=40.1(c) t=42.2 (d) t=43.3 Figure 3: Evolution of the f = 0 . he field f could, in principle, change everywhere and at any time. Finallywe remark that the growing interface frequently branches and merges as it isshown in Fig. 3. Coherently with the linear analysis, merging and branchingepisodes cover a few time scales.Secondly we observe that in the growing pattern of Fig. 2 a characteristiclength-scale appears which is the typical distance between two walls. We reportthat a similar result was obtained for all the initial conditions and boundaryconditions we considered. In Fig. 4 (top) we show a two-dimensional sliceof the f field for two different simulations (a,b) where (a) has triply periodicboundary conditions and was initiated with white noise everywhere, while (b) isthe one already presented in Fig. 2. Both simulations were considered at largetime, namely at t = 64 time scales for simulation (a) and t = 93 time scales forsimulation (b). To track the existence of a characteristic length-scale in ourobjects we perform a three-dimensional fast Fourier transform (FFT) which isshown in Fig. S2. The dominant frequency of the FFT was then obtained andsuperimposed (squares) to the patterns of Fig. 4 as a comparison. We remarkthat the dominant length indicated by the FFT is consistent with the observedpatterns. Finally we note that the emerging lengths are fairly consistent withthe most unstable length indicated by linear analysis as shown in Fig. S1.Nonetheless the observed patterns significantly differ from the planar solutionswe consider in section 2.1, and appears to be modulated in multiple directions.To quantify this feature we compute the auto-correlation functions C f ( x ), C f ( y ), C f ( z ) in the three Cartesian directions (symbols) and the function C f ( | x | ) averaged over all directions (solid line), with all the distances rescaledby the typical length given by the FFT. The results are shown at the bottomof Fig. 4. First, one observes that in all the diagrams C f ( | x | ) shows a firstmaximum around d = 1 which is consistent with FFT analysis. Then weremark that in simulation (a) all the directional correlation functions presentthe same peak as C f ( | x | ) which indicates that the observed pattern is highlyisotropic, while in simulation (b) there is no peak in the vertical correlationfunction C f ( z ), which seems to indicate that there is low modulation in thevertical direction. Indeed, all the simulations we performed ultimately showa similar sheet-like anisotropy. This happens after a very slow process ofrearrangement which usually takes tens or even hundreds of time scales afterthe pattern has invaded the entire domain. Our interpretation is that thisparticular pattern is always selected by nonlinear interactions or secondaryinstabilities similar to those observed in elastic plates [27, 2], but differentinitial conditions can enhance or not this effect, determining at which timethis pattern takes over. Comparison with real nests
We collected nest fragments of two different species of
Nasutitermes , namely
N.walkeri from the Sydney area (Australia) and
N.ephratae from Guyana; twoof these are shown in Fig. 1. The walls of these nests are sometimes infra-millimetric thus a correct reconstruction of their structure was possible onlywith the use of a micro-CT scanner whose resolution was set to 0 . f that we consider in our model and simulations. We repeated our scale analysison two ct-scans of fragments belonging to real nests of the species N. walk-eri (c) and
N. ephratae (d) which are reported in Fig. 4. As for simulatednests, the FFT of real nests indicates a dominant length-scale which is consis-tent with the observed distance between walls (see Fig. S2). Looking at theauto-correlation functions one sees that fragment (c) is poorly modulated inthe horizontal direction, which is consistent with the observed pattern, whilefragment (d) appears fully isotropic. Note that for these fragments verticaland horizontal directions were completely arbitrary. Indeed different degreesof isotropy or different local structures can be encountered in different speciesand even within the same nest at different depths, as in the case of
N. ephratae [37]. Interestingly [21] reported that in
N. costalis , if the nest is cut open alonga plane parallel to the growth direction one observes a tall columnar structurewhose description highly resembles our simulation (b). As a summary ourmodel appears complex enough to reproduce the variety of structures encoun-tered in real nests as a response to different initial conditions, which suggeststhat sensitivity to initial conditions also matters in the real growth process. Asa second comparison we report that whenever the growing interface approachesa boundary where f is set to 0, the growing tips tend to connect forming asort of roof scattered with holes, as it is shown in Fig. 5. Waiting even longerthe holes tend to be filled and the connections between the roof and the restof the structure are progressively eroded. Very interestingly in all the speciesof arboreal Nasutitermes , nests are covered with a uniform thin crust whichin species like
N. ephratae happens also to be poorly connected to the bulkstructure [20, 37].In the present literature there is no comprehension of what could controlthe activation/inhibition of a construction stage, which is usually extremelyrapid compared to the lifetime of a colony [21]. Establishing a direct correspon-dence between real nests and simulations on this question may be premature,however we consider extremely encouraging that such a simple model incorpo-rates a boundary sensitivity that can mimic the real process of finishing up thenest, which likely involves several other external (weather conditions, material a) z (b) (c) (d) C f Figure 4: Top: Cross section of the nest volume for two simulations (a,b) and twoct-scans of real samples (c,d). The walls are in orange and the empty space is inblue. The black square indicates the peak of the 3D fast Fourier transform over thewhole volume. Bottom: Auto-correlation functions C f , each one corresponding tothe nest above it. Symbols correspond to the value of C f when computed only inone Cartesian direction, while the solid lines correspond to the average value of C f over all possible directions. 12igure 5: Effect of a f = 0 boundary on the nest growth. Left and center panelcorrespond to a 2D vertical cut and a 3D rendering of simulation (b) at t = 128.The top face of the simulation domain was maintained constant at f = 0. Oneobserves, that approaching the forbidden boundary f = 0 the simulated nest closesitself forming a roof which resembles the thin external layer of Nasutitermes nests(right). availability, etc) ad internal (size and health of the colony) variables, that havebeen completely ignored at this stage. Moreover, we could observe in the fieldthat whenever the external layer of a thriving colony is disrupted, termitesrapidly construct a patch that seals the internal structure from the exterior.This suggests that boundary conditions could mimic the maintaining of a nonexpanding stage.At this point we have already established a strong correspondence betweenour model and real nests, which are: the very convoluted structure (frequentmerging and branching), the existence of a typical scale, the sensitivity to ex-ternal constraints both in space (boundary conditions) and time (initial condi-tions). However all these are mostly qualitative comparisons. Below we definea protocol to make our comparison quantitative. The difficulty of this task isdue to the fact that our objects are fundamentally disordered so that a statis-tical approach is the only possible route. The choice of the relevant quantityto compare may also be discussed, but coherently with our observations of realnests and with our growth model, the local curvature is the best candidate.The local curvature of a 3D surface in one point P can be entirely defined bytwo principal curvatures k and k which are defined as the minimum and themaximum among the curvatures k of all the normal sections, i.e. the planarcurves given by the intersection between the surface and a normal plane. Giventhe osculating circle to such a curve, the value of k is the reciprocal of its radius H ) - Gaussian curvature (Γ). Eachcontour corresponds to the successive 10 percentile of points from more frequentarea (black) to less frequent area (light gray). Dots represent surface elements (oneover four dots is shown). The small square indicates the bin size for KDE.14 hile the sign is positive or negative depending on whether the curve turns inthe same or the opposite sense of the local normal vector ˆ n . Equivalently onecan define the mean curvature H (that we already encountered in sec. 2) andthe Gaussian curvature Γ: H = ( k + k ) / k k (5)At this stage one may simply keep in mind that the sign of Γ tell us if thesurface is more similar to a sphere (Γ > < H tells us if we are at the exterior ( H >
H < H = 0 characterize peculiar saddle-like surfaces which are known as minimalsurfaces where no interior or exterior can be defined, the two sides beingindistinguishable. To define our statistical sample we take advantage of thediscrete representation of our objects. As we defined them, both simulationsand real nests are scalar fields defined on a grid of voxels, from where we extractan iso-surface in the form of three-dimensional triangular meshes obtainedvia the Lewiner Marching Cube [26] algorithm implemented in Python . Oursample space will then be constituted by the set of all surface elements, towhich we associate two independent attributes which are the mean curvature H and the Gaussian curvature Γ. These are obtained as a function of the meshvertices coordinates as described in section S.II; importantly these coordinatesare priorily normalised by the typical length given by FFT in order to make thecomparisons consistent. In Fig. S3 we report a histogram of the frequency ofour two observables H and Γ for a collection of CT-scans of real nest fragments,simulations and synthetic data. The distribution of H is always peaked andquite symmetric around H = 0 while the values of Γ are more frequent atΓ < H and Γ at once.In Fig. 6 we have reported the cumulative Gaussian kernel density estima-tion KDE (see equation S2) of the surface elements in the space ( H, Γ). Theregions between two contours indicate where ten percent of the total countsoccur, from the most frequent area, in black, to the most rare area, in lightgray. In the top and central rows we have reported the two simulations (a-b)and two fragments of real nests (c-d) already analysed in Fig. 4. At the bottomwe considered two synthetic benchmark surfaces which are a random surface(e) obtained from a volume of white noise which has been band-pass filteredaround a certain length, and a gyroid (f) which is a minimal surface, i.e H = 0everywhere. One recognises a strong resemblance in the distribution betweensimulation (b) and nest (d). Compared to all the other diagrams one observesthat the curvature distribution is pushed toward the center H = Γ = 0 and its hape is squeezed in the direction of the Γ axis. This distinction is coherentwith the one we have made in Fig. 4 in term of isotropy/anisotropy, sincein simulation (b) and nest (d) we observe an higher proportion of quasi-flatregions which is consistent with higher anisotropy. Conversely simulation (a)and nest (c) are more isotropic and share a broader distribution along the Γaxis which indicates an higher abundance of saddle-shaped regions and is closerto the curvature distribution of a gyroid (f). One also remarks that in nest(c) curvature distribution is more spread along the H axis, and it is extremelysimilar to the curvature distribution of the random surface (e). These resultssuggest that the construction process is intrinsically noisy and that probabilis-tic instead of deterministic rules should control the building behaviour [19, 8].Interestingly our deterministic model manages to capture this feature whenwe let the system to invade the empty space from a growing boundary as insimulation (b) but the noise is partially removed if we assign a random ini-tial condition to all the volume and let the system evolve for long time as insimulation (a). Construction, as performed by humans or other animals, usually involves theinteraction between a building agent and some building material under theguidance of innate or learned rules of building. In contrast, growth - as itis known in physics or biology - is an autonomous process where a substrategrows by itself, depending on its interaction with the external environment.Arboreal nests built by
Nasutitermes termites lie in the middle between thesetwo categories because they originate from the cooperation of thousands ofindividual workers but the interactions between them is largely guided bythe shape of the substrate itself which makes the nest construction a self-organized process. This justifies our approach of formulating a model wherethe agents’ behaviour is incorporated in a classical accretion process where thenest growth is focused in the region of high curvature. In spite of its simplicity,our model is sufficient to reproduce some of the global and local properties ofarboreal
Nasutitermes nests, namely the intricate structure, the existence ofa typical length-scale, the sensitivity to boundary and initial conditions, andthe abundance of saddle-shaped structure of zero mean curvature. We modelthe nest as a continuous phase endowed with a scalar field f whose intensitydelimits the nest walls or cavities according to a threshold value. The nestgrowth can then be described by a single nonlinear equation that includes agrowing term proportional to local curvature and a dissipative term which issensitive to small spatial scales, their relative importance being fixed by one djustable parameter d .First we report that when d is large enough our minimal equation suc-cessfully reproduces a growth process where the growing substrate progres-sively invades the available space, through a rich dynamics of branching andmerging which produces intricate structures similar to Nasutitermes nests. Westress here that while branching is commonly observed in many gradient drivengrowth phenomena both in inert [31] and biological matter [23, 9] merging isimpossible in such simple models [25]. Only recently when the gradient is rein-troduced on both sides some merging was observed [6]. With regard to ourmodel, we report that merging occurs only in 3D, while in 2D the symmetry f → (1 − f ) must be broken in order to observe it. We also report that, farfrom the front of freshly grown substrate, rearrangement of the interface isnegligible, although it is permitted. This is consistent with the empirical ob-servations by [21, 38] which report that in Nasutitermes , nest walls are barelyrearranged after the initial construction.All the simulated patterns clearly show the emergence of a characteristiclength which corresponds to the one predicted by linear analysis and scales as1 / √ d . This feature constitutes a strong similarity with nest samples collectedin the field which show a characteristic distance between neighboring walls.In real nests such a distance is comparable with termites body size, whichsuggest that termites may use their body as a template as recently proposedfor ants [24]. However our model shows that there is no need to enforce theuse of a template to observe the appearance of a typical distance. Interest-ingly at the linear order our problem is analogous to the buckling instabilityof a thin elastic plate under a symmetric in-plane compression of magnitude d , when d is large enough the plate buckles and the compression energy isreleased in the form of bending energy of the plate. This problem was origi-nally addressed more than one century ago but has received renewed attentionin recent studies that showed how in this case the bending patterns usuallydiffer from the planar waves predicted by linear analysis [2]. Instead richerpatterns appear which are the combination of multiple planar waves orientedin different directions, because of secondary instabilities and nonlinear inter-actions. A direct comparison between our model and the thin plate problemis not possible because of the different dimensionality and initial conditions.Yet, it is interesting to notice how, similarly to the thin plates problem, inour simulations we never obtain the planar solution but rather observe trulythree-dimensional patterns with a variable degree of isotropy depending onthe initial conditions. As a general rule our equation tends to select structureswhich have a narrow wavelength in one direction and are slightly modulatedin the two others. Incidentally we observe that nest fragments from differenttermites species, or even within the same nest, also show a distinct variability n term of isotropy, thus we speculate that initial conditions should be relevantin the growth of real nests as well, but further investigations will be necessaryin this regard. Generated patterns are also very sensitive to boundary con-ditions and in particular we observe that forbidding the nest expansion (e.g.imposing f = 0) at one boundary induces the formation of a uniform layer thatisolates the nest from the forbidden boundary. Interestingly all Nasutitermes nests are also covered with an external thin layer which seals the nest fromthe outside world and it is constantly repaired whenever damaged by wet-ting/drying episodes or other external mechanical stresses, as it was observedboth in the field and in the laboratory. A forbidden boundary does not nec-essarily correspond to the existence of a solid boundary but could correspondfor instance to a threshold-level of pheromone emanating from the colony, asin the recent modelling study by [29]. In any case we find relevant that oursimple construction model responds in a consistent way in term of expressedmorphology. Future empirical studies would be needed to identify the releaserstimuli that trigger the initiation and termination of nest expansion.Both the real nests and the simulated patterns include many saddle-shapedregions where the mean curvature is zero. This is also shown by our quantita-tive comparison based on the curvature distribution, that maps our objects toa specific distribution in the two-dimensional space of mean ( H ) and Gaussiancurvature (Γ). Both the simulated and scanned surfaces show large portions ofzero mean curvature and non positive Gaussian curvature, which corroboratesthe hypothesis that curvature is the main driving mechanism in construction.In fact if growth happens in the region of high mean curvature, we expect thata stationary or quasi-stationary configuration should correspond to one wherelocal mean curvature is low or null. Consistently, whenever the structure isnon trivial (e.g. different from a plane), highly connected and isotropic, wefind that Gaussian curvature is negative Γ <
0. Alternatively we observe bothin simulations and real nests the existence of sheet-like structures of the type( H ∼
0, Γ .
0) which corresponds to the emergence of a clear anisotropy.We also remark that curvature distribution of isotropic nests resembles thatof a random synthetic surface with a typical length-scale, which confirms theabsence of a global predetermined architecture. This is consistent with a con-struction process driven primarily by local information, as we implemented inour model.While nest structure is primarily random on the local scale, we alreadyobserved that nest fragments can have different degrees of internal anisotropyand distinctive features such as the external envelope. Due to the relativelysmall size of our samples, we do not have a detailed characterisation of nest-topology and connectivity on the global scale and the question remains open,if the observed topology and connectivity are an emerging property of a same ocal construction mechanism - as in the present model - or they must beenforced through additional global rules. On another side, we plan to deeperinvestigate the reasons that trigger the beginning or the end of a constructionstage. Previous literature is all but conclusive in this regard but we expectthat activation and inhibition should be related to the colony size, which canbe possibly translated in the abundance of some auxiliary field, like relativehumidity, CO or a pheromone. To this aim we will couple our model withsome auxiliary field that is constrained to diffuse or convect in the nest galleriesand explore whether we can observe significant time modulation in the nestexpansion.In the general context of accretion processes we believe that our equationmay be adapted to other phenomena where a growth rate sensitive to curvatureand a local smoothing process coexist, which is encouraged by the simplicityof the equation and the richness of the expressed patterns. While our modelwas primarily formulated to describe the morphogenesis of arboreal termitenests, saddle-shaped structures with a characteristic scale are also observed inother biological systems, ranging from the nests of other insects not directlyrelated to termites (e.g. Lasius fuliginosus ), to the micro-structure of butterflywings [28]. It is possible that our model could provide insight also on theformation of these other structures. The scientific literature on social insectshas often focused on agent-based models as a tool to describe the behaviourof individual insects and the emerging organisation of the colony. We areconfident that simple “stigmergic” models like the one developed here havea great potential to explain not only how structures are built, but also howtransitions between different nest plans can be triggered by small changesin the building parameters or environmental conditions. In relation to realnests, such changes could correspond to transitions between nest plans thatemerge across phylogenies and could also reflect adaptive changes in responseto different habitats.
References [1] B. Audoly. Stability of Straight Delamination Blisters.
Physical Re-view Letters , 83(20):4124–4127, Nov. 1999. ISSN 0031-9007, 1079-7114.doi: 10.1103/PhysRevLett.83.4124. URL https://link.aps.org/doi/10.1103/PhysRevLett.83.4124 .[2] B. Audoly and A. Boudaoud. Buckling of a stiff film bound to a compli-ant substrate—Part I:: Formulation, linear stability of cylindrical pat-terns, secondary bifurcations.
Journal of the Mechanics and Physics f Solids , 56(7):2401 – 2421, 2008. ISSN 0022-5096. doi: https://doi.org/10.1016/j.jmps.2008.03.003. URL .[3] P. Bardunias and N.-Y. Su. Dead reckoning in tunnel propagation of theformosan subterranean termite (isoptera: Rhinotermitidae). Annals ofthe Entomological Society of America , 102(1):158–165, 2009.[4] P. Bourgine and A. Lesne.
Morphogenesis: origins of patterns and shapes .Springer Science & Business Media, 2010.[5] N. Bowden, S. Brittain, A. G. Evans, J. W. Hutchinson, and G. M.Whitesides. Spontaneous formation of ordered structures in thin filmsof metals supported on an elastomeric polymer.
Nature , 393(6681):146–149, May 1998. ISSN 0028-0836, 1476-4687. doi: 10.1038/30193. URL .[6] A. Budek, K. Kwiatkowski, and P. Szymczak. Effect of mobility ratioon interaction between the fingers in unstable growth processes.
Phys.Rev. E , 96:042218, Oct 2017. doi: 10.1103/PhysRevE.96.042218. URL https://link.aps.org/doi/10.1103/PhysRevE.96.042218 .[7] D. S. Calovi, P. Bardunias, N. Carey, J. Scott Turner, R. Nagpal, andJ. Werfel. Surface curvature guides early construction activity in mound-building termites.
Philosophical Transactions of the Royal Society B:Biological Sciences , 374(1774):20180374, 2019. doi: 10.1098/rstb.2018.0374. URL https://royalsocietypublishing.org/doi/abs/10.1098/rstb.2018.0374 .[8] S. Camazine, J.-L. Deneubourg, N. R. Franks, J. Sneyd, G. Theraulaz,and E. Bonabeau.
Self-organization in biological systems . Princeton stud-ies in complexity. Princeton University Press, Princeton, N.J, 2001. ISBN0691012113 (cloth : alk. paper). URL .[9] R. Clément, P. Blanc, B. Mauroy, V. Sapin, and S. Douady. Shape Self-Regulation in Early Lung Morphogenesis.
PLoS ONE , 7(5):e36925, May2012. ISSN 1932-6203. doi: 10.1371/journal.pone.0036925. URL https://dx.plos.org/10.1371/journal.pone.0036925 .[10] Y. Couder, J. Maurer, R. González-Cinca, and A. Hernández-Machado.Side-branch growth in two-dimensional dendrites. i. experiments.
Physicalreview. E, Statistical, nonlinear, and soft matter physics , 71:031602, 032005. doi: 10.1103/PhysRevE.71.031602.
11] J.-L. Deneubourg. Application de l’ordre par fluctuations a la descrip-tion de certaines étapes de la construction du nid chez les Termites.
In-sectes Sociaux , 24(2):117–130, June 1977. ISSN 0020-1812, 1420-9098.doi: 10.1007/BF02227166. URL http://link.springer.com/10.1007/BF02227166 .[12] A. Fasano and M. Primicerio.
Fix, G.J. in Free Boundary Problems:Theory and Applications , volume II. Pitman Advanced Pub. Program,1983. ISBN 0273085905,9780273085904.[13] A. Föppl. Vorlesungen uber technische mechanik.
B.G. Teubner, Leipzig,Germany , page 132, 1907.[14] D. Fouquet, A. M. Costa-Leonardo, R. Fournier, S. Blanco, and C. Jost.Coordination of construction behavior in the termite Procornitermesaraujoi: structure is a stronger stimulus than volatile marking.
In-sectes Sociaux , 61(3):253–264, Aug. 2014. ISSN 0020-1812, 1420-9098.doi: 10.1007/s00040-014-0350-x. URL http://link.springer.com/10.1007/s00040-014-0350-x .[15] R. Goldman. Curvature formulas for implicit curves and surfaces.
Com-puter Aided Geometric Design , 22(7):632–658, Oct. 2005. ISSN 01678396.doi: 10.1016/j.cagd.2005.06.005. URL https://linkinghub.elsevier.com/retrieve/pii/S0167839605000737 .[16] P. P. Grasse. La reconstruction du nid et les coordinations interindi-viduelles chez bellicositermes natalensis et cubitermes sp. la théorie dela stigmergie: Essai d’interprétation du comportement des termites con-structeurs.
Insectes Sociaux , page 40, 1959.[17] B. Green, P. Bardunias, J. S. Turner, R. Nagpal, and J. Werfel. Exca-vation and aggregation as organizing factors in de novo construction bymound-building termites.
Proceedings of the Royal Society B: BiologicalSciences , 284(1856):20162730, 2017. doi: 10.1098/rspb.2016.2730.[18] W. Hangartner. Carbon dioxide, a releaser for digging behavior in solenop-sis geminata (hymenoptera: Formicidae).
Psyche , 76:58–67, 1969. URL http://psyche.entclub.org/76/76-058.html .[19] E. Invernizzi and G. D. Ruxton. Deconstructing collective building insocial insects: implications for ecological adaptation and evolution.
In-sectes Sociaux , 66(4):507–518, Nov. 2019. ISSN 0020-1812, 1420-9098.doi: 10.1007/s00040-019-00719-7. URL http://link.springer.com/10.1007/s00040-019-00719-7 .
20] R. J. Jones.
The building behavior of arboreal Nasutitermes . PhD thesis,Harvard University, 1977.[21] R. J. Jones. Expansion of the nest of nasutitermes costalis.
InsectesSociaux , 26(4):322–342, Dec. 1979. ISSN 1420-9098. doi: 10.1007/BF02223552. URL https://doi.org/10.1007/BF02223552 .[22] R. J. Jones. Gallery construction by nasutitermes costalis: Polyethismand the behavior of individuals.
Insectes Sociaux , 27(1):5–28, Mar. 1980.ISSN 0020-1812, 1420-9098. doi: 10.1007/BF02224518. URL http://link.springer.com/10.1007/BF02224518 .[23] J. A. Kaandorp. Morphological analysis of growth forms of branch-ing marine sessile organisms along environmental gradients.
MarineBiology , 134(2):295–306, July 1999. ISSN 0025-3162, 1432-1793. doi:10.1007/s002270050547. URL http://link.springer.com/10.1007/s002270050547 .[24] A. Khuong, J. Gautrais, A. Perna, C. Sbaï, M. Combe, P. Kuntz, C. Jost,and G. Theraulaz. Stigmergic construction and topochemical informa-tion shape ant nest architecture.
Proceedings of the National Academy ofSciences , 113(5):1303–1308, 2016. ISSN 0027-8424. doi: 10.1073/pnas.1509829113. URL .[25] E. Lajeunesse and Y. Couder. On the tip-splitting instability of viscousfingers.
Journal of Fluid Mechanics , 419:125 – 149, 09 2000. doi: 10.1017/S0022112000001324.[26] T. Lewiner, H. Lopes, A. W. Vieira, and G. Tavares. Efficient implemen-tation of marching cubes’ cases with topological guarantees.
Journal ofGraphics Tools , 8(2):1–15, 2003. doi: 10.1080/10867651.2003.10487582.URL https://doi.org/10.1080/10867651.2003.10487582 .[27] L. Mahadevan. Self-Organized Origami.
Science , 307(5716):1740–1740, Mar. 2005. ISSN 0036-8075, 1095-9203. doi: 10.1126/science.1105169. URL .[28] K. Michielsen and D. Stavenga. Gyroid cuticular structures in butterflywing scales: biological photonic crystals.
Journal of The Royal SocietyInterface , 5(18):85–94, 2008.
29] S. A. Ocko, A. Heyde, and L. Mahadevan. Morphogenesis of termitemounds.
Proceedings of the National Academy of Sciences , 116(9):3379–3384, 2019. ISSN 0027-8424. doi: 10.1073/pnas.1818759116. URL .[30] M. Ortiz and G. Gioia. The morphology and folding patterns of buckling-driven thin-film blisters.
Journal of the Mechanics and Physics ofSolids , 42(3):531 – 559, 1994. ISSN 0022-5096. doi: https://doi.org/10.1016/0022-5096(94)90030-2. URL .[31] P. Pelcé.
Théorie des formes de croissance . CNRS editions, 2000.[32] A. Perna and G. Theraulaz. When social behaviour is moulded in clay: ongrowth and form of social insect nests.
Journal of Experimental Biology ,220(1):83–91, 2017.[33] D. Ruiz-Reynés, F. Schönsberg, E. Hernández-García, and D. Gomila.A simple model for pattern formation in clonal-growth plants. arXiv:1908.04603 [physics, q-bio] , Aug. 2019. URL http://arxiv.org/abs/1908.04603 . arXiv: 1908.04603.[34] P. G. Saffman and G. I. Taylor. The penetration of a fluid into aporous medium or hele-shaw cell containing a more viscous liquid.
Pro-ceedings of the Royal Society of London. Series A. Mathematical andPhysical Sciences , 245(1242):312–329, 1958. doi: 10.1098/rspa.1958.0085.URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1958.0085 .[35] G. Theraulaz and E. Bonabeau. A Brief History of Stigmergy.
Arti-ficial Life , 5(2):97–116, Apr. 1999. ISSN 1064-5462, 1530-9185. doi:10.1162/106454699568700. URL .[36] D. W. Thompson et al. On growth and form.
On growth and form. , 1942.[37] B. L. Thorne. Differences in nest architecture between the neotropi-cal arboreal termites nasutitermes corniger and nasutitermes ephratae(isoptera: Termitidae).
Psyche: A Journal of Entomology , 87, 01 1980.doi: 10.1155/1980/12305.
38] B. L. Thorne, M. S. Collins, and K. A. Bjorndal. Architecture andNutrient Analysis of Arboreal Carton Nests of Two Neotropical Nasu-titermes Species (Isoptera: Termitidae), with Notes on Embedded Nod-ules.
The Florida Entomologist , 79(1):27, Mar. 1996. ISSN 00154040.doi: 10.2307/3495751. URL .[39] S. P. Timoshenko and J. N. Goodier.
Theory of elasticity . McGraw-Hill,New York, 1970. ISBN 978-0-07-064270-6 978-0-07-085805-3 978-0-07-064720-6. OCLC: 45153.[40] S. P. Timoshenko, J. M. Gere, and W. Prager. Theory of Elastic Stability,Second Edition.
Journal of Applied Mechanics , 29(1):220, Jan 1962. doi:10.1115/1.3636481.[41] E. Toffin, D. Di Paolo, A. Campo, C. Detrain, and J.-L. Deneubourg.Shape transition during nest digging in ants.
Proc. Natl. Acad. Sci. USA ,106(44):18616–18620, Nov 2009. ISSN 1091-6490. doi: 10.1073/pnas.0902685106. upplementary material: a growth model driven bycurvature reproduces geometric features of arborealtermite nests G. Facchini , A. Lazarescu , A. Perna , and S. Douady Life Sciences Department, University of Roehampton, London, UK Institut de Recherche en Mathématique et Physique, UC Louvain,Louvain-la-Neuve, Belgium Laboratoire Matière et Systèmes Complexe, Université ParisDiderot, Paris, France
Figure S1: Dispersion relation S1 for small-amplitude perturbations, here ˜ f = 0 . d = 4. The solid line indicates the most unstable length scale while dashedand dotted lines indicate the typical length that appears in simulations (a) and (b)reported in figure 4. Black dots refer to lengths that can fit in our simulation box.1 a r X i v : . [ n li n . AO ] F e b
2 4 6 8 10 | k | FF T
0 2 4 6 8 10 | k | FF T
0 2 4 6 8 10 | k | FF T
0 2 4 6 8 10 | k | FF T Figure S2: Amplitude spectrum of the fast fourier transform (FFT) of the scalarfield associated to simulations (a,b) and CT-scans of real nests (c,d) shown at the topof figure 4. The local maxima correspond to the emergent characteristic length scaleswhich are reported in figure 4. Wave numbers k are normalized by the characteristiclength scale. S.I Linear and scale analysis
In its simplified form the growth equation (3) can be easily linearised around f = 0 . f = f + ˆ f exp [ i k · x + σt ]which gives the following dispersion relation: σ ( k ) = ˜ f dk − k (S1)where for simplicity we wrote | k | = k and ˜ f = f (1 − f ). The dispersion relationis shown in figure S1. As d is strictly positive σ ( k ) changes sign in k = ˜ f d ) / and has a local maximum at k max = ( ˜ f d/ / : equation (3) is linearly unstablewith the most unstable length scale at λ max = 2 π/k max while all the small scaledisturbances beyond the cutoff length λ = 2 π/k are damped. Note that k > d , thus the system is formally always unstable. Nonetheless athreshold is fixed by the smallest k which can appear in a given domain of size L which is k = 2 π/L . The instability threshold d c is then fixed by L the size of thedomain, with d c = (2 π/L ) / ˜ f the marginal condition. In figure S2 we report the fastFourier transform performed on the scalar field associated to two simulations andthe computer-tomography images of two nests collected from the field. One observesthe presence of a local maximum which denotes the emergence of a characteristiclength scale. 2 .II Segmentation and curvature statistical distri-bution In this section we explain how we obtain surfaces from a scalar field and construct thecurvature distributions shown in figure S3 and 6. The protocol includes five steps: (i)smoothing and segmentation of the original scalar field, (ii) isosurface extraction, (iii)smoothing of the surface and (iv) computation of curvature, (v) statistical analysisof curvature distribution.The first step only concerns nest samples collected in the field, in fact CT imagesmust be processed to remove some of the small scale roughness of the nests and thenoise introduced by the CT-scanner. First, we filter the original greyscale imagesby eliminating 10 to 20% of the highest Fourier components, which convenientlyreduces the size of our data and removes isolated noise voxels due to acquisitionnoise or small spurious particles inside the nests (dust, dead termites, etc.). Wefollow this operation by applying a Gaussian filter ( σ = 1), mainly to remove theripples introduced by the previous sharp low-pass filtering operation. After thisinitial blurring operation, we binarise our stack of images and we apply a secondGaussian filter ( σ = 1).This final smoothing operation allows having a smooth transition between wallsand galleries while having a uniform value f = 1 deep in the walls and f = 0 inthe middle of a gallery. Note that the threshold for segmentation must be chosen byhand and adjusted in order to preserve the topology of the original scanner which isinspected comparing slices of the original stack to the processed ones.The second step is the isosurface extraction, that results in the production ofa three-dimensional triangular mesh which represent the iso-contour surface of thescalar field at a given threshold f . We use the efficient Lewiner implementation ofthe Marching Cubes algorithm [ S1 ] provided by Python scikit-image library [ S4 ].This essentially performs an iteration over all the voxels of the considered volume inwhich a new surface element is drawn for all the voxels where there is at least onevertex with f > f and one with f < f .As a third step we compute the mean curvature H and the Gaussian curvatureΓ at each vertex of the mesh using the algorithm of [ S2 ]. Mean curvature basicallycorresponds to the sum of all edges connecting the considered vertex to its neighbors,while Gaussian curvature is essentially the sum of the angles between these edgesand those connecting the neighbor vertices weighted by the area of the correspondingface. Finally both H and Γ are averaged over the three vertices of a surface element.Prior to this computation the coordinates of all vertices are normalised by the typicallength given by the FFT of the original scalar field: with this choice a surface element3a)(b)(c)(d)(e)(f) H ΓFigure S3: Normalised frequency of mean curvature H (left) and Gaussian curvatureΓ (right) for different surfaces, (a): simulation initiated with white noise, (b): sim-ulation initiated with a nest-like seed, (c): nest fragment from Guyane- N.ephratae ,(d): nest fragment from New South Wales (AUS)-
N. walkeri , (e): pass-band filterednoise, (f): gyroid minimal surface. The 3D rendering at the top represent the surfaceof nest (c) coloured according to the local value of H (left) and Γ (right).4hose curvature is ( H, Γ) = (1 ,
1) corresponds to an element of spherical surface ofradius equal to the typical distance between walls. In addition, as we are using a first-neighbors algorithm, we expect the curvature computation to be extremely sensitiveto the mesh roughness which results from the discrete nature of the original scalarfield (sampled in discrete voxels). In order to remove this small scale noise (which isnot relevant to the actual characteristic scale of the structure), each triangular meshwas smoothed replacing each vertex coordinates with the average of the coordinatesof neighboring vertices. This step is crucial, since on one side we want to get rid of thissensitivity and on the other side we ultimately want surface curvature to characterizeour objects, so that we must avoid modifying the surface too much. Nonetheless, weshow that, by repeating the smoothing operation several times, curvature propertiesseem to have converged after at most 20 iterations, as it is shown in figure S4 andS5. Moreover, we observe that if we add spurious small scale noise to a smoothedmesh, the curvature distribution significantly changes, as it is shown in figure S6.The last steps consist in producing the curvature distribution of the analysedsurfaces. In figure S3 we report the normalised frequency of occurrence for H andΓ in the form of histograms for a sample of two different simulations (a,b), twodifferent nest fragments (c,d), and two benchmark synthetic surfaces (e,f) where(e) is a random surface obtained by segmentation of a volume of band-pass filteredwhite noise and (f) is a minimal surface (i.e. H = 0 everywhere) called gyroid.Subsequently, we present the same data in the two dimensional diagram ( H, Γ) andestimate the probability density ρ kde with a Gaussian kernel which is defined asbelow: ρ kde ( x ) = 1 nh n X i =1 e | x − x i | h , (S2)where x is one generic point in the plane ( H, Γ) while x i =1 ,n are all the surfaceelements. This tool is provided by the Python scipy library and the bandwidth h isestimated with the rule-of-thumb [ S3 ]. The cumulative probability density associatedwith ρ kde is shown in figure 6, S4, S5 and S6 where each contour corresponds to thesuccessive ten percentiles from the most frequent area (black) to the least frequentarea (grey). Dots correspond to the raw data before the kernel density estimationis realised. One observes that raw data are confined to remain below the parabolaΓ = H which is coherent with the definition of H and Γ given in equation (5). Notethat kernel density estimation smoothens this boundary which explains why contourlines (but no dots) extend to the forbidden region.5igure S4: Curvature distribution for a Nasutitermes walkeri nest sample, whereLaplacian smoothing algorithm was applied 0, 5, 10, and 20 times. The distributionremains relatively stable after 10 iterations.6igure S5: Curvature distribution for a
Nasutitermes ephratae a sample, whereLaplacian smoothing algorithm was applied 0, 5, 10, and 20 times. One observesthat the distribution remains relatively stable after 10 iterations.7igure S6: Top: Curvature distribution for synthetic noise where Laplacian smooth-ing algorithm was applied for 0 (left) and 20 (right) iterations. Bottom: Same asabove but a white noise was added to the initial surface. The addition of noiseinitially introduces a significant change in the curvature distribution (left top andbottom) but this can easily be removed by applying the Laplacian smoothing algo-rithm. 8 eferences [S1] T. Lewiner, H. Lopes, A. W. Vieira, and G. Tavares. Efficient implementation ofmarching cubes’ cases with topological guarantees.
Journal of Graphics Tools , 8(2):1–15, 2003. doi: 10.1080/10867651.2003.10487582. URL https://doi.org/10.1080/10867651.2003.10487582 .[S2] M. Meyer, M. Desbrun, P. Schröder, and A. H. Barr. Discrete Differential-Geometry Operators for Triangulated 2-Manifolds. In G. Farin, H.-C. Hege,D. Hoffman, C. R. Johnson, K. Polthier, H.-C. Hege, and K. Polthier, editors,
Visualization and Mathematics III , pages 35–57. Springer Berlin Heidelberg,Berlin, Heidelberg, 2003.[S3] D. Scott.
Multivariate Density Estimation: Theory, Practice, and Visualization .A Wiley-interscience publication. Wiley, 1992. ISBN 9780471547709. URL https://books.google.fr/books?id=7crCUS_F2ocC .[S4] S. van der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner,N. Yager, E. Gouillart, T. Yu, and the scikit-image contributors. scikit-image:image processing in Python.
PeerJ , 2:e453, 6 2014. ISSN 2167-8359. doi:10.7717/peerj.453. URL https://doi.org/10.7717/peerj.453https://doi.org/10.7717/peerj.453