Attraction, Dynamics, and Phase Transitions in Fire Ant Tower-Building

Many insect species, and even some vertebrates, assemble their bodies to form multi-functional materials that combine sensing, computation, and actuation. The tower-building behavior of red imported fire ants, Solenopsis invicta, presents a key example of this phenomenon of collective construction. While biological studies of collective construction focus on behavioral assays to measure the dynamics of formation and studies of swarm robotics focus on developing hardware that can assemble and interact, algorithms for designing such collective aggregations have been mostly overlooked. We address this gap by formulating an agent-based model for collective tower-building with a set of behavioral rules that incorporate local sensing of neighboring agents. We find that an attractive force makes tower building possible. Next, we explore the trade-offs between attraction and random motion to characterize the dynamics and phase transition of the tower building process. Lastly, we provide an optimization tool that may be used to design towers of specific shapes, mechanical loads, and dynamical properties, such as mechanical stability and mobility of the center of mass.


INTRODUCTION
Collective aggregation is a prevalent behavior among social animals, where many individuals cluster together while feeding, defending against predators, or as a thermoregulation strategy, effectively reducing the exposed surface area per individual. Examples of species that aggregate include vertebrates, such as penguins (Waters et al., 2012) and bats (Roverud and Chappell, 1991;Kerth, 2008) as well as insects, such as beetles (Deneubourg et al., 1990), ants Reynaert et al., 2006), and cockroaches (Ame et al., 2004;Jeanson et al., 2005). While these aggregations are often planar, eusocial insects, such as honey bees (Seeley, 2010;Kastberger et al., 2011), army ants (Franks, 1989), and fire ants (Mlot et al., 2011) extend this strategy and create three-dimensional assemblages. These self-assemblages are composed of a multitude of individuals who link their bodies, doing so without a global overseer and with limited cognitive abilities (Anderson et al., 2002).
Ants in particular are capable of a wide variety of self-assemblages and collective behavior. For example, ants of the genus Oecophylla build chains for gap crossing and during nest construction (Lioni et al., 2001). In addition, army ants are known for their construction of bivouacs (Franks, 1989), and are also capable of building bridges out of their bodies to cross gaps along a foraging trail (Reid et al., 2015). Finally, as we will discuss further, fire ants gather together to form rafts and towers when their habitat floods (Mlot et al., 2011(Mlot et al., , 2012Phonekeo et al., 2017).
The structures that these insects create are, in essence, autonomous materials that embed sensing, computation, and actuation. These properties are some of the long-standing aspirations in the fields of multi-functional materials and robotic materials (Şahin, 2004;Hughes et al., 2019). Self-assembling agents have already begun to inspire robotic applications (Bonabeau et al., 1999;Brambilla et al., 2013;Hamann, 2018). For example, Del Dottore et al. (2018) have described the concept of "growing robots, " which are systems of a large number of individual robots working together to mimic biological growth in plants or groups of molecules or cells. Other collective robots are directly inspired by eusocial insects, such as the S-bots (Şahin et al., 2002;Groß et al., 2006), which form chains to collectively move larger payloads, just like ants working together to move larger food (Buffin and Pratt, 2016). Also inspired by ants, Swissler and Rubenstein (2018) have developed robots with a new docking mechanism to form self-assembling structures. Another class of robots, inspired by termites (Werfel et al., 2014), build three-dimensional structures out of external building materials. Finally, the cube-shaped M-Blocks (Romanishin et al., 2015) construct aggregations out of their own bodies, using magnetism and angular momentum to climb on top of neighbors. These works represent examples from the emerging field of multi-agent robotic systems built out of many inexpensive individual robots and utilizing control strategies that may include redundancies to overcome individual malfunctions. While much of the focus in robotics has been on developing the hardware, the algorithmic development of assembling processes has often been overlooked. We address this gap by borrowing tools from computational material science and characterize the dynamics of 3-dimensional aggregation formation inspired by fire ant towers.
In nature, red imported fire ant (Solenopsis invicta) towers tend to occur in the event of flooding. Initially, fire ants gather together to form hydrophobic rafts (Mlot et al., 2011(Mlot et al., , 2012 to float above the water surface. When the rafts approach vegetation emerging from the surface, they may attach to the vegetation and form towers on top of their floating rafts, as pictured in Figure 1A. In a recent, Phonekeo et al. (2017) described an experimental assay of the tower-building process in fire ants. The experimental setup involved fire ants constructing towers around a vertical rod to represent the emergent vegetation. In their analysis, the authors propose four rules which allow ants to build towers: 1. Do not move if ants are on top of you. 2. If atop other ants, repeatedly move a short distance in a random direction. 3. Upon reaching available space adjacent to non-moving ants, stop and link with them. 4. The top layer of the tower is not stable unless there is a complete innermost ring of ants gripping each other around the rod.
Note that the "available space adjacent to non-moving ants" is primarily discussed by Phonekeo et al. (2017) in the context of a ring around the vertical rod or vegetation. We will take a more general definition of an available space in the present study, discussed below in section 2.1. The work of Phonekeo et al. (2017) shows an agreement between the resulting tower shapes in the long-timescale limit; however, it does not explore the time dynamics and parameter space systematically. This is what the present work aims to do, since local rules, such as these provide a systematic way of analyzing collective behavior through agent-based modeling, and importantly, they are directly implementable in swarm robotic systems. By simulating the behavior of individuals following a set of local rules, it is possible to investigate how local interactions between agents lead to global emergent behavior and explore the space of possible behavior beyond what is possible with experiments.
Modeling efforts of collective behavior using local behavioral rules include the boids model (or Vicsek model) (Reynolds, 1987;Vicsek et al., 1995), which simulates agents moving under attraction, repulsion, and alignment as well as more complicated models (Couzin and Krause, 2003;Mishra et al., 2012;Wilensky and Rand, 2015). However, boids-type models best describe the behavior of more sparse collectives, such as flocks of birds or schools of fish. To model ants building a tower, we must account for dense aggregations where the interaction range is limited to a short length scale, preferably defined by the size of an individual agent. Models of more dense collective assemblies include aggregation in slime mold based on chemical signal amplification (Levine et al., 1997;Umeda and Inouye, 1999), and nest building in wasps using an agent-based model in which swarms of builders deposit bricks and build up a nest (Theraulaz and Bonabeau, 1995;Bonabeau et al., 2000). Agentbased modeling has been successfully applied to studies of ant collective behavior as well (Dorigo et al., 2000) to modeling traffic organization in ant foraging (Goss et al., 1989;Couzin and Franks, 2003;Strömbom and Dussutour, 2018), bridge and chain formation (Lioni et al., 2001;Garnier et al., 2013), and trail clearing (Bochynek et al., 2017). However, for the present study we must consider both moving ants as in traffic organization and trail clearing, which climb the tower and form the shape, as well as stationary ants that support the structure as in bridge and chain formation. Based on the similarity to the aggregation of inanimate systems, such as colloids Vernerey et al., 2018), we reason that ant tower building would experience dynamic phase separation processes including nucleation (Vlasov, 2019), jamming (Bak, 1996), and ripening (Voorhees, 1985). These phase transitions are also observed at the thermodynamic transition between phases of matter, which have been studied experimentally (Panagiotou et al., 1984) as well as computationally (Rovere et al., 1990;Navarro and Fielding, 2015). Hence, we formulate an agent-based model with a set of behavioral rules that lead to aggregation and experience dynamical phase transitions.
Section 2 describes the details of the model we study in the present work and lays out the modifications to the local rules (presented above) that we introduce to achieve tower-building. In section 3, we explore the parameter space of the local rules to identify the impacts of each component: locking, unlocking, and FIGURE 1 | (A) Formation of towers by fire ants around vegetation, reproduced with permission from Phonekeo et al. (2017). (B) Schematic of computational model considered in this study, as described in section 2.2. Individual free agents (black) move to an adjacent square at every time step, while locked agents (red) remain stationary. The side view below the grid shows the 3-dimensionality of the model arena. The blue region represents the neighborhood of adjacent locations which an agent considers when searching for locked neighbors and for determining its attractive force. (C) Schematic representation of the three states an agent can assume in the model. Free agents move with constant velocity, locked agents stop to build towers, while covered agents cannot unlock.
attraction. We find that towers undergo a phase transition when varying the attraction parameter, and explore how this phase transition changes across various densities. Finally, we introduce an optimization algorithm to generate the largest possible tower for a given density of agents in the system. In section 4, we discuss the significance of the results and talk about implications for both the understanding of collective biological systems and the design of multi-agent robotic control strategies.

AGENT-BASED MODEL
We consider a system of N individual agents simulated to move in a L × L × ∞ arena, discretized into a cubic lattice made of voxels of volume ℓ × ℓ × ℓ. The volume of an individual agent is set to the be volume of a voxel, where ℓ ≡ 1. At each time step, an individual agent can move into one of its 26 neighboring voxels: 9 above, 9 below, and 8 on the same level. A schematic of agents moving within the arena is shown in Figure 1B. In the present work, we will not consider the effects of solid wall boundaries and will instead implement periodic boundaries. The horizontal plane of the arena, therefore, contains periodic boundary conditionswhen an agent leaves the right side of the arena, for example, it re-enters the left side. Periodic boundaries are also taken into account when distances between agents are calculated. The equations that define the periodic boundary conditions are given in (S1) and (S2) in Appendix 2. The vertical direction of the arena is semi-infinite, extending upward from a solid floor.
Agents move horizontally and climb up if the voxel they intended to move into is occupied by a locked agent. Note that the local rules described above, from Phonekeo et al. (2017), refer to agents "linking" with one another, while in the present work we will refer to an agent that stops to support tower building as "locked." Each pixel along the horizontal plane has an associated height equal to the number of locked agents on top of each other in that location. The free agents, therefore, are moving under 2-dimensional rules along the surface defined by locked agents, which is embedded in 3-dimensions. If an agent attempts to climb on top of neighbors to a voxel that is more than ℓ higher, it does not move at this time step. Agents may move down any distance but never move below the floor. Agents in the model may take on three different states, depicted in Figure 1C: free, locked, or covered. A free agent may move around the arena according to a specific set of behavioral rules with a constant velocity of one voxel per time step. All agents determine their intended movement before moving, and movement order is chosen randomly at each time step. To prevent two individuals from occupying the same position, if two free agents attempt to move to the same voxel, the second agent to arrive randomly chooses an unoccupied voxel horizontally adjacent to the target voxel. Locked agents are those which have decided to become a part of a tower and allow their neighbors to climb on top of them. We explore different schemes for the decision to "lock" as defined in sections 2.1 and 2.2. Covered agents are locked agents with at least one other agent on top of them. Each time step consists of first evaluating movement for all individuals and then evaluating locking decisions for all individuals based on their new configuration. We will not consider the effect of stability and assume that each agent has infinite strength to support neighbors.
It is likely that pheromones play a role in fire ant tower building, but for the present study, we consider whether this behavior can arise from solely physical proximity to neighbors.
Hence, an agent can sense which of its surroundings 26 voxels are occupied by another agent. This local model will allow for easier implementation by collective robotic systems, as it merely requires local sensing.
Unless specified otherwise, all simulations contain N = 1, 000 agents moving in a 100 × 100 × ∞ arena, corresponding to a density of ρ = N L 2 = 0.1. Based on preliminary simulations, we have chosen to evaluate each trial for 500,000 time steps, for which 97.8% of all simulations considered reached a steady state, where the largest tower size remained approximately constant (±5% of N) for at least the last 100,000 time steps of the simulation. Exceptions will be discussed below in section 3.2.

Diffusion-Limited Aggregation
We start by investigating whether a dynamic simulation of the proposed local rules above can lead to tower-building. As we are not considering effects of stability, we will ignore rule (iv) in the present study. We simulate the rules (i)-(iii) from Mlot et al. (2012) and Phonekeo et al. (2017) with a naive approach to what constitutes an available space adjacent to non-moving agents, assuming no direct knowledge of the agents about where they are relative to the rest of the tower. At each time step, each individual agent randomly chooses an adjacent square to move into, performing a random walk and fulfilling rule (ii). When an agent arrives in a voxel with at least one locked neighbor sharing a corner, edge, or side, it decides to lock, fulfilling rule (iii). Locked agents remain in place, and allow others to move on top of them. Finally, when agents climb on top of locked agents, the locked agent's status changes to covered, fulfilling rule (i). For the sake of simple implementation, we also allow agents to start tower building with a constant probability of spontaneous locking, P sl = 1 20,000 . This model leads to aggregations which grow horizontally rather than upward. An example of a final configurations from one such simulation is shown in Appendix 1 and correspond to the boxed-in panel of Figure S1. This is illustrated in Supplementary Video S1, where each tower growing outward in fractal shapes from a center point. This behavior arises due the higher likelihood of an agent performing a random walk to find other agents near the outer edge of the aggregation.
These results closely resemble a phenomenon known as diffusion-limited aggregation (DLA) (Witten and Sander, 1981). DLA was developed to model the aggregation of metal particles which gather in wispy, fractal shapes, similar to the simulated agent aggregation in Figure S1 for P u = 0, k nl = 1. DLA has also been observed in experimental colloidal aggregation systems, as in Reynaert et al. (2006). Without any rule modifications, DLA is unable to form dense aggregations of agents, because agents on the edge of the aggregation shadow those closer to the center. Hence, we propose several modifications to the behavioral rules which are necessary to mimic the time dynamics of tower shapes experimentally observed by Phonekeo et al. (2017).

Probability of Unlocking
First, we allow locked agents to unlock with a constant probability, as long as they are not covered by other agents. This allows individuals past the first locked neighbor they encounter and move further in toward the center of an aggregation. To model this, we introduce a constant probability of unlocking P u which applies equally to all uncovered locked agents. This rule introduces a distinction between locked agents and covered agents-covered agents cannot unlock.

Neighbor-Influenced Locking Probability
Second, we loosen the requirement that agents must lock upon encountering another locked agent, and instead allow for their probability of locking to increase with an increasing number of locked neighbors. This new rule (ii) replaces the previously discussed rule that individuals lock immediately upon finding a locked neighbor. Instead, an individual has a probability to lock based on the number of locked agents in its neighborhood. We define this probability of neighbor-influenced locking as P nl = k nl N n , with N n representing the number of locked agents in an individual's neighborhood and k nl specifying the increase in probability for each additional neighbor. The neighborhood is defined as a distance of one above below, or horizontally adjacent to the agent's location, highlighted by the blue region in Figure 1B.
The overall probability that a free agent chooses to lock is given by, where P sl is the probability of spontaneously locking. Note that the model allows for up to 26 neighbors, so the value of P sl + k nl N n may be >1. In this case, locking is guaranteed. Therefore, a min function is used to state that when P sl + k nl N n > 1, the locking probability is P l = 1. Additionally, the inverse of the neighbor locking factor, 1 k nl , may be thought of as the number of neighbors required to guarantee locking.
The probability of spontaneous locking provides a baseline probability of locking, to allow for individuals to randomly seed towers. In our simulations, we keep this probability small and set it to P sl = 1 20,000 . The neighbor-influenced locking factor provides the urgency with which an agent locks next to locked neighbors.

Attraction Forces
As we show below in section 3.1 and Appendix 1, the two rule modifications above are unable to reproduce large towerlike structures. Therefore, we extend the random walk model discussed above, and add an attractive "force" representing a behavioral tendency to cluster together. Under this effect, individual agents search their immediate local neighborhood for other agents, and move toward the center of all neighbors. This motion is then perturbed by the randomness associated with a simple random walk model. The resulting velocity is given by, where n i is the number of neighbors in the agent's immediate neighborhood sharing at least one corner, edge, or side with the agent's current position, and c is ratio of the magnitude of attraction relative to the magnitude of randomness. Each agent moves toward the available voxel most closely aligned with the direction of v i . The resulting normalized velocity,v i is defined in (S3) in Appendix 2. Each agent moves into the voxel defined by the surface height at the resulting pixel. Software that simulates agents following these modified rules in MATLAB is provided in Supplementary Code S8.

Set of Modified Behavioral Rules
With the three modifications mentioned above, we modify the first three rules of Phonekeo et al. (2017) and Mlot et al. (2011) into four new local rules: 1. Do not move if agents are on top of you. 2. Upon reaching available space adjacent to non-moving agents, stop and lock with them with probability P l = k nl N n . 3. If stopped and locked, but uncovered by other agents, spontaneously unlock with probability P u . 4. If free to move, move generally toward your neighbors with some random noise as defined by equation (2).

Measurements of Tower Geometry
Each simulation is post-processed to measure the geometry of each tower in order to determine how tower-like the aggregation is. For the final configuration of each simulation, a 2-dimensional height map is constructed by assigning each pixel in the 2D projection of the arena a value equal to its maximum height (Figure 2A). We treat the resulting L × L array of pixel values as an image and apply connected-component analysis (Shapiro, 1996) to identify different towers-a labeled image is generated where any two non-zero pixels that share a corner or edge have the same label. Each agent in the simulation is then given the label corresponding to its horizontal position within the labeled image.
As we are interested in building a single large tower, properties for the tower containing the largest number of agents from each simulation are reported. Three tower properties are considered: number of individuals per tower, maximum tower height, and the ratio of the tower height to its equivalent diameter. Equivalent diameter is defined as the diameter of a circle with area equivalent to the tower's base (Figure 2A).

RESULTS
To gain an intuition for the effects of the modifications to the tower-building rules discussed in section 2.2, simulations were run over a range of locking and unlocking parameters, k nl and P u , across multiple attraction parameters, c, and in section 3.3, across varying densities of agents in the system, ρ. We begin with a parameter sweep across the locking and unlocking parameters and attraction parameter in section 3.1. Then, selecting a pair of locking and unlocking parameters, we systematically vary attraction c to show a rapid phase transition, and investigate the time dynamics of tower properties, both near and far from the phase transition in section 3.2. In section 3.3, we vary the density of agents along with attraction, and observe, in section 3.4, that the center of mass of the towers may continue to move. Finally, we optimize for tower size and height in section 3.5 and identified a set of parameters where a tower formed of nearly all individuals in the simulation.

Tower Geometry
To explore the range of possible tower shapes in the model, we sweep the parameter space of the three rule modifications, including probability of unlocking P u , neighbor-locking factor k nl , and attraction factor c. Resulting tower properties and example final configurations from these simulations are shown in Figure 2. Every data point represents the mean of the largest tower's properties for each of 10 simulations. The left column of the array of tower properties, representing simulations with c = 0, shows that without attraction, towers tend to contain a small number of agents, a small height, and an especially low height-diameter ratio. These simulations with c = 0 represent the first two rule modifications-individual unlocking and neighborinfluenced locking-alone. From the measured tower properties in Figure 2B, we see the effects of the first two rule modifications without attraction. The aggregations with the largest number of agents are found in the simulations with parameters k nl = 1 and P u = 0, representing the case of no rule modifications at all. These aggregations lead to diffusion-limited aggregation as discussed above and shown in Supplementary Video S1. The locking and unlocking rule modifications, therefore, decrease the number of agents in the largest aggregation. They do provide an increase in tower height and the height-diameter ratio. This increase is modest, however, with the tallest average tower height reaching 3.4 agents tall for P u = 0.002, k nl = 1 12 , corresponding to a height-diameter ratio of 0.314. The largest height-diameter ratio occurs for the parameters P u = 0.02, k nl = 1, reaching a value of 0.49, with a corresponding average height of 2.2 and 19.9 agents in the largest tower for each simulation. Finally, the simulations with c = 0 and P u = 0.2 with a small lock factor k nl ≤ 1 8 finish the simulations without forming aggregations. Supplementary Video S2 and the c = 0 configuration snapshot in Figure 2C show the dynamics and final configuration, respectively, of one such simulation which is unable to form aggregations, with parameters P u = 0.2, k nl = 1 26 , c = 0. These tower measurements show that without attraction, all of the tested parameter sets produce aggregations that remain small in number of individuals, do not reach average heights more than 3.4 layers, and remain wide and shallow.
When an attractive force is added, larger aggregations form, as shown by the center column of Figure 2B for an attraction ratio of c = 1. As unlock probability P u increases and lock factor k nl decreases, larger aggregations form, with the largest reaching over 500 individuals. However, these largest aggregations have the smallest height-diameter ratios of this set, showing that these large aggregations are particularly wide, as is visible in the c = 1 example in Figure 2C and Supplementary Video S3. Increasing the attraction ratio to c = 2 finally reveals a more typical towerlike shape, with taller aggregations, even reaching a height of 11 agents. Interestingly, these taller towers contain fewer agents than the c = 1 case. The reason for this is clear in the snapshots shown in Figure 2C and Supplementary Video S4 -stronger attraction yields more densely-packed towers with larger height-diameter ratios-the towers are so dense that multiple, smaller towers form instead of most individuals aggregating into a single tower.

Phase Transition and Time Dynamics of Tower-Building
The example configurations shown in Figure 2C represent the same set of locking and unlocking parameters, P u = 0.2, k nl = 1 26 across c = {0, 1, 2}. These locking and unlocking parameters give the largest towers for both c = 1 and c = 2, but no aggregations at all for c = 0. To investigate the effects of the attraction ratio c further, we selected a fixed pair of locking and unlocking parameters and explored both the height and number of agents in the largest tower in the system for a densely sampled range of the attraction parameter c. The results of these simulations are shown in Figures 3A,C. The presence of a phase transition occurs between c = 0.92 to c = 1.06, where the number of agents in the largest tower climbs from close to 0 to over 700 agents. The results show that as c increases beyond that critical value, the number of individuals in the largest tower decreases (Figure 3A) while the height of the largest tower increases (Figure 3C).
In Figure 3B, we show the time dynamics of the number of agents for tower for two cases, close to the phase transition and further from it. To illustrate tower growth further from the phase transition, Figure 3B shows the time histories of all 10 simulations for P u = 0.2, k nl = 1 26 , c = 2.0 in green and the mean of all simulations in black. One of these simulations is shown in Supplementary Video S4. The tower formation process in this model demonstrates two time scales: the time scale of initial nucleation, and the time scale of growth. Nucleation generally occurs within the first 5,000 time steps, the first 1% of each simulation. After nucleation, towers often continue to grow slowly through the rest of the simulation. Occasionally, two towers will merge into one, which manifests as a sharp jump in the time histories of Figure 3B. Some of these tower collisions last through the rest of the simulation, while others briefly merge and then separate again, which shows up as a sharp peak in the time history of tower size. The fast nucleation followed by slow growth seen for c = 2.0 is typical for most simulations in the present work.
However, there are some examples, particularly within the phase transition regime, for which a critical slowing down occurs. Trajectories close to the phase transition are shown in Figure 3D. Two trajectories are shown for each of c = {0.96, 1.0, 1.3} with P u = 0.2, k nl = 1 26 . The critical slowing down is particularly evident for the c = 0.96 trajectories, where agents aggregate into a tower after 250,000 time steps while the other never transitions out of the disordered state. The c = 1.0 trajectories also show variation in nucleation time, although in this case, all simulations have transitioned to their aggregated state, in which the largest tower contains at least 100 agents. The variation in tower size is highest for these examples, varying in size by 200 or more individuals. There are also cases where the towers continue to grow in size, even after 500,000 time steps, which can also be seen in the case of c = 1.3.
As discussed in section 2, the simulation time of 500,000 was chosen because nearly all simulations have reached a steady state. The cases highlighted in Figure 3D represent the exceptions, and there is no guarantee that these simulations will ever converge. The figure shows that c = {0.96, 0.98} are the only cases that give a mixture of aggregated and non-aggregated results.

The Effect of Density
The parameters varied up to this point in the model represent entirely behavioral parameters, that is, those associated with the decision-making of individuals. While these parameters are testable within multi-agent robotic examples, they do not represent a variable that can be systematically changed in experiments with live fire ants, or robots, in order to test the predictions of the model. To develop a set of testable predictions, we turn to explore the parameter of density of agents, ρ.
In our model, density is varied by changing the number of individuals in a fixed arena size. The computational complexity of the model is O(N 2 ), so practical limits of computational time place an upper bound on density we explore here. In a 100 × 100 × ∞ arena, our test set is N = {200, 500, 750, 1000, 1500, 2000} which corresponds to the densities ρ = {0.02, 0.05, 0.75, 0.1, 0.15, 0.2}. We will use unlocking and locking parameters of P u = 0.2, k nl = 1 26 for consistency with section 3.2.
The results of these simulations are presented in Figure 4, showing several key differences and similarities across densities. As density increases, less attraction is required for tower formation. The data points highlighted by circles in Figure 4B show the critical attraction ratio c * , which represents the minimum value of c for which the largest aggregation is at least 100 individuals, representing the onset of the phase transition. This result also implies that there exists a critical density across a range of attraction factors, below which no tower formation occurs. Another key result is that the largest towers, in terms of number of agents per tower, occur shortly after the transition from no towers at all. agents, height, and height-diameter ratio over a range of density ρ and attraction factor c. Note that the vertical axis is not linear. Each data point represents the average properties of the largest tower from 10 simulations after 500,000 time steps. The circles on the density plot represent the onset of phase transition. The circles indicate the minimum attraction coefficient c for each density at which the largest tower contains at least 100 agents.
Beyond this point, tower height and ratio increase while number of agents decreases. Finally, tower shape remains close to constant across densities, particularly in heightdiameter ratio.

Moving Towers
The introduction of unlocking probability effectively adds noise to the system (equivalent to higher temperatures in thermodynamic systems), which allows towers to move. Agents locking on one side of the tower while others unlock on the other side can lead to tower motion. Traces of the center of area for each tower from two example simulations may be seen in Figure 5B. To quantify this phenomenon, we consider the motion of towers as a Brownian random walk and investigate the diffusion coefficient of each tower. The diffusion coefficient (D) for a Brownian random walk follows the relationship, for each trajectory of length T. Therefore, we measure the mean square displacement (MSD) of each tower in each simulation over a variety of times, t = {0, 250, 500, ..., 12, 500}, and perform a linear fit for each tower trajectory. The average slope of these lines is then twice the diffusion coefficient ( Figure 5A). These results show that the maximum diffusion occurs in the highest density regime, and for the lowest attraction parameters that generate aggregations, particularly for ρ = 0.2 at c = 0.75 and c = 1.0. These towers have lower height-diameter ratios, as seen in Figure 4B, which leads to a larger proportion of agents on the surface of the tower, and therefore a higher probability that individuals on the surface will be unlocking. The towers at c = 0.75 have a smaller number of agents than those of c = 1.0, which leads to an even higher proportion of individuals on the surface. This is illustrated in Figure 5C, showing the timeevolution of tower configurations for two example simulations (see also Supplementary Videos S5, S6).

Tower Optimization
One question that still remains is, what parameter values are optimal for tower building? To answer this, we need to think about what may constitute optimal. It may be that the optimal tower reaches as high as possible, which would, in practice, allow as many agents as possible to attach to a support structure. Or, for robotics applications, this would allow the tower to reach higher heights. On the other hand, it may be best to include as many individuals as possible in the tower, and the optimal tower would be the one that includes every single agent in the tower. As observed in Phonekeo et al. (2017), fire ants built towers that equally distribute load among the individuals. Therefore, an optimal tower from their perspective may be one that optimizes for load distribution. In this section, we use a genetic algorithm to explore optimal tower building considering each of these optimization targets.
To search for an optimal tower, we employ the Covariance Matrix Adaptation-Evolutionary Strategy (CMA-ES) algorithm developed by Hansen et al. (2003). This algorithm randomly generates parameter sets within the search space and evaluates a cost function for each parameter set. From the results of this function evaluation, it updates the covariance matrix to expand in the direction of the most optimal value. Using the updated covariance matrix, the algorithm generates new parameter sets and repeats the process until convergence, generally defined as finding a parameter set with a cost function below some threshold. We applied the CMA-ES algorithm to the tower-building model introduced above, using the average final properties of three trials for each parameter set across 10 parameter sets per iteration. For the optimization, we choose a cost function defining the optimal tower as the largest tower, both in terms of tower height and number of individuals within the tower. Therefore, the cost function is given by, where N tower and h tower represent the number of individuals and height of the largest tower, N max is the number of individuals in the simulation, and h max is a prescribed maximum height. The height term is included to ensure that the results are effectively tower-like, preventing the optimal tower from simply achieving a large, wide aggregation. From the results of the attraction sweep in Figure 3A, we observe that h max = 14 is an approximate upper bound on tower height, so it is therefore chosen as h max for the purpose of this optimization. Note that a tower height of h tower ≥ h max results in a zero second term, and the simulation therefore allows for a taller tower. For the purposes of optimization, we reduce the simulation time to 50,000 time steps. This serves the practical role of making iterated simulation possible, but also places an effective minimization of convergence time. Therefore, we are optimizing for a tower that maximizes both height and number of agents quickly (within 50,000 time steps). Figure 6 shows the progression of the minimum cost at each iteration of the CMA-ES algorithm along with snapshots of intermediate results to show the algorithm's progress. The optimal tower occurs for the parameters, P u = 0.938, k nl = 0.029, and c = 2.56, which led to a tower of 993 agents reaching 16 agents tall after 50,000 time steps. The final cost function, averaged over three trials, was f = 0.01.
A video of one simulation with these parameters is shown in Supplementary Video S7.
The CMA-ES optimization code of Hansen et al. (2003) applied to the present model may allow future research and consideration of other conditions of optimal tower-building. For example, when designing a robotic system where each individual robot has a maximum load capability, it may be necessary to calculate the maximum load experienced by an individual in the tower and add that term to cost function.

DISCUSSION
In this work, we have extended a previously proposed set of local rules to replicate the tower-building behavior of red imported fire ants, Solenopsis invicta. This model and its insights will allow for the design of control strategies for tower-building swarm robotics and greater insight into the collective behavior of social insects. The results presented above show that individuals moving under the influence of local attraction are able to form large towers. We find that an attractive force is necessary for significant tower-building and show the impacts of this attractive force over a range of locking and unlocking parameters as well as a range of densities. We find that the system contains a sudden phase transition as the attraction parameter is varied, and that this phase transition is density-dependent. Finally, the largest towers, in both height and number of individuals, occur with a combination of very strong attraction and highly probable unlocking.
On the other hand, without attraction, no towers form, as shown in the c = 0 case of Figure 2 and Supplementary Video S2 and discussed further in the Appendix 1 and Figure S1. The effective force of attraction may also be thought of as a desire of the ants to climb, because the tallest available square to move toward will also have the most neighbors.
Near the phase transition, a critical slowing down occurs, and there are parameter sets that do not result in tower formation within a simulation time of 500,000 time steps. This critical slowing down is reminiscent of other examples of systems with phase transitions, such as the spin-glass model, the Ising model, and molecular dynamics models (Dasgupta et al., 1979;Hu, 2013). Further from the phase transition (c ≫ c * ), towers form rapidly, but the possibility exists for these towers to encounter one another and merge into larger towers. The number of individuals in the tower and tower motion are largest just after the phase transition, but the largest height occurs with stronger attraction. Phase transitions have previously been observed experimentally and computationally in other ant and insect systems, such as in Pharaoh's ant foraging (Beekman et al., 2001) and in marching desert locusts (Buhl et al., 2006).
Our results also illustrate the exploration-exploitation tradeoff, which balances attraction forces with random movement and unlocking events. Following this trade-off, stronger attraction may lead to higher towers with fewer individuals, as the attraction rapidly draws individuals from the edge of the aggregation toward the center of the tower, and therefore upward. This balance of unlock probability and attraction is found through the combined optimization of number of individuals and tower height, which discovered that with an unlock probability of P u = 0.938 and an attraction of c = 2.56, it is possible to include nearly all of the individuals in a simulation, with a tower reaching a height of 16 layers. This large unlock probability of the largest towers in our simulations connects with the observation from Phonekeo et al. (2017) that, in the experimental system, the fire ants are constantly rebuilding their tower and circulating ants throughout the tower. The work of Phonekeo et al. (2017) showed that fire ants build towers of constant load, and future optimization work could incorporate the load experienced by each individual to achieve towers that prioritize stability. More refined ant models may also incorporate the mechanics and viscoelastic properties of fire ant aggregations (Tennenbaum et al., 2016), which are observed to change depending on the number of active ants, such as the free ants included in the present model (Tennenbaum and Fernandez-Nieves, 2017).
The results of the parameter sweep in density values showed both similarities and differences across densities. In general, for a fixed attraction ratio c, the tower height-diameter ratio remains fairly constant, even as the numbers of agents per tower and tower height vary. The biggest difference across densities is the change in critical attraction parameter, c * . These observations lead to testable hypotheses for animal experiments. Below a certain density threshold, tower formation should cease, due to the move past the critical attraction. Additionally, the heightdiameter ratio should remain constant across a large range of densities. Finally, we have shown that the towers built in our simulations move over time, with a diffusion coefficient that is dependent on both attraction and density, and should be taken into account when considering practical application to robotics.
This work also lays the groundwork for future robotic studies, where robots are able to built towers out of themselves in a manner similar to, for example, the M-blocks of Romanishin et al. (2015) or the Roombots of Spröwitz et al. (2014), which have also been proposed for bridge-building applications (Nguyen-Duc et al., 2019). The tower is a ubiquitous structure in building, and designing rigorous control strategies for towerbuilding represents a fundamental starting point toward fully autonomous, locally-sensed swarm building applications. In practice, a tower of robots could be useful in the case of, for example, seeing over obstacles, providing scaffolding for climbing, or clearly marking a location of interest. Robotic towerbuilders would need to be have the following capabilities: sense neighbors, climb onto and off of one another, and support appropriate loads. At the moment, there is no robot with all of these capabilities, and we believe that this would be a promising avenue for future robotics research. The control strategies introduced in the present study could also be further modified to more closely replicate experimentally-observed fire ant behavior, developing a control strategy for interacting with a support structure.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are available upon request.

AUTHOR CONTRIBUTIONS
NM, JC, TS, and JL developed the initial version of the present model as a project in a course taught by OP. GN refined the model, conducted the simulations, and wrote and edited the manuscript. OP supervised the research and edited the manuscript.

FUNDING
The study was conducted with institutional funding.

ACKNOWLEDGMENTS
We thank Prof. David L. Hu, and members of the Peleg lab for insightful discussions, and the BioFrontiers Institute for the utmost support.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frobt. 2020.00025/full#supplementary-material Supplementary Video S1 | The diffusion-limited aggregation case of the model, with P u = 0, k nl = 1, c = 0. The video is shown at a speed of 120 time steps per second.
Supplementary Video S2 | Simulation in which no aggregations form, with P u = 0.2, k nl = 1 26 , c = 0. The video is shown at a speed of 10,000 time steps per second.
Supplementary Video S3 | Simulation in which large, wide aggregations form, with P u = 0.2, k nl = 1 26 , c = 1. The video is shown at a speed of 10,000 time steps per second.