# Modeling Multi-Agent Self-Organization through the Lens of Higher Order Attractor Dynamics

- Department of Psychology, University of Utah, Salt Lake City, USA

Social interaction occurs across many time scales and varying numbers of agents; from one-on-one to large-scale coordination in organizations, crowds, cities, and colonies. These contexts, are characterized by emergent self-organization that implies higher order coordinated patterns occurring over time that are not due to the actions of any particular agents, but rather due to the collective ordering that occurs from the interactions of the agents. Extant research to understand these social coordination dynamics (SCD) has primarily examined dyadic contexts performing rhythmic tasks. To advance this area of study, we elaborate on attractor dynamics, our ability to depict them visually, and quantitatively model them. Primarily, we combine difference/differential equation modeling with mixture modeling as a way to infer the underlying topological features of the data, which can be described in terms of attractor dynamic patterns. The advantage of this approach is that we are able to quantify the self-organized dynamics that agents exhibit, link these dynamics back to activity from individual agents, and relate it to other variables central to understanding the coordinative functionality of a system's behavior. We present four examples that differ in the number of variables used to depict the attractor dynamics (1, 2, and 6) and range from simulated to non-simulated data sources. We demonstrate that this is a flexible method that advances scientific study of SCD in a variety of multi-agent systems.

## Introduction

For many animals and humans, social interaction is pervasive in daily life. Social interaction occurs across many time scales and varying numbers of agents; from one-on-one to large-scale coordination in organizations, crowds, cities, and colonies. Since social interactions occur at different scales, and in ways that change dynamically over time, they can be quite a complex phenomenon to study without appropriate guiding theoretical and methodological frameworks.

In dynamical systems theory, complexity arguably occurs due to the emergent, self-organizational nature of the system. Emergent, self-organization here implies that there are higher order macroscopic patterns occurring over time that are not necessarily due to the actions of any particular controlling agents or components, but rather due to the collective ordering that occurs from the individual interactions of the agents or components of the system (Halley and Winkler, 2008). Taken together these diffuse interactions contribute to more macro-scale phenomena that are observed over time. Common examples of this type of emergent, self-organization of social behavior occurs in flocking birds and in schools of fish that appear to move in a highly coordinated fashion (e.g., Couzin and Krause, 2003).

Because of the multitude of agents (or system components) that give rise to emergent patterns, it is often difficult to determine how one should depict the resultant system. In line with approaches to *social coordination dynamics* (SCD), we aim to uncover the dynamic processes that underlie the ways in which agents are able to organize their behavior and change together in time (Oullier and Kelso, 2009). This emergence is a form of c*oordination* that specifically implies the occurrence of a functional ordering of components that interact across spatial and temporal dimensions, often with multi-directional relationships (Kelso, 2009; Butner et al., 2014a). We aim to model this emergent, multi-agent coordination through attractor dynamics depictions (which we discuss in detail in the next section).

From the SCD perspective, two or more agents are able to coordinate their behavior based on some form of mutual information exchange. This information exchange generates coordinative structures with higher order patterns not easily identifiable from the lower order interactions. The resultant higher order patterns are then depicted through attractor dynamics in which the patterns are attributed with stability properties implied by an underlying topology (Kelso, 1995). SCD is thus consistent with notions of weak emergence (Bedau and Humphreys, 2008), but the resultant patterns have then been modeled using attractor dynamics descriptions depicting patterning over time that have stable properties under perturbations. One major distinction between SCD and examples of weak emergence (usually through agent-based or cellular automata models) is in the scale of the social systems involved.

Agent-based models are usually quite large-scale social systems, while SCD has often focused on a dyadic scale of analysis. SCD has excelled in generating models of intentional and spontaneous dyadic interpersonal rhythmic behavior such as finger or limb oscillations (e.g., Haken et al., 1985; Schmidt et al., 1990; Oullier et al., 2008), swinging pendula (Schmidt and O'Brien, 1997), and rocking in chairs (Richardson et al., 2007). Some recent work has provided ways to assess social interactions in larger scales such as coordination of groups bigger than dyads (e.g., Richardson et al., 2012; Duarte et al., 2013). One challenge is to generate models of emergent, multi-agent coordination in social systems where the agents may not behave rhythmically *per se*, but are following some organizing rules or structures that give rise to coordinated behavior serving a functional purpose.

In the current paper we build on SCD approaches, by modeling the results from large-scale agent-based systems as a function of attractor dynamics. Our chosen technique utilizes mixture modeling in conjunction with topological equations to represent attractor dynamics. This approach is particularly attractive in that topological representations of phase space can yield more qualitative information in comparison to other time series approaches (Strogatz, 2014), generating a more complete picture of the underlying system dynamics. Specifically, we examine a series of agent-based examples, and model each set of time series as a function of their changes through time. We show how sets of linear equations can depict the higher order emergent patterns in ways consistent with attractor dynamics. The advantage of this approach is that we are able to quantify the self-organized dynamics that agents exhibit, link these dynamics back to activity from individual agents, and relate it to other variables central to understanding the coordinative functionality of a system's behavior. Our goal is to exemplify the strategy. In all, we present four examples that differ in the number of variables used to depict the attractor dynamics (i.e., the dimensionality of the systems) and range from simulated to non-simulated data sources.

### Attractor Dynamics

In dynamical systems theory, the concern is often placed on what states a system is drawn toward, or away from, as it changes over time (e.g., Richardson et al., 2014). This epitomizes the notion of attractor dynamics. *Attractor dynamics* are merely a mathematical way of expressing repetitious behavior in the face of constant disruptions to those repetitions. The constant disruptions are inherently part of the system in that open systems are dissipative and function far from equilibrium to maintain patterns (Prigogene and Stangers, 1984). By only examining a portion of a system, as is common in empirical research, the unexamined variables are treated as constant disruptions or perturbations to those patterns. These repetitious behaviors describe the most probable system states and their ability to remain in these states (while facing perturbations) conveys the inherent stability of those states. These attractor dynamics can then be modeled using differential/difference equations, allowing for exploration of the inferred dynamics and theorizing the manifolds in which the system functions (Differential equations are based on idealized models for when change in time approaches zero while difference equations estimate models using the observed discrete differences; Butner et al., 2015).

Assuming a system exhibits stability, the emergence of a limited set of patterns, which can be described in terms of topological features, are plausible. These topological features can be described using map analogies, because there is a strong tie between topology and maps. In fact, differential topology *is* the math behind maps. Traditionally, topographical maps convey elevation of a landscape. But, the notion of topologies can also be applied as a graphical representation of how data are changing over time.

To ease the interpretation of differential topology, we will temporarily link movement on maps to different topological features. This interpretation is directly relevant to several of the examples (although more general definitions are extant; Butner et al., 2015). An *Attractor* is when the agents are drawn toward a particular coordinate over time or a particular directional heading. This is akin to a topographical valley. A *Repeller* is a coordinate that agents move away from. These would be reflected topographically as a mountain peak. A *Saddle* occurs when agents are attracted in one dimension and repelled in another. It is analogous to a topographical ridgeline because it can separate different patterns such as two attractors (Abraham and Shaw, 1983; Butner et al., 2015). A *Cycle* corresponds to a push/pull of two dimensions on one another. Combined versions of patterns described can also be observed such as spiral attractors where there are circling movements for how agents converge toward an attractor. Saddles and cycles require at least two dimensions and thus will only be possible in the later examples (not merely with heading as it is a one dimensional example).

To continue with the link to maps, we will begin with agent-based models that function spatially. As a simplification, we can reduce their behavior to movement along an X and Y axis or merely the directional heading of agents (when we only require a single dimension to depict the system). We can then model the simultaneous change of these variables over time. In this way, we capture the movement of many agents and can characterize them with attractor descriptions. With this information it is possible to examine and identify patterns of change for the overall system using the particular topological features defined above to describe how multiple agents are moving over time (Butner et al., 2015). It is in these terms that we gain an understanding of the emergent, coordination of many agents.

As a beginning example, consider a Flocking agent based model (Wilensky, 1998) in NetLogo v5.2.1 (Wilensky, 1999) designed to emulate the self-organized behavior of how flocks of birds might come to match one another's movements creating complex group behavior. Agents start with a random heading and constant velocity in a wrapped environment (makes a torus). The heading for each agent is determined by three rules: (1) *alignment* states that each agent tends to turn to be moving in the same direction as nearby agents; (2) *separation* states that each agent will turn to avoid an agent when it gets too close; and (3) *cohesion* states that agents tend to move toward other agents. As the agents “fly” through the two dimensional environment they update their headings over time. Figure 1 shows time series of the headings for all (300) agents simultaneously over one thousand iterations. It is clear that early in the simulation the full range of headings are observed yet, in later times the range of headings become more restricted and shared by the agents. This is an example of the emergent coordination that occurs within the Flocking model.

**Figure 1. Time series of the headings for 300 flocking agents**. Note that the flock moves toward a very restricted heading.

To depict these results topologically requires identifying the underlying map in which the agents are interacting. An attractor, in this case, would be the heading(s) in which the agents move toward and the stability would be the resistance exhibited in the system when an agent begins to diverge from this attractive heading and pulled back thusly. The map is not one of actual hills and valleys, but instead the resultant decrease in heading directions, from the emergent self-organization between agents. Thus, the trajectories for agents imply an underlying pattern that we can infer. One assumption of dynamical systems theory is that there is one—or perhaps multiple—underlying patterns emergent from the interactions of individual agents over time. Interactions result in a consistent pattern, that the system flexibly returns to when interactions or outside forces briefly move the system away from its primary pattern. This notion of consistency in the face of perturbations is stability. While it is easy to observe the convergence of heading amongst the agents in Figure 1, little information can be drawn in regards to the number of underlying patterns and their inherent stabilities.

The flocking example is a useful one in that the implied map is not a map of X and Y coordinates, but one of heading—it is a one dimensional map. One dimensional maps are not very interesting to draw; they are a line showing where the data converges over time. That is, attractor dynamics are time implicit models rather than time explicit ones and thus, are akin to collapsing the X axis in Figure 1, while adding in notions of *where each agent goes next* to determine the map.

Dynamical systems theory has long provided the theoretical framework and terminology for describing multi-agent self-organized patterning. Returning to Figure 1, an apt depiction of the Flocking simulation is one that begins with many attractors that cease to exist over time, which produce a limited set of stable attractors. This qualitative description captures the evolving process, without any of the quantitative dynamics. By quantifying them through topological equation representations, we can further differentiate aspects of the system and specify the strength of the attractors. We therefore next cover the step of quantifying the dynamic.

#### A Vector Based Approach

There are several ways to estimate differential topological equations. In all cases, we must first express the data in terms of data vectors rather than values. For the heading data illustrated in Figure 1, the data is structured such that two or more points in time are used to define a data vector, known as a time delay or Toeplitz data structure (e.g., Boker and Laurenceau, 2006). The data is structured so that a value at time t and a value at time *t*+1 are two variables within the model. Further, our models are all estimated in structural equation modeling wherein change was built into the models themselves as latent variables (McArdle, 2009). One can also estimate change directly through a discrete difference or various methods for estimating derivatives and thus while we use structural equations to build our models, this is far from a necessity (Boker et al., 2010).

Different attractor dynamics are then captured through expressions of change predicted by value. For example, Equation (1) expresses the potential dynamics for the headings of the various agents.

Current heading of a given agent at time *t* is *x, x*-dot is an estimate of its derivative with respect to time, *b*_{0} is the intercept, *b*_{1}, the slope with respect to *x* and *e*_{t} is error. For clarity, this equation is written in regression form where velocity in heading at each point in time is treated as the criterion and position (current heading) is the predictor. When the slope in Equation 1 is negative, we observe an attractor where the time series are attracted toward a value of −*b*_{0}/*b*_{1} known as the set point (Butner et al., 2015). A repeller occurs when the slope is positive instead of negative. The strength of attraction/repulsion is defined by the steepness of the slope relative to zero.

Equation (1) is limited in that it can only capture a single topological feature (Butner et al., 2015). While the system may converge to a single heading, this convergence is developed over time. Consistent with the qualitative description of the flocking model, we should observe several patterns that cease to exist as time continues. This results in a much more limited set of dynamic patterns that occur at later times. We therefore expand our approach to allow for multiple sets of Equation (1). We did this through an analytic technique known as mixture modeling.

#### Mixture Modeling Methods

Mixture modeling is a taxonomic approach that can be combined with structural equation modeling (Enders, 2006) as an alternative way to capture interactions (Jung and Wickrama, 2008). Non-linear dynamical systems allow for multiple topological patterns by building non-linear transformations, such as the interaction and therefore mixture modeling can be used as a way to capture the different topological features by slicing up the overarching state space under an assumption that each dynamic is locally linear. One description of mixture modeling is as a multiple group analysis (stacked model), where assignment to group is unknown (Muthen, 2001). Multiple group models allow for different parameters across groups. We can extract different equation sets by allowing key parameters to differ across these groups while equating others. Specifically, we allowed the slope coefficients characterizing how position predicted each velocity, the intercepts for the velocity factors, the means for the position factors, the residual variances for the velocity factors, and the variances for the position factors to vary across sets of equations [see Appendix A (Supplementary Material) for an example in Mplus (Muthén and Muthén, 1998–2012)].

As previously described, the sign of the slope coefficients capture the type (e.g., attractor, repeller, limit cycle) and strength of attraction for the dynamic implied by the equations (see also Butner et al., 2015). In addition, the velocity intercepts help determine the set point, or relative position to which the dynamics can be described (e.g., the location of the attractor). Following logic laid out under notions of centering and simple slopes analysis (Cohen et al., 2003), the means and variances for the position factors help depict common trajectories implied by the pattern and thus help identify the basin of attraction. By allowing for variation in these parameters across latent classes, we can infer a number of varying topological features, as opposed to a single feature.

Mixture modeling can be used as a confirmatory or exploratory method. In either case, there must be established criteria for fit. The current preferred methods are through forms of the Bayesian Information Criterion (BIC) or through forms of model testing such as log likelihood or chi-square comparisons to see if the current number of extracted groups improves description of the data beyond the previous number of groups. Specifically, the BIC and sample size adjusted BIC tend to minimize when the proper number of mixture groups has been extracted (Sclove, 1987; Nagin and Tremblay, 1999; Nylund et al., 2007) and both have been used under different circumstances usually relating to the sample size (sample size adjusted BIC is preferred when *n* < Bauer and Curran, 2003; Lubke and Neale, 2006; Enders and Tofighi, 2007).

Model identification can also be informed by various likelihood ratio tests (LRT), which are used to test relative model fit by testing the null hypothesis that competing models demonstrate comparable fit (Vuong, 1989). Within latent variable models such as the present one, the Vuong-Lo-Mendell-Rubin test (Lo et al., 2001) is an accepted methodology for testing the equivalence of two associated probability density functions (Henson et al., 2007). Simulation studies have indicated that the VLMR test favors selection of more components when used with small samples, resulting in increased Type I error rates; this suggests the need for an adjusted test (aVLMR) with samples less than 300 (Lo et al., 2001). For our purposes, we chose to rely on the BIC.

Note that our data had an inherent dependency—the nesting of multiple measures through time within each agent. Ignoring a data dependency is known to produce biased standard errors with large alpha inflation as the common result (Cohen et al., 2003). However, current mixture modeling practices that incorporate methods for accounting for the dependency preclude any descriptions of predictors. In this case, that would result in the loss of the means and variances for the position factors that detail key information about the basins of attraction. We therefore chose to temporarily ignore the dependency, recognizing that the standard errors for each coefficient may be biased toward Type 1 errors.

To better understand the extracted equation groups, we saved out the *posterior probabilities* for each data vector. This is the probability that each instance in time for a given agent belonged to one of the classes characterized by a particular equation set where the set of posteriors for a given vector sum to one. It is the equivalent of factor scores if mixture groups as likened to a categorical latent variable. The value of the posterior probabilities is that they allow us to specifically link each agent to the various attractor dynamics at each point in time. Through the combination of the description of each attractor dynamic and the posterior probabilities linking the agents to the topologies, we are able to traverse between the observed vectors from the agents to the underlying topology.

### One Dimensional Systems

What follows is an illustration of the analytic strategy for the flocking example using the headings from all agents. Fit indices of the 300 flocking agents over 1,000 iterations resulted in sixteen unique attractors (as indicated by the BIC at its lowest value). Table 1 contains the estimated parameters for each of the sixteen equations. All sixteen patterns are attractors as indicated by the negative slopes. They vary in their stability, indicated by the range of slopes. The headings to which each pattern indicates a point of attraction is identified by converting the intercepts and slopes into the set point (−*b*_{0}/*b*_{1}). In essence, the flock example is characterized by a total of sixteen unique attractors.

**Table 1. Unstandardized coefficients from the sixteen attractor solution for the Flocking model of headings**.

We can link the attractors back to the individual agents through the posterior probabilities. For purposes of relating to the initial assessment of the many unique patterns dying off, we chose to illustrate the average posterior probabilities (the average likelihood a given agent is depicted by a given attractor) as a function of time. Figure 2 shows the average posterior probabilities for each attractor dynamic. The legend shows the heading attracted to (set point) and level of attraction (slope) as a function of time. Consistent with Figure 1 (and expectations), initially there were many attractors, but somewhere around iteration 300, two specific attractors started to dominate (dotted lines in Figure 2).

**Figure 2. Time series of the average posterior probabilities for each attractor dynamic pattern**. All patterns were attractors as indicated by their negative slopes, but differed in terms of their set points (SP; the heading in which the agents were attracted toward) and the attractor strengths (the steepness of their slopes). At around iteration 300, the system shows a phase transition wherein two patterns with the same heading begin to dominate.

Notice that they share the same heading of 273 degrees, but with slightly different degrees of attraction. Recall the three rules that constitute the changes in heading over time: alignment, separation, and cohesion. Alignment and cohesion drive the agents toward a single heading, but separation instead evokes divergence when agents become too close (and specifically overrides the other two rules). What distinguishes the patterns is not the heading they are drawn toward, but in the divergences themselves due to separation that produces a weaker attractor. Note that agents can be switching between the two attractors over time moving to the slightly weaker attractor, as they need to avoid collisions.

We gain additional information from the quantitative attractor dynamic description as illustrated in Figure 2 when compared to Figure 1. Each data vector is now depicted not only in terms of its vector, but also the likely attractor in which it is drawn (through the posterior probabilities). Further, the description is now in terms of the underlying system forces that depict the type of pattern (all attractors since all the slopes were negative), the location to which the patterns are relative (the set points), and their stability under perturbations (the deviation of the slopes from zero). However, thinking topologically becomes even more beneficial as we move toward systems with more dimensions.

### Two Dimensional Systems

Modeling a two dimensional system can be captured through two simultaneous equations.

These equations represent two variables measured simultaneously (x and y) at time t, x-dot and y-dot are their estimated derivatives at time *t*, *b*_{0} and *b*_{3} are intercepts, *b*_{1} and *b*_{5} are each variable predicting its own derivative, *b*_{2} and *b*_{4} are crossover or coupling relationships, and *e*_{xt} and *e*_{yt} are errors in equation. Notice that Equation (2) is identical to Equation (1) with the addition of the other changing variable also predicting velocity in x (or x predicting velocity in y). By having both variables changing simultaneously, we generate a two dimensional depiction. The emergent dynamic (attractor, repeller, etc.) is a function of all the *b* coefficients in the equations (Gottman et al., 2002). Common interpretation is that the *own effects* (i.e., x predicting change in x and y predicting change in y) depict the stability properties of the dynamic pattern (attractor, repeller, or saddle) such that negative coefficients are indicative of attractive behavior and positive coefficients are indicative of repulsive behavior in the respective dimensions. The crossover relationships (also known as *coupling effects*) are commonly interpreted to represent the push-pull of variables that constitute cycles and swirling qualities graphically. The set point is a function of both equations. And as noted earlier, two-dimensional systems can include saddles and cycles, which are topological features that are not possible in one-dimensional systems.

While many cases can be interpreted as described in the previous paragraph, some cases do not always conform to the conventional interpretations (and we include some examples of this below). A common violation relates to the notion of collinearity. If all variation in both x and y perfectly map onto one another, then x and y are essentially a single dimension. Under this circumstance the coefficients can be misrepresentative of the dynamic pattern. In our spatial movement circumstance, agents will sometimes capitalize on diagonal movement as a primary, singular dimension. Assessment of the eigenvalues and eigenvectors of the coefficients (treated as a Jacobian matrix of partial derivatives for estimating local Lyapunov exponents; Arabanel et al., 1992) is a method for verifying and determining whether to follow the classic interpretation or whether the interpretation should be modified.

#### The Ants Model

Consider the Ants model (Wilensky, 1997) in NetLogo v5.2.1 (Wilensky, 1999). This agent-based model was designed to simulate ant colony foraging behavior. The simulation consists of 125 ants each with the same instructions, starting at a nest in the center of a two-dimensional space. Ants are released one at a time from the nest, moving at a constant velocity. Three food sources are placed within the two-dimensional space each with a finite quantity of food supply. The ants search the environment for food (following a random direction algorithm) and upon locating and collecting food, return it to the nest. The primary mechanism for the emergent foraging behavior involves the ants releasing digital pheromones while carrying food and that the ants are attracted to this pheromone. This is much like how stigmergy, a form of environmental modification by individual social animals that affords collective coordination, is proposed to work in live ant populations (Theraulaz and Bonabeau, 1999). The nest also releases a pheromone signal so that the ants can find the nest. The simulation allows for the manipulation of the evaporation and diffusion rates of the pheromones, which we left at default settings. Figure 3 shows the standard placement of food sources in the environment in relation to the nest at the center.

From visual inspection, several emergent colony behaviors can be observed. Ants will search the environment until a critical threshold of ants find a given food source. At this point the ants will form a trail between the food source and the nest. There are sometimes congestion-like behaviors that occur in the middle of the trail or near the nest as more ants converge toward the strongest pheromone locales. Once the food source is used up, the ants once again spread out into a search pattern until a new food source is found. In this case, we will depict the attractor dynamics of the ant movement in two dimensions as a way to characterize the different ant behavioral patterns.

We extracted the horizontal (x) and vertical (y) coordinate position of every ant from the beginning of the simulation until the last food pile was fully exhausted, totaling 1,080 iterations. Figure 4 is a kernel density plot of the ant positions, collapsed across all ants and all iterations. This shows the regions where ants spent most of their time and can be thought of as the probability density function of the data (under the assumption of two dimensions)—a graphical illustration of the integral of the dynamics. The density plot is read in the same fashion as a topographical map, where the lines illustrate more density. Note that the greatest density is at the nest (0,0). This was likely a function of all the ants starting at the nest, including the dispersion algorithm of only a single ant leaving the nest per iteration. It is also a function of all the ants returning to the nest to deliver food. Each branch of the density plot corresponds to one of the food sources, consistent with a trail between the nest and the food source. The densest part for each of the branches was, however, closer to the nest than the food source.

**Figure 4. Kernel density plot of where the ants spent most of their time during the simulation**. Note that the highest densities correspond to the three food source locations and the nest.

Figure 5 contains trails of three exemplar ants as vector plots in time to help illustrate the link between individual agents and the model estimated from all agents. Figure 5A shows the trail of an ant that helped collect food from all of the food piles. However, it also shows searching behavior in some of the areas of the world where food did not reside. Figure 5B illustrates an ant that only helped collect food from two piles and also participated in searching behavior in empty quadrants of the world. Figure 5C shows an ant that participated minimally in food collection instead spending more time searching. As a whole, these illustrate that the emergent behavior is not from any one ant. Instead, it is through their interactions with one another (through pheromones) and the environment (food resources relative to the nest) that their behavior becomes emergently coordinated.

**Figure 5. (A-C)** Three example ant trails that illustrate how the ant behavior is shared across all the ants while each ant had unique behavior.

Our mixture model identified a total of 7 different patterns in the example ant model (minimized BIC at 7 groups). Table 2 shows all the coefficients for the seven different patterns, labeled by their colors from Figure 6. The last two columns are the eigenvalues wherein we built matrices of the own and coupling effects in the same order as Equations (2) and (3) (First row: own predicting x, coup predicting x and second row: coup predicting y, own predicting y). The eigenvalue procedure allows us to account for when the coefficients do not directly represent the type of attractor dynamic due to the primary axes for the dynamic depictions being different from the variables used in the equations. When the eigenvalues are both real numbers and negative, the system depicts an attractor. When the eigenvalues are both real and positive, the system depicts a repeller. When one is positive and one is negative, the system depicts a saddle. Imaginary numbers instead depict cyclic behavior with complex numbers being a combination of cyclic and attractive/repulsive at the same time (Abraham and Shaw, 1983).

**Table 2. Unstandardized coefficients (and standard errors) for the seven group solution along with eigenvalues for the Ants model**.

Figure 6 is a topographical representation of the seven attractor dynamics patterns emergent in the ant behavior. Figure 6 was generated by using the estimated equations from the mixture model in conjunction with the adaptive Runge-Kutta algorithm from the deSolve package (Soetaert et al., 2010) in R (R Core Team, 2016) to estimate example trajectories iterated over time. In each case, values were chosen using the position means and variances extrapolating in all possible combinations of one standard deviation in X and Y and iterating the trajectories forward in time. Details on each pattern follow.

The blue, brown, and green patterns correspond to the food piles while the red pattern corresponds to the ant nest. The yellow pattern corresponds to searching an area where no food existed. The light blue captured the pattern of the ants converging in the middle of the trail as the pheromones were most intense there and the purple captured the dispersal after the food pile in the upper left had been fully collected (it was the first pile found in the simulation).

Notice how each pattern is captured through a different attractor dynamic. For example, the red nest pattern shows a repeller in which ants leave the location. If we capture each ant trail of food collection through the other patterns, then what primarily remains is the initial leaving from the nest. The blue and brown patterns, both corresponding to food piles, show cyclic properties (they have imaginary components to their eigenvalues). This is capturing the pattern of getting the food from the pile, bringing it to the nest and returning. The pattern corresponding to the lower left food pile was a saddle, however—attractive in one dimension and repulsive in the other. By having the set point far from the dynamic pattern, the saddle generated curved trails that could then be completed by feeding into other, already established, patterns.

Now, we link the agents to these patterns and to key system descriptions—in this case food depletion. Figure 7 shows the decline in the food piles as a function of time. Notably, the ants found the pile in the upper left first, followed by the lower left and then finally the middle right. We ran seven multilevel models treating the posterior probabilities of each pattern as the outcome as a function of the proportion of food remaining in each pile (a three predictor MLM). The fixed and random effects along with intraclass correlations (ICC) are in Table 3. All random effects were significantly non-zero suggesting that there was variability in their likely pattern as a function of the remaining food piles among the individual ants. The fixed effects can be interpreted as whether or not the likelihood of being in a pattern occurred where a positive sign meant that declines in a food pile corresponded to declines in the pattern and a negative sign meaning that declines in a food pile corresponded to increases in the pattern. Given the order of the food pile depletions, the pattern of effects can also roughly determine when the pattern was more prevalent.

**Figure 7. Time series plot of the amount of food available in each of the three food piles**. The legend describes where in the coordinate space a given food pile was located (see also Figure 3).

**Table 3. Unstandardized coefficients (and standard errors) and intraclass correlations from multilevel models predicting the posterior probabilities of being in each of the seven groups for the Ants models**.

The red pattern at the nest was likely when all the food sources were untouched and declined in likelihood as all the food piles declined, consistent with the ants initially leaving the nest to search. The blue, light blue, and purple patterns all associated with the upper left quadrant were all less likely when the last food pile was untouched, but only the purple (the theoretical dispersion after the food pile was depleted) was contingent upon the corresponding upper left food pile. The negative sign was indicative that declines in the first pile increased the likelihood of the purple dispersion pattern consistent with leaving the trail to find another food source once the food in the first pile was depleted. The green (corresponding to the lower left food pile) and brown (corresponding to the middle right food pile) patterns were predicted by all three food piles with negative coefficients suggesting that as any food depleted, these became more likely—consistent with these food piles being found later. Finally, the yellow pattern was only uniquely predicted by the middle right food pile depletion such that as the food pile declined, so did the likelihood of being in the search pattern. Given that as more ants found the last food pile, more converged on it. Once it depletes, however, fewer ants would be in this search pattern.

#### Baboons Navigation Data

So far, we have relied on simulations to illustrate how one can depict higher order emergent coordination for agent interactions using attractor dynamics. Our next two examples are derived from observed data. Figure 8 represents a solution from global positioning system (GPS) data collected from a troop of baboons at the De Hoop Nature Reserve in South Africa. Table 4 contains the coefficients and eigenvalues, again using colors to indicate correspondence. To collect this data, researchers recorded the positions of 14 adult baboons by holding a GPS device over or very close to each animal at different points over a 74 day period (data was made available by Bonnell et al., 2016; and further details of the original study can be found at Bonnell et al., 2017). Consistent, with the ants data, this example data is in an x/y coordinate space, but now in longitude and latitude. To facilitate estimation due to variability occurring in small decimal places, longitude and latitude were mean-centered and multiplied by 1,000.

**Table 4. Unstandardized coefficients (standard errors) and eigenvalues for the 10 pattern solution from the Baboon GPS data**.

Figure 8 illustrates that several of the patterns show cyclic behaviors. In fact, all the eigenvalues were negative with 7 of the 10 showing imaginary eigenvalues consistent with cyclic behaviors. Further, all patterns had at least one negative real eigenvalue suggesting that they all were attractive indicating a pattern of convergence for baboons. Figure 8 clearly shows that the patterns were not equally attractive, however, in that vector length differed dramatically when example trajectories were estimated. This can also be seen by the size of the eigenvalues where some were quite close to zero in their real number portion(s) while others were much smaller numbers approaching and surpassing negative one. Thus, some of these patterns were more stable clusters for the baboons while others were more loose associations around the shared longitude/latitude set point.

In their original work, Bonnell et al. (2017) evaluated whether the movement patterns of a focal individual baboon was influenced by the location of the troop as a collective or by the locations of specific influential members of the troop. Ultimately, their results showed evidence for both of these patterns. In some cases, the focal baboon's movement was highly influenced by the average movement location of the entire troop. In other cases, the focal baboon's movement was quite sensitive to the movements of the alpha female (F1) and the alpha male (M1). To link back to individual baboons, our results suggest a consistent pattern as illustrated in Figure 9 wherein we show the average posterior probabilities for each baboon illustrating which pattern would arguably influence a given baboon the majority of the time (again, colors correspond). Few distinctions existed between the female and dominant male baboons showing preference for the light green (cyclic attractor) and yellow (attractor) patterns. The dominant male (M1) showed slightly more preference for the magenta pattern (also an attractor). Thus, there is evidence of following the primary male baboon, but also one of a female majority. And yet in both cases these most common patterns represent the least attractive patterns (eigenvalues closest to zero) in that there is lots of wandering in comparison to the other patterns inferred from the GPS data.

**Figure 9. Average posterior probabilities associated with each equation group, by baboon**. F is for female and M is for male.

### Beyond Two Dimensions

As we move beyond two dimensions, it is difficult to make easy to read and meaningful maps of the data. However, our approach is not limited to two dimensions. By relying on the eigenvalues presented earlier, one can derive the higher order patterns to illustrate what is occurring without a means to draw them. Further, it also allows us to point out that any variables can be captured as attractor dynamics—they do not inherently need to be spatial, as illustrated by our next example.

Each new dimension corresponds to an additional equation. In the six-dimensional case that follows, we model six simultaneous equations where change in each variable is treated as the outcome from each equation. Further each variable at a given point in time is allowed to freely predict the changes in each equation. The matrix used to generate the eigenvalues is based on the coefficients where, as before, the main diagonal are the own effects and the off diagonals are the coupling relationships. Each matrix row corresponds to a different equation.

#### Affect in Families Data

To show a non-spatial example with more than 2 simultaneous change equations, we modeled positive and negative affect from the PANAS (Watson et al., 1988) taken from mothers, fathers, and adolescents from 252 families where the adolescent has type 1 diabetes. The data are taken from the Adolescents with Diabetes and Parents Together study where each family member completed a daily diary for 14 days (further study details can be found at Berg et al., 2009).

We extracted two stable patterns (three patterns would not properly converge and fit indices supported the two pattern solution). Table 5 provides the estimated coefficients. Notably, the eigenvalues were quite different between the two patterns. The first pattern generates all negative eigenvalues indicating that it forms one large six dimensional attractor (−0.709, −0.522, −0.455, −0.428, −0.285, −0.257). The second pattern, on the other hand had complex numbers for the first two eigenvalues suggesting cyclic behavior as a primary component (−0.613+0.027i, −0.613 −0.027i, −0.420, −0.343, −0.180, −0.157).

**Table 5. Unstandardized coefficients (and standard errors) from the two pattern solution for the Affect Daily Diary**.

Though we cannot draw a map to represent this higher order pattern, one way to represent the changes in the system is through a network diagram. Figures 10A,B shows only significant (alpha = 0.05, two-tailed) pathways between affect variables. The beginning of an arrow is value and the end of an arrow is change. Blue arrows represent negative relationships and brown ones are positive. Note between Figures 10A,B the connections between individuals breaks down substantially with the cyclic nature relating to the less connected network. The most noteworthy is the changing connections of father's affect to the mother and adolescent. It is noteworthy that these coefficients merely indicate prediction and thus any interpretation of causality would overstate the relationship. That said, fathers were clearly showing less connection in the second pattern.

**Figure 10. (A,B)** Two network diagrams that illustrate the two different equations. Beginning of arrows represent value at time t. Arrow heads represent change in value.

To link back to individual families, we built a multilevel model predicting the posterior probabilities for the first pattern as a function of diabetes risk for the adolescent on a given day. We use the variable *risk* as an easy to interpret indicator as to how well the adolescent was managing their diabetes on a given day. Risk is a rescaled version of daily blood glucose variability and level such that zero indicates perfect maintenance at doctor recommended levels and 100 indicates either going too high or too low repeatedly (both of which can be quite dangerous; see Kovatchev et al., 2006). Since posterior probabilities for a given data vector add to one, high probabilities of being in the first pattern inherently implies a low probability of being in the second. Table 6 contains the coefficients. At zero risk on a given day, families were equally likely to be in each pattern (the intercept is the posterior probability when risk was zero). As risk increased, however, families were more likely to fall into the second pattern. That is, on good days we see the more connected attractor pattern and on bad days the father appears less connected and the family affect adopts a cyclic attraction pattern instead.

**Table 6. Unstandardized coefficients (and standard errors) from multilevel model predicting the posterior probability of the first pattern as a function of Diabetes risk**.

## Discussion

Kelso (2009) posited that SCD “unites the spontaneous, self-organizing nature of coordination and the obviously directed, agent-like properties characteristic of animate nature into a single framework” (p. 1540). This logic matches with self-organization from agent-based models, and cases where many agents engage in social coordination, more generally. By connecting attractor dynamics modeling with cases where there are a range of agents and a range of outcomes allows for a generalized approach to quantifying the emergent patterns.

Through various examples, we illustrated that the attractor dynamics can be captured using a combination of difference/differential equation modeling and mixture modeling. Further, we showed that these attractor patterns and their occurrence could be linked with different outcomes. For the flocking model, we found sixteen attractor patterns of the agents' heading that converged on fewer attractors over time. For the ants model, we found seven dynamic patterns to depict their motion in a two-dimensional x/y space that roughly corresponded to qualitative depictions of rules the ants follow. For the baboon navigation data, we found ten patterns in two-dimensional longitudinal and latitudinal space in which the probability of exhibiting a particular attractor was contingent upon influential baboons in the troop (e.g., an alpha male). For families where an adolescent has type 1 diabetes, we found two patterns in a six dimensional affect space that corresponded to higher and lower levels of risk from the disease. By using the data from all the agents, the underlying topology is inclusive of all the agents. In the ants model, for example, not all ants illustrated being influenced by every pattern. Instead, ants can exist in a single pattern their entire time or move between them. Thus, the underlying map implied by the set of dynamic patterns generates an inclusive generalization both within and between agents that capitalizes on the most probable systems states over the duration of the observation period.

In each circumstance, the technique depicts the topological feature in terms of the implied patterns and the stability of those patterns. Whereas, the flocking model only contained attractors that varied in their set points (attractive headings) and their stabilities, the ants model illustrated all the common possible attractor dynamic patterns including attractors, repellers, saddles, and cycles.

The complexity of the underlying pattern is directly related to the number of dimensions. With a single dimension, attractor dynamics may only convey attractors and/or repellers. With two dimensions, cycles and saddles can be inferred. Beyond two dimensions, chaotic (strange) attractors are possible, though all currently known chaotic attractors require non-linear equation forms and the equations herein were restricted to linearity within each equation group. Thus, this is a limitation of the technique provided.

In each case, we then linked the quantification back to the individual agents. Through mixture modeling we did this by outputting the posterior probabilities. These probabilities are the probability that a given data vector is under the influence of a given dynamic pattern, the probabilities for a given vector sum to one across all the possible patterns. Therefore, these probabilities maintain the data dependency we inherently ignored in the estimation for the dynamic patterns themselves. We therefore always either examined the probabilities at a collapsed agent level (e.g., averages) or through multilevel modeling wherein the dependencies could be properly taken into account. In each case, it could be linked to possible variables of interest used to depict the system. For the headings, this was illustrated with time in that attractors should collapse as time goes on. For the ants model, this was illustrated through food supply. For the baboons, this was illustrated through the location of the alpha male and the females. For affect in families where the adolescent had type 1 diabetes, it was illustrated with the diabetes risk exhibited that day. In all, this allows one to link the higher order patterns back to meaningful outcomes that characterize when agents behave in certain ways or exhibit theoretically important states.

In the spatial examples, we utilized variables that depicted the spatial movement. As an initial foray into understanding attractor dynamics, thinking spatially helps make the concepts more intuitive. But, ultimately, these concepts can be applied in many contexts where relationships are not inherently spatial. Being able to think about the spatial analogs helps ground what is being observed, but does not inherently limit the domains in which attractor dynamics can be examined.

Further, individual equation parameters do not always align with the system depiction graphically or through the eigenvalue procedure. In the ants example, this had to do with the reliance on diagonal movement of the ants. By depicting the system through an equation of x and an equation of y, we mask diagonal movement—it is really a straightforward combination of the two dimensions rather than showing some independence. More generally, the coefficients are under an assumption that the dimensions chosen are the primary dimensions for depicting the changes occurring in the system. The eigenvalue procedure bypasses this assumption by instead capitalizing on dimensions that maximize the strength of the attractor dynamics.

Once we moved beyond two dimensions, the eigenvalue procedure becomes even more valuable. There is no easy way to graphically “see” the implied dynamic, but the sign and distinctions between real and imaginary portions elucidate the attractor pattern. In practice, anytime we model a system with two or more equations we should adopt the eigenvalue procedure as a means to understand the higher order pattern in addition to any interpretations applied to the individual coefficients themselves. For example, it is common to interpret coupling coefficients as the push/pull of one variable upon another. However, this fails to capture what pattern the push-pull creates as their interpretation is under an assumption that we somehow picked ideal dimensions to represent them. Locally, the coefficients maintain their meaning, but we cannot extrapolate the more global pattern of which they are a part.

In regards to equation identification, the technique is not without its limitations. The choice of slicing up the data into a series of locally linear equations is an imperfect method for capturing non-linear dynamic models. Specifically, non-linear dynamic models can have both multistability in which more than one pattern is stable simultaneously and cases where variables differentiate when one pattern is or is not accessible. By slicing up the data into a series of locally linear equations through mixture modeling, these two circumstances are difficult to distinguish. One can begin to distinguish these circumstances by attempting to predict the posterior probabilities. However, ultimately multistability is distinguished by states being probable despite nothing differentiating them (or when the dimensions being examined are all that differentiate them). That is, multistability would occur under a lack of being able to predict differences of when agents would be in one or the other. Thus, this approach provides a limited potential for knowing when multistability exists as opposed to having some variable differentiate them. We may never examine the “right” variable or are instead in the situation of arguing a null finding to support the multistable case.

In contrast, it is possible through a cusp catastrophe model in conjunction with multilevel modeling, for example, to allow for differentiating variables (also known as control parameters) without their identification (Butner et al., 2014b), though knowing which scenario you are observing requires examination of many more qualities than discussed herein (Gilmore, 1981). Further, manifolds (the surfaces implied by topological equations) are smooth, while the mixture modeling approach is more patchwork. We do not know the reach of a given attractor dynamic—we chose to represent each dynamic through one standard deviation in each direction from the means when we utilized the Runge-Kutta algorithm to graph plausible trajectories. Notably the means and standard deviations are specific to each dynamic pattern (allowing some to be large and others to be smaller). However, the boundaries of one pattern to another are truly unknown, requiring some inference.

Notably, SCD has tended to rely on cyclical descriptions to model the rhythmic coordination of social agents. While the modeling approach illustrated herein allows for cycles, it does not assume their existence. The direct equation link is that SCD generally functions with second order equations where the second derivatives (acceleration or change in velocity) are treated as the outcomes. Within our structural equation model, it would be analogous to building a quadratic growth model on Toeplitz data where the quadratic growth latent variable would be the second derivative predicted by the other two latent variables (Butner and Story, 2010). Moving to a second order model automatically implies two dimensions and thus generates cycles. However, it is not without a cost. Specifically, second order modeling in this form assumes that the set point of the cycles must equal zero. Overcoming this assumption is currently something under consideration for modeling dynamic patterns and once resolved will unite these approaches more generally.

## Conclusion

Understanding how large-scale, multi-agent social systems coordinate is challenging and complex. In part, the challenge is due to the fact that there are so many agents, system components, and potential system states that can become coordinated; all of which may change over time (Van Orden et al., 2003). These many components interact generating higher order system behavior that is emergent and dynamic. However, knowing the “Dynamics demystifies…emergence” and it can also provide “basic laws for a quantitative description of phenomena that are observed” (Kelso, 2009; p. 1540). As such, we have expanded on work in SCD by demonstrating the utility of modeling the attractor dynamics of several systems to characterize their higher-order behavioral patterns and showed how these patterns varied over time and could be linked to meaningful aspects of the systems.

Within domains, such as agent based modeling, qualitative depictions of higher-order patterns are often known, but not quantitatively modeled. In SCD, phenomena can be non-rhythmic, and yet dynamically coordinated. They can exhibit stability and multistability. Thus, using attractor dynamic descriptions along with statistical innovations, such as mixture modeling, provide a reasonable solution to understanding the large-scale, multi-agent social coordination. Characterizing the higher order properties of the system in this way forms a foundation for examining the emergent patterns through time in either a confirmatory or exploratory manner. This same technique, as we have shown, can be utilized with simulated as well as observational data.

It is our aim that we recognize that we study systems that are inherently open systems (even though simulations are often closed). By examining part of the system (the variables we measure), unobserved aspects of the system function as perturbations to the system. Thus, a system depicting families is open because we are only examining some of the variables involved. To understand how agents exhibit emergent self-organization and coordination, we have advanced a general quantification that can be applied to a range of social systems, such as two individuals that form a couple up to a crowd's behavior. We hope that the widely applicable techniques will be adopted to advance scientific understanding of SCD.

## Ethics Statement

The Adolescents with Diabetes and Parents Together project was conducted in accordance with the recommendations of the University of Utah Institutional Review Board. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the University of Utah Institutional Review Board.

## Author Contributions

All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

This study was partially supported by grant R01 DK063044-01A1 from the National Institute of Diabetes and Digestive Kidney Diseases awarded to Deborah Wiebe (PI) and Cynthia Berg (co-PI). We also would like the thank Tyler Bonnell for making the baboon data publicly available and for answering our inquiries regarding the data set.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fpsyg.2017.00380/full#supplementary-material

## References

Abraham, R. H., and Shaw, C. D. (1983). *Dynamics-the Geometry of Behavior*. Santa Cruz, CA: Aerial Press.

Arabanel, H. D. I., Brown, R., and Kennel, M. B. (1992). Local Lyapunov exponents computed from observed data. *J. Nonlin. Sci.* 2, 343–365. doi: 10.1007/BF01208929

Bauer, D. J., and Curran, P. J. (2003). Distributional assumptions of growth mixture models: implications for overextraction of latent trajectory classes. *Psychol. Methods* 8, 338–363. doi: 10.1037/1082-989X.8.3.338

Bedau, M. A., and Humphreys, P. E. (2008). *Causation and Autonomy in Weak Emergence, Emergence: Contemporary Readings in Philosophy and Science*. Cambridge, MA: MIT Press.

Berg, C. A., Skinner, M., Ko, K., Butler, J. M., Palmer, D. L., Butner, J., et al. (2009). The fit between stress appraisal and dyadic coping in understanding perceived coping effectiveness for adolescents with type 1 diabetes. *J. Family Psychol.* 23, 521–530. doi: 10.1037/a0015556

Boker, S. M., Deboeck, P. R., Edler, C., and Keel, P. K. (2010). “Generalized local linear approximation of derivatives from time series,” in *Statistical Methods for Modeling Human Dynamics: An Interdisciplinary Dialogue*, eds S.-M. Chow, and E. Ferrar (Boca Raton, FL: Taylor and Francis), 161–178.

Boker, S. M., and Laurenceau, J. P. (2006). “Dynamical systems modeling: an application to the regulation of intimacy and disclosure in marriage,” in *Models Intensive Longitudinal Data*, eds T. A. Walls and J. L. Schafer (New York, NY: Oxford University Press), 195–218.

Bonnell, T. R., Henzi, S. P., and Barrett, L. (2016). Data from: direction matching for sparse movement data sets: determining interaction rules in social groups. *Dryad Digital Repository.* doi: 10.5061/dryad.kv2kh

Bonnell, T. R., Henzi, S. P., and Barrett, L. (2017). Direction matching for sparse movement data sets: determining interaction rules in social groups. *Behav. Ecol.* 28, 193–203. doi: 10.1093/beheco/arw145

Butner, J. E., Berg, C. A., Baucom, B. R., and Wiebe, D. J. (2014a). Modeling coordination in multiple simultaneous latent change scores. *Multivar. Behav. Res.* 49, 554–570. doi: 10.1080/00273171.2014.934321

Butner, J. E., Dyer, H. L., Malloy, T. S., and Kranz, L. V. (2014b). Uncertainty in cost performance as a function of the cusp catastrophe in the NASA program performance management system. *Nonlin. Dyn. Psychol. Life Sci.* 18, 397–417.

Butner, J. E., Gagnon, K. T., Geuss, M. N., Lessard, D. A., and Story, T. N. (2015). Utilizing topology to generate and test theories of change. *Psychol. Methods* 20:1. doi: 10.1037/a0037802

Butner, J., and Story, T. N. (2010). “Oscillators with differential equations,” in *Nonlinear Dynamical Systems Analysis for the Behavioral Sciences: Real Data*, eds S. Guastello and R. Gregson (Boca Raton, FL: CRC Press), 367–400.

Cohen, J., Cohen, P., West, S. G., and Aiken, L. S. (2003). *Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences*. Mahwah, NJ: Lawrence Erlbaum Associates.

Couzin, I. D., and Krause, J. (2003). Self-organization and collective behavior in vertebrates. *Adv. Study Behav.* 32, 1–75. doi: 10.1016/S0065-3454(03)01001-5

Duarte, R., Araújo, D., Correia, V., Davids, K., Marques, P., and Richardson, M. J. (2013). Competing together: assessing the dynamics of team–team and player–team synchrony in professional association football. *Hum. Mov. Sci.* 32, 555–566. doi: 10.1016/j.humov.2013.01.011

Enders, C. K. (2006). “Analyzing structural equation models with missing data,” in *Structural Equation Modeling: A Second Course*, eds G. R. Hancock and R. O. Mueller (Greenwich, CT: Information Age Publishing), 313–342.

Enders, C. K., and Tofighi, D. (2007). Centering predictor variables in cross-sectional multilevel models: a new look at an old issue. *Psychol. Methods* 12:121. doi: 10.1037/1082-989X.12.2.121

Gilmore, R. (1981). *Catastrophe Theory for Scientists and Engineers*. New York, NY: John Wiley and Sons.

Gottman, J. M., Murray, J. D., Swanson, C. C., Tyson, R., and Swanson, K. R. (2002). *The Mathematics of Marriage: Dynamic Nonlinear Models*. Cambridge, MA: MIT Press.

Haken, H., Kelso, J. S., and Bunz, H. (1985). A theoretical model of phase transitions in human hand movements. *Biol. Cybern.* 51, 347–356. doi: 10.1007/BF00336922

Halley, J., and Winkler, D. A. (2008). Classification of emergence and its relation to self-organization. *Complexity* 13, 10–15. doi: 10.1002/cplx.20216

Henson, J. M., Reise, S. P., and Kim, K. H. (2007). Detecting mixtures from structural model differences using latent variable mixture modeling: a comparison of relative model fit statistics. *Struct. Equation Model. Multidiscipl. J.* 14, 202–226. doi: 10.1080/10705510709336744

Jung, T., and Wickrama, K. A. S. (2008). An introduction to latent class growth analysis and growth mixture modeling. *Soc. Pers. Psychol. Compass* 2, 302–317. doi: 10.1111/j.1751-9004.2007.00054.x

Kelso, J. A. (1995). *Dynamic Patterns: The Self Organization of Brain and Behavior*. Cambridge: MIT Press.

Kelso, J. A. (2009). “Coordination dynamics,” in *Encyclopedia of Complexity and Systems Science*, ed R. A. Meyers (New York, NY: Springer), 1537–1565

Kovatchev, B. P., Otto, E., Cox, D., Gonder-Frederick, L., and Clarke, W. (2006). Evaluation of a new measure of blood glucose variability in diabetes. *Diabetes Care* 29, 2433–2438. doi: 10.2337/dc06-1085

Lo, Y., Mendell, N. R., and Rubin, D. B. (2001). Testing the number of components in a normal mixture. *Biometrika* 88, 767–778. doi: 10.1093/biomet/88.3.767

Lubke, G., and Neale, M. C. (2006). Distinguishing between latent classes and continuous factors: resolution by maximum likelihood?. *Multivariate Behav. Res.* 41, 499–532. doi: 10.1207/s15327906mbr4104_4

McArdle, J. J. (2009). Latent variable modeling of differences and changes with longitudinal data. *Annu. Rev. Psychol.* 60, 577–605. doi: 10.1146/annurev.psych.60.110707.163612

Muthen, B. (2001). “Latent variable mixture modeling,” in *New Developments and Techniques in Structural Equation Modeling*, eds G. A. Marcoulides and R. E. Schumacker (Mawwah, NJ: Lawrence Erlbaum Associates), 1–33.

Muthén, L. K., and Muthén, B. O. (1998-2012). *Mplus User's Guide, 7th Edn*. Los Angeles, CA: Muthén and Muthén.

Nagin, D., and Tremblay, R. E. (1999). Trajectories of boys' physical aggression, opposition, and hyperactivity on the path to physically violent and nonviolent juvenile delinquency. *Child Dev.* 70, 1181–1196. doi: 10.1111/1467-8624.00086

Nylund, K. L., Asparouhov, T., and Muthén, B. O. (2007). Deciding on the number of classes in latent class analysis and growth mixture modeling: a Monte Carlo simulation study. *Struct. Equation Model.* 14, 535–569. doi: 10.1080/10705510701575396

Oullier, O., De Guzman, G. C., Jantzen, K. J., Lagarde, J., and Scott Kelso, J. A. (2008). Social coordination dynamics: measuring human bonding. *Soc. Neurosci.* 3, 178–192. doi: 10.1080/17470910701563392

Oullier, O., and Kelso, J. A. S. (2009). “Social coordination, from the perspective of coordination dynamics,” in *Encyclopedia of Complexity and Systems Science*, ed R. A. Meyers (New York, NY: Springer), 8198–8213

Prigogene, I., and Stangers, I. (1984). *Order out of Chaos: Man's New Dialog with Nature*. New York, NY: Bantam Books.

R Core Team (2016). *R: A Language and Environment for Statistical Computing*. Vienna: R Foundation for Statistical Computing

Richardson, M. J., Dale, R., and Marsh, K. L. (2014). “Complex dynamical systems in social and personality psychology,” in *Handbook of Research Methods in Social and Personality Psychology*, eds H. T. Reis and C. M. Judd (New York, NY: Cambridge University Press), 253.

Richardson, M. J., Garcia, R. L., Frank, T. D., Gergor, M., and Marsh, K. L. (2012). Measuring group synchrony: a cluster-phase method for analyzing multivariate movement time-series. *Front. Physiol.* 3:405. doi: 10.3389/fphys.2012.00405

Richardson, M. J., Marsh, K. L., Isenhower, R. W., Goodman, J. R., and Schmidt, R. C. (2007). Rocking together: dynamics of intentional and unintentional interpersonal coordination. *Hum. Mov. Sci.* 26, 867–891. doi: 10.1016/j.humov.2007.07.002

Schmidt, R. C., Carello, C., and Turvey, M. T. (1990). Phase transitions and critical fluctuations in the visual coordination of rhythmic movements between people. *J. Exp.Psychol. Hum. Percept. Perform.* 16:227. doi: 10.1037/0096-1523.16.2.227

Schmidt, R. C., and O'Brien, B. (1997). Evaluating the dynamics of unintended interpersonal coordination. *Ecol. Psychol.* 9, 189–206. doi: 10.1207/s15326969eco0903_2

Sclove, S. L. (1987). Application of model-selection criteria to some problems in multivariate analysis. *Psychometrika* 52, 333–343. doi: 10.1007/BF02294360

Soetaert, K., Petzoldt, T., and Setzer, R. W. (2010). Solving differential equations in R: package deSolve. *J. Stat. Softw.* 33, 1–25. doi: 10.18637/jss.v033.i09

Strogatz, S. H. (2014). *Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering*. Boulder, CO: Westview Press.

Theraulaz, G., and Bonabeau, E. (1999). A brief history of stigmergy. *Artif. Life* 5, 97–116. doi: 10.1162/106454699568700

Van Orden, G. C., Holden, J. G., and Turvey, M. T. (2003). Self-organization of cognitive performance. *J. Exp. Psychol. Gen.* 132:331. doi: 10.1037/0096-3445.132.3.331

Vuong, Q. H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. *Econometrica* 57, 307–333. doi: 10.2307/1912557

Watson, D., Clark, L. A., and Tellegen, A. (1988). Development and Validation of brief measures of positive and negative affect: the PANAS scales. *J. Pers. Soc. Psychol.* 54, 1063–1070. doi: 10.1037/0022-3514.54.6.1063

Wilensky, U. (1997). *NetLogo Ants Model*. Evanston, IL: Center for Connected Learning and Computer-Based Modeling, Northwestern University. Available Online at: http://ccl.northwestern.edu/netlogo/models/Ants

Wilensky, U. (1998). *NetLogo Flocking Model.* Evanston, IL: Center for Connected Learning and Computer-Based Modeling, Northwestern University.

Wilensky, U. (1999). *NetLogo*. Evanston, IL: Center for Connected Learning and Computer-Based Modeling, Northwestern University. Available online at: http://ccl.northwestern.edu/netlogo/

Keywords: dynamical systems, social coordination dynamics, multi-agent coordination, attractors, agent-based modeling

Citation: Butner JE, Wiltshire TJ and Munion AK (2017) Modeling Multi-Agent Self-Organization through the Lens of Higher Order Attractor Dynamics. *Front. Psychol.* 8:380. doi: 10.3389/fpsyg.2017.00380

Received: 08 August 2016; Accepted: 27 February 2017;

Published: 20 March 2017.

Edited by:

Joanna Raczaszek-Leonardi, University of Warsaw, PolandReviewed by:

Robert J. Lowe, University of Skövde, SwedenHecke Schrobsdorff, Max Planck Institute for Dynamics and Self Organization, Germany

Copyright © 2017 Butner, Wiltshire and Munion. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Jonathan E. Butner, jonathan.butner@psych.utah.edu