Systems with non-linear dynamics frequently exhibit emergent system behavior, which is important to find and specify rigorously to understand the nature of the modeled phenomena. Through this analysis, it is possible to characterize phenomena such as how systems assemble or dissipate and what behaviors lead to specific final system configurations. Agent Based Modeling (ABM) is one of the modeling techniques used to study the interaction dynamics between a system's agents and its environment. Although the methodology of ABM construction is well understood and practiced, there are no computational, statistically rigorous, comprehensive tools to evaluate an ABM's execution. Often, a human has to observe an ABM's execution in order to analyze how the ABM functions, identify the emergent processes in the agent's behavior, or study a parameter's effect on the system-wide behavior. This paper introduces a new statistically based framework to automatically analyze agents' behavior, identify common system-wide patterns, and record the probability of agents changing their behavior from one pattern of behavior to another. We use network based techniques to analyze the landscape of common behaviors in an ABM's execution. Finally, we test the proposed framework with a series of experiments featuring increasingly emergent behavior. The proposed framework will allow computational comparison of ABM executions, exploration of a model's parameter configuration space, and identification of the behavioral building blocks in a model's dynamics.

Agent based modeling (ABM) is a popular tool to study the nature of emergent phenomena with complex, stochastic behavior. The analysis of ABM model behavior is often left to a human since the agents have stochastic, non-linear behavior similar to the dynamics of the phenomena the model was built to understand. We present a statistically based methodology to analyze agent behaviors to identify common, system-wide patterns of behavior using only the agent behaviors recorded during the model execution. The methodology is automatic, tunable, and problem-independent. In addition, no knowledge of the model parameter space is assumed, and the methodology provides multi-scale information about the model's behavior. It can be used as a fitness function in the stochastic search algorithms, as heuristics for tuning the model parameters, to analyze emergent behavior in the physical systems, or to find a high diversity set of the model parameter configurations that satisfy model objectives.

I. INTRODUCTION Section: Choose Top of page ABSTRACT I. INTRODUCTION << II. BACKGROUND III. METHODOLOGY IV. RESULTS AND ANALYSIS V. CONCLUSION, DISCUSSION... CITING ARTICLES systems have been studied in many fields in an attempt to understand how these systems assemble, function, transition, and scale. 14,33,38 14. S. Kauffman, The Origins of Order: Self Organization and Selection in Evolution ( Oxford University Press , 1993). 33. S. H. Strogatz, Nonlinear Dynamics and Chaos: With applications to Physics, Biology, Chemistry, and Engineering ( Westview Press , 2014). 38. A. Wuensche, Exploring Discrete Dynamics ( Luniver Press , 2011). 34 34. A. Tkachenko, “ Inverse problem in self-assembly ,” in APS Meeting Abstracts (2012), p. 46002. 1,10,24,25 1. R. Axelrod, The Complexity of Cooperation: Agent-based Models of Competition and Collaboration. Princeton Studies in Complexity ( Princeton University Press , 1997). Modeling civil violence: An agent-based computational approach ,” Proc. Natl. Acad. Sci. 99, 7243– 7250 (2002). 10. J. M. Epstein, “,” Proc. Natl. Acad. Sci., 7243–(2002). https://doi.org/10.1073/pnas.092080199 24. S. Ortman, Winds from the North: Tewa Origins and Historical Anthropology ( The University of Utah Press , Salt Lake City, UT , 2012). Settlement scaling and increasing returns in an ancient society ,” Sci. Adv. 1, e1400066 (2015). 25. S. G. Ortman, A. H. F. Cabaniss, J. O. Sturm, and L. M. A. Bettencourt, “,” Sci. Adv., e1400066 (2015). https://doi.org/10.1126/sciadv.1400066 dynamics of a software product as it is being built; 13 13. V. Honsel, “ Statistical learning and software mining for agent based simulation of software evolution ,” in 2015 IEEE/ACM 37th IEEE International Conference on Software Engineering (ICSE) (2015), Vol. 2, pp. 863– 866 . 5,26 5. S. Chandrasekaran and D. Hougen, “ Swarm intelligence for cooperation of bio-nano robots using quorum sensing ,” in Bio Micro and Nanosystems Conference, BMN'06 (2006), p. 104. Nanorobot: A versatile tool in nanomedicine ,” J. Drug Targeting 14(2), 63– 67 (2006). 26. G. M. Patel, G. C. Patel, R. B. Patel, J. K. Patel, and M. Patel, “,” J. Drug Targeting(2), 63–(2006). https://doi.org/10.1080/10611860600612862 Emergent processes in complexhave been studied in many fields in an attempt to understand how theseassemble, function, transition, and scale.For example, physicists working on soft-condensed matter rely on particle interactions to build arbitrary tertiary structures;anthropologists and sociologists study the rise and fall of cultures and social movements;managers of large software projects study the structure, evolution, and personalof a software product as it is being built;and engineers design biologically inspired nano-scale swarm robots for applications in medicine. Modeling (ABM) is a common approach to study systems with emergent properties. 11–13 11. G. Gilbert, Agent-Based Models ( SAGE Publications , 2008), Vol. 153. A standard protocol for describing individual-based and agent-based models ,” Ecol. Modell. 198(12), 115– 126 (2006). 12. V. Grimm, U. Berger, F. Bastiansen, S. Eliassen, V. Ginot, J. Giske, J. Goss-Custard, T. Grand, S. K. Heinz, G. Huse, A. Huth, J. U. Jepsen, C. Jrgensen, W. M. Mooij, B. Mller, G. Peer, C. Piou, S. F. Railsback, A. M. Robbins, M. M. Robbins, E. Rossmanith, N. Rger, E. Strand, S. Souissi, R. A. Stillman, R. Vab, U. Visser, and D. L. DeAngelis, “,” Ecol. Modell.(12), 115–(2006). https://doi.org/10.1016/j.ecolmodel.2006.04.023 13. V. Honsel, “ Statistical learning and software mining for agent based simulation of software evolution ,” in 2015 IEEE/ACM 37th IEEE International Conference on Software Engineering (ICSE) (2015), Vol. 2, pp. 863– 866 . ABM is used to study how system-wide properties emerge from local interactions of a system's components. The general principles of how to select, describe, and parametrize the significant phenomena of a system have been widely studied and have resulted in a number of model construction protocols. 9,12,13,17,27 9. J. M. Epstein, Agent-Based Computational Models and Generative Social Science ( John Wiley & Sons, Inc. , 1999), Vol. 4(5), pp. 41– 46 . A standard protocol for describing individual-based and agent-based models ,” Ecol. Modell. 198(12), 115– 126 (2006). 12. V. Grimm, U. Berger, F. Bastiansen, S. Eliassen, V. Ginot, J. Giske, J. Goss-Custard, T. Grand, S. K. Heinz, G. Huse, A. Huth, J. U. Jepsen, C. Jrgensen, W. M. Mooij, B. Mller, G. Peer, C. Piou, S. F. Railsback, A. M. Robbins, M. M. Robbins, E. Rossmanith, N. Rger, E. Strand, S. Souissi, R. A. Stillman, R. Vab, U. Visser, and D. L. DeAngelis, “,” Ecol. Modell.(12), 115–(2006). https://doi.org/10.1016/j.ecolmodel.2006.04.023 13. V. Honsel, “ Statistical learning and software mining for agent based simulation of software evolution ,” in 2015 IEEE/ACM 37th IEEE International Conference on Software Engineering (ICSE) (2015), Vol. 2, pp. 863– 866 . On legitimacy feedback mechanisms in agent-based modeling of civil violence ,” Int. J. Intelligent Syst. 31(2), 106– 127 (2016). 17. C. Lemos, R. J. Lopes, and H. Coelho, “,” Int. J. Intelligent Syst.(2), 106–(2016). https://doi.org/10.1002/int.21747 27. S. F. Railsback and V. Grimm, Agent-Based and Individual-Based Modeling: A Practical Introduction ( Princeton University Press , 2011). ABM is, the execution of the model often has non-linear dynamics that are as difficult to understand as the studied physical phenomena itself. The current lack of computational tools for ABM analysis explains why model analysis is frequently simplified to either binary decision or a simple statistical distribution used to measure if the simulation goal was reached. How the simulation goal was reached or what the model's specific behavior was, on the other hand, is usually delegated to a human observing the model's execution. The manual analysis of the model's behavior is also not feasible when evaluating the model's behavior on a combinatorial number of all possible parameter configurations that control the model. Independent of the problem domain, Agent Basedis a common approach to studywith emergent properties.is used to study how system-wide properties emerge from local interactions of a system's components. The general principles of how to select, describe, and parametrize the significant phenomena of ahave been widely studied and have resulted in a number ofconstruction protocols.Regardless of how simple or sophisticated anis, the execution of theoften hasthat are as difficult to understand as the studied physical phenomena itself. The current lack of computational tools forexplains whyis frequently simplified to either binary decision or a simple statistical distribution used to measure if the simulation goal was reached.the simulation goal was reached orthespecific behavior was, on the other hand, is usually delegated to a human observing theexecution. The manualof thebehavior is also not feasible when evaluating thebehavior on a combinatorial number of all possible parameter configurations that control the 4 4. M. Cenek and S. Dahl, see https://github.com/mcenek/BehavioralSpacesFramework for Geometry of behavioral spaces framework. analyze agents' behaviors and identify common patterns of system-wide behavior that are present during the model's execution. This analysis can be used as a tool to identify the behavioral building blocks of individual agent that led to the emergent system dynamics, to use as a fitness function for the stochastic search algorithms to find high diversity configuration of the model's parameters that satisfy the simulation goals, 32 32. F. J. Stonedahl, “ Genetic algorithms for the exploration of parameter spaces in agent-based models ,” Ph.D. thesis ( Northwestern University, Evanston, IL, USA , 2011), AAI3489404. model's executions on different parameter configurations. 3 3. B. Calvez and G. Hutzler, “ Automatic tuning of agent-based models using genetic algorithms ,” in Proceedings of the 6th International Workshop on Multi-Agent-Based Simulation MABS 2005 (Springer, 2006), pp. 41– 57 . In this paper, we present a statistically based methodologyto computationallyagents' behaviors and identify common patterns of system-wide behavior that are present during theexecution. Thiscan be used as a tool to identify the behavioral building blocks of individual agent that led to the emergentto use as a fitness function for thesearch algorithms to find high diversity configuration of theparameters that satisfy the simulation goals,and to use as a comparison metric of theexecutions on different parameter configurations.

II. BACKGROUND Section: Choose Top of page ABSTRACT I. INTRODUCTION II. BACKGROUND << III. METHODOLOGY IV. RESULTS AND ANALYSIS V. CONCLUSION, DISCUSSION... CITING ARTICLES ABMs are used to understand the connection between an individual agent's actions and a system's overall behavior in order to solve problems such as system-level control, analysis of coupled agent-to-agent or agent-to-environment system dynamics, or analysis of a system's sensitivity to parameter(s) changes. To achieve these goals, (1) a model cannot be examined a priori without running it and (2) the model's dynamics must be analyzed by tools that are multi-scale, quantitative, and problem domain-independent. model's parameter space and its optimization are often based on measuring the divergence between the initial model and its behavior after certain parameters are altered. Miller uses a set of non-linear equations as tests to explore parameter configurations that result in significant changes to a system's dynamics. 22 Active nonlinear tests (ants) of complex simulation models ,” Manage. Sci. 44(6), 820– 830 (1998). 22. J. H. Miller, “,” Manage. Sci.(6), 820–(1998). https://doi.org/10.1287/mnsc.44.6.820 et al. summarize the mathematical and statistical tools that describe the dynamics of the collective behavior of multi-agent and swarm systems. 18,19 18. K. Lerman and A. Galstyan, “ A general methodology for mathematical analysis of multi-agent systems ,” in ISI-TR-529, USC Information Sciences Institute, Marina del Rey, CA (2001). 19. K. Lerman, A. Martinoli, and A. Galstyan, “ A review of probabilistic macroscopic models for swarm robotic systems ,” in Swarm Robotics ( Springer , 2005), pp. 143– 152 . stochastic search tools that explore the ABM's parameter search-space. 3,32 3. B. Calvez and G. Hutzler, “ Automatic tuning of agent-based models using genetic algorithms ,” in Proceedings of the 6th International Workshop on Multi-Agent-Based Simulation MABS 2005 (Springer, 2006), pp. 41– 57 . 32. F. J. Stonedahl, “ Genetic algorithms for the exploration of parameter spaces in agent-based models ,” Ph.D. thesis ( Northwestern University, Evanston, IL, USA , 2011), AAI3489404. Existing macro-level analytic tools used to examine aparameter space and its optimization are often based on measuring the divergence between the initialand its behavior after certain parameters are altered. Miller uses a set ofequations as tests to explore parameter configurations that result in significant changes to a system'sLerman. summarize the mathematical and statistical tools that describe theof the collective behavior of multi-agent and swarmFinally, Stonedahl, Calvez, and others extend and implement existing research as automaticsearch tools that explore theparameter search-space. analyses that use agent based behavior to control and analyze system dynamics are either domain specific, relying on parameter experimentation, or only apply to models with a tipping-point behavior. Spears et al. build a physics-based model for control and analysis of vehicle swarms. 31 Distributed, physics-based control of swarms of vehicles ,” Auton. Robots 17(2–3), 137– 162 (2004). 31. W. Spears, D. Spears, J. Hamann, and R. Heil, “,” Auton. Robots(2–3), 137–(2004). https://doi.org/10.1023/B:AURO.0000033970.96785.f2 et al. generalized existing analytic framework and described agents as finite state machines and their interactions as Markov processes to characterize the behavior of multi-agent systems with the threshold behavior. 23 23. D. Miner and M. Desjardins, Predicting and controlling system-level parameters of multi-agent systems, in Complex Adaptive Systems and the Threshold Effect: View from the National and Social Sciences: AAAI Fall Symposium, 5–7 November 2009 (2009), Vol. FS-09-04, pp. 92–95. system dynamics from individual agent behaviors, but the resulting methodology is too lossy for analysis of both macro/system and micro/agent level dynamics. 20 20. J. T. Lizier, The Local Information Dynamics of Distributed Computation in Complex Systems, Springer Theses ( Springer , Berlin/Heidelberg , 2013). analysis tools of Van et al. measure the expected behavior of agents in Kugler–Turvey's model of an ant colony. 15,35 15. P. N. Kugler and M. T. Turvey, Information, Natural Law, and the Self-Assembly of Rhythmic Movement ( Lawrence Erlbaum Associates, Inc. , 1987). 35. H. Van Dyke Parunak and S. Brueckner, “ Entropy and self-organization in multi-agent systems ,” in Proceedings of the Fifth International Conference on Autonomous Agents (2001), pp. 124– 130 . et al. used information storage and information modification filters from Lizier et al. to describe the system-wide information dynamics in swarm robotics. 20,21,36,37 20. J. T. Lizier, The Local Information Dynamics of Distributed Computation in Complex Systems, Springer Theses ( Springer , Berlin/Heidelberg , 2013). Local measures of information storage in complex distributed computation ,” Inf. Sci. 208, 39– 54 (2012). 21. J. T. Lizier, M. Prokopenko, and A. Y. Zomaya, “,” Inf. Sci., 39–(2012). https://doi.org/10.1016/j.ins.2012.04.016 36. X. R. Wang, J. M. Miller, J. T. Lizier, M. Prokopenko, and L. F. Rossi, “ Measuring information storage and transfer in swarms ,” in Proceedings of the Eleventh European Conference on the Synthesis and Simulation of Living Systems (2011), pp. 838– 845 . Quantifying and tracing information cascades in swarms ,” PLoS One 7(7), e40084 (2012). 37. X. R. Wang, J. M. Miller, J. T. Lizier, M. Prokopenko, and L. F. Rossi, “,” PLoS One(7), e40084 (2012). https://doi.org/10.1371/journal.pone.0040084 Existing micro-levelthat use agent based behavior to control andare either domain specific, relying on parameter experimentation, or only apply towith a tipping-point behavior. Spears. build a physics-basedfor control andof vehicle swarms.Miner. generalized existing analytic framework and described agents as finite state machines and their interactions as Markov processes to characterize the behavior ofwith the threshold behavior.Shannon's entropy based tools were proposed to provide information aboutfrom individual agent behaviors, but the resulting methodology is too lossy forof both macro/system and micro/agent levelInformationtools of Van. measure the expected behavior of agents in Kugler–Turvey'sof an ant colony.Similarly, Wang. used information storage and information modification filters from Lizierto describe the system-wide informationin swarm robotics. We present a methodology that complements previous research by (1) providing a domain-independent approach that only analyzes the recorded agent behaviors without knowledge of the model's parameter space and (2) builds a network abstraction of common agent behaviors that captures both micro- and macrolevel dynamics. ABM arrives at its final configuration or how it solves a given problem, we define information in an ABM as the behavioral patterns of the system's agents, and computation as a collection of behaviors along with the transition probabilities of the system's agents changing from one behavior to another. The model's computational mechanics can then be viewed as a state-space transition graph of system's constituent behaviors that result in a model achieving its goal. 28,30 Computational mechanics: Pattern and prediction, structure and simplicity ,” J. Stat. Phys. 104(3–4), 817– 879 (2001). 28. C. Shalizi and J. Crutchfield, “,” J. Stat. Phys.(3–4), 817–(2001). https://doi.org/10.1023/A:1010388907793 30. C. R. Shalizi, “ Optimal nonlinear prediction of random fields on networks ,” in Discrete Models for Complex Systems, DMCS (2003), pp. 11– 30 . system dynamics of agents' collective behavior, oscillation, or a phase transition. The proposed computational framework consists of two phases. First, we design a statistical basis for the analysis of agent behaviors. This step is loosely based on Shalizi and Crutchfield's ϵ–machine reconstruction that finds and groups common patterns of system behavior into what we will define as behavioral primitives. 6,7,28,30 Thermodynamic depth of causal states: Objective complexity via minimal representations ,” Phys. Rev. E 59, 275– 283 (1999). 6. J. P. Crutchfield and C. R. Shalizi, “,” Phys. Rev. E, 275–(1999). https://doi.org/10.1103/PhysRevE.59.275 Inferring statistical complexity ,” Phys. Rev. Lett. 63, 105– 108 (1989). 7. J. P. Crutchfield and K. Young, “,” Phys. Rev. Lett., 105–(1989). https://doi.org/10.1103/PhysRevLett.63.105 Computational mechanics: Pattern and prediction, structure and simplicity ,” J. Stat. Phys. 104(3–4), 817– 879 (2001). 28. C. Shalizi and J. Crutchfield, “,” J. Stat. Phys.(3–4), 817–(2001). https://doi.org/10.1023/A:1010388907793 30. C. R. Shalizi, “ Optimal nonlinear prediction of random fields on networks ,” in Discrete Models for Complex Systems, DMCS (2003), pp. 11– 30 . network of how agents' behavior changes over time. The network is used to visualize and analyze the landscape of the behavioral primitives of a model's dynamics. To understand how anarrives at its final configuration or how it solves a given problem, we definein anas the behavioral patterns of the system's agents, andas a collection of behaviors along with the transition probabilities of the system's agents changing from one behavior to another. Thecan then be viewed as a state-space transition graph of system's constituent behaviors that result in aachieving its goal.Here, a problem-solution has a loose definition that refers to even the simpleof agents' collective behavior, oscillation, or a phase transition. The proposed computational framework consists of two phases. First, we design a statistical basis for theof agent behaviors. This step is loosely based on Shalizi and Crutchfield's–machine reconstruction that finds and groups common patterns ofbehavior into what we will define as behavioral primitives.Next, we construct a state-space transitionof how agents' behavior changes over time. Theis used to visualize andthe landscape of the behavioral primitives of a

III. METHODOLOGY Section: Choose Top of page ABSTRACT I. INTRODUCTION II. BACKGROUND III. METHODOLOGY << IV. RESULTS AND ANALYSIS V. CONCLUSION, DISCUSSION... CITING ARTICLES Our methodology generalizes the behavior of an ABM and analyzes it as a network of behavioral state-space transitions (also referred to as a behavioral space) in four distinct steps: data collection, creation of a co-occurrence matrix from agent behaviors, clustering of behaviors into common behaviors, and the creation of a final transition matrix that describes the behavioral space. This methodology generalizes small scale, high frequency stochastic behavior while using the system-wide behaviors to describe information and its transfer among the model's agents. A. Data collection models, a randomly seeded Voronoi regions or a regular grid can be used to simplify the simulation space. Figure 1(a) First, the agents' simulation space is tessellated into discrete regions. Space tessellations can have regular shape of an orthogonal (x/y), honeycomb, triangular grids, or irregular patterns such as geo-socio-political boundaries, Voronoi tiling, or watershed boundaries. For spatially non-explicita randomly seeded Voronoi regions or a regular grid can be used to simplify the simulation space. Figureshows an example of a regular grid used to record agent histories. Next, the direction of each agent's movement is recorded, assigning values 1–4 to agent's northern, western, southern, and eastern crossing of grid boundaries, respectively. This spatially non-explicit simplification allows agents in different areas of the model with similar behavior to appear as having the same overall behavior pattern (translational invariance). Increasing the precision of the recorded histories will increase the fidelity of the analysis which can be achieved by using the tessellations with higher number of cardinal directions. The space division eliminates high frequency noise by ignoring any agent behavior that is smaller than the size of a grid cell and records only the low frequency behavioral patterns. B. Co-occurrence matrix C C in three steps. First, at each agent history step x t , we extracted a vector of agent's previous δ − and future δ + behaviors (Figures 1(a) 1(b) v ( δ − ) and v ( δ + ) v ( δ − ) = 〈 c o u n t ( 1 ∈ δ − ) , c o u n t ( 2 ∈ δ − ) , c o u n t ( 3 ∈ δ − ) , c o u n t ( 4 ∈ δ − ) 〉 v ( δ + ) = 〈 c o u n t ( 1 ∈ δ + ) , c o u n t ( 2 ∈ δ + ) , c o u n t ( 3 ∈ δ + ) , c o u n t ( 4 ∈ δ + ) 〉 , (1) and finally the number of previous v ( δ − ) and future v ( δ + ) four-tuple pairs is logged into the co-occurrence matrix C for all agents for all times x t + 1 . The agent's past and future behaviors are recorded from its histories into the co-occurrence matrixin three steps. First, at each agent history step, we extracted a vector of agent's previousand futurebehaviors (Figuresandshow the history vector length of range = 5). These vectors are then compressed into a past and future four-tupleandand finally the number of previousand futurefour-tuple pairs is logged into the co-occurrence matrixfor all agents for all times δ are compressed into a four-tuple. Each tuple element has the count of agent's directional moves in a given history vector δ. The example in Figures 1(a) 1(b) 2 , 1 , 2 , 1 , 2 , 1 is compressed as 〈 3 , 3 , 0 , 0 〉 , as would be the vector 1 , 2 , 1 , 2 , 1 , 2 . This compression reduces the effective dimension of the problem and provides invariance of recorded behaviors. In other words, regardless of where in time the record of agent's histories starts, the histories are compressed into the same tuple. The past and future history vectorsare compressed into a four-tuple. Each tuple element has the count of agent's directional moves in a given history vector. The example in Figuresandrecorded agent moves in four directions, thus the four-tuple compression of each history vector. For example, a history vector length range = 6, the agent's history vectoris compressed as, as would be the vector. This compression reduces the effective dimension of the problem and provides invariance of recorded behaviors. In other words, regardless of where in time the record of agent's histories starts, the histories are compressed into the same tuple. C C ( v ( δ x t − ) , v ( δ x t + ) ) = c o u n t ∀ k ( v ( δ x t + k − ) : v ( δ x t + k + ) ) (2) records each adjacent past and future history tuples δ x t + k − and δ x t + k + for all agents, for all sliding history vector windows x t (also see Figures 1(b) 1(c) C(i, j) represents how many times agent's past behavior i was followed by the future behavior j. The matrix row C ( i , : ) is then a frequency distribution (or a likelihood) of how often agent's past behavior i resulted in any of the future behaviors j. The co-occurrence matrixrecords each adjacent past and future history tuplesandfor all agents, for all sliding history vector windows(also see Figuresand). The matrix rows have the past history tuples, and the columns have the future history tuples, so the matrix cell) represents how many times agent's past behaviorwas followed by the future behavior. The matrix rowis then a frequency distribution (or a likelihood) of how often agent's past behaviorresulted in any of the future behaviors C. Clustering past behaviors into behavioral primitives ϵ are created by clustering the rows of co-occurrence matrix C into significant behaviors, also referred to as behavioral clusters. Each cluster ϵ i is defined as ϵ i = if ( χ 2 ( C ( i , : ) , C ( j , : ) ) < α ϵ ) then ϵ i ∪ { C ( j , : ) } . (3) The behavioral primitivesare created bythe rows of co-occurrence matrixinto significant behaviors, also referred to as behavioral clusters. Each clusteris defined as clustering implementation first selects a random row C ( i , : ) that has not been assigned to a behavioral cluster to represent a new leader or representative of the behavioral cluster – ϵ i . A selected leader is then compared with all other matrix rows C ( j , : ) that have not been assigned to a behavioral cluster. The row C ( j , : ) is added to the behavioral cluster ϵ i represented by a leader C ( i , : ) if a two sample hypothesis test χ 2 χ 2 ( C ( i , : ) , C ( j , : ) ) = ∑ k = 1 | C ( i , : ) | ( K 1 C ( i , k ) − K 2 C ( j , k ) ) 2 C ( i , k ) + C ( j , k ) , where K 1 = ∑ k = 1 | C ( j , : ) | ( C ( j , k ) ) ∑ k = 1 | C ( i , : ) | ( C ( i , k ) ) K 2 = ∑ k = 1 | C ( i , : ) | ( C ( i , k ) ) ∑ k = 1 | C ( j , : ) | ( C ( j , k ) ) (4) is accepted. Theimplementation first selects a random rowthat has not been assigned to a behavioral cluster to represent a newof the behavioral cluster –. A selected leader is then compared with all other matrix rowsthat have not been assigned to a behavioral cluster. The rowis added to the behavioral clusterrepresented by a leaderif a two sample hypothesis testis accepted. 1(c) clustering of matrix rows into behavioral clusters ϵ i . The process of picking a leader from unassigned rows, finding all similar rows, and adding them to the cluster represented by a leader is repeated until all rows have been assigned to a behavioral cluster. Figure(right) shows possibleof matrix rows into behavioral clusters χ 2 to calculate the similarity between two multidimensional row vectors in the co-occurrence matrix C. 29 Automatic filters for the detection of coherent structure in spatiotemporal systems ,” Phys. Rev. E 73, 036104 (2006). 29. C. R. Shalizi, R. Haslinger, J.-B. Rouquier, K. L. Klinkner, and C. Moore, “,” Phys. Rev. E, 036104 (2006). https://doi.org/10.1103/PhysRevE.73.036104 χ 2 tests the statistical independence between two past behaviors if they are sufficiently close to be generalized as one prototypical behavioral primitive. Basically, χ 2 performs a hypothesis test if two past behaviors are statistically equivalent in terms of their future behaviors. Semantically, if a row is similar to another row, it means that those agents' past behaviors have a similar enough distribution of future behaviors to be considered the same behavioral primitive. We usedto calculate the similarity between two multidimensional row vectors in the co-occurrence matrixtests the statistical independence between two past behaviors if they are sufficiently close to be generalized as one prototypical behavioral primitive. Basically,performs a hypothesis test if two past behaviors are statistically equivalent in terms of their future behaviors. D. The state-space transition matrix T records how likely agents are to change their behavior from one behavioral primitive to another in two consecutive time-steps. The square matrix T has a size of the number of behavioral primitives | ϵ | × | ϵ | . The recorded history of agent movement (described in Section T. An agent's past history vector δ − is extracted, and encoded into a four-tuple v ( δ − ) , and then, the four-tuple's membership in a behavioral cluster ϵ i is looked up (Figure 1(d) T T ( p , r ) = count ∀ k ( ϵ p : ( v ( δ x t − ) ∈ ϵ p ) , ϵ r : ( v ( δ x t + 1 − ) ∈ ϵ r ) ) (5) records how many times an agent changed its behavior from one behavioral cluster to another in two subsequent time steps x t and x t + 1 (also see Figure 1(e) δ x t − and δ x t + 1 − belong to the same behavioral cluster ϵ p = ϵ r , then the agent's behavior remains unchanged. If the agent's compressed history tuples belongs to two different behavioral clusters, then the agent's behavior transitioned from one pattern of behavior to another ( ϵ p ≠ ϵ r ). A transition matrixrecords how likely agents are to change their behavior from one behavioral primitive to another in two consecutive time-steps. The square matrixhas a size of the number of behavioral primitives. The recorded history of agent movement (described in Section III A ) is used for the second time to create the matrix. An agent's past history vectoris extracted, and encoded into a four-tuple, and then, the four-tuple's membership in a behavioral clusteris looked up (Figure). The transition matrixrecords how many times an agent changed its behavior from one behavioral cluster to another in two subsequent time stepsand(also see Figure). If the agent's past vectorsandbelong to the same behavioral cluster, then the agent's behavior remains unchanged. If the agent's compressed history tuples belongs to two different behavioral clusters, then the agent's behavior transitioned from one pattern of behavior to another (). ϵ × ϵ square transition matrix T preserves and further compresses the co-occurrence matrix C (see Section model's behavior using the statistical analysis (see Section network's nodes) is a constant and only dependent on the agent's attraction range (see Thesquare transition matrixpreserves and further compresses the co-occurrence matrix(see Section III B ) by generalizing thebehavior using the(see Section III C ). The upper-bound for the number of behavioral primitives (thenodes) is a constant and only dependent on the agent's attraction range (see Appendix A ).

IV. RESULTS AND ANALYSIS Section: Choose Top of page ABSTRACT I. INTRODUCTION II. BACKGROUND III. METHODOLOGY IV. RESULTS AND ANALYSIS << V. CONCLUSION, DISCUSSION... CITING ARTICLES model of Newton's universal gravitation law, where the force F acting on two bodies with the mass m 1 , m 2 with distance between objects d is F = G * ( m 1 * m 2 ) / ( d 2 ) . We implemented the model using two agent types: stationary gravitational attractors and moving agents that were drawn towards the attractors with variable strength. Instead of using multiple parameters of the Newton's equation, the model implements a single parameter a t t r a c t i o n = [ 0 , 7 ] for the attractors that acts as a gravitational force on the moving agents. No interactions among moving agents were assumed. We chose this model because varying single attraction parameter results in random, complex, and stable oscillating system regimes shown in Figures 2(aa) 2(cc) 2(dd) model snapshots). We applied the proposed methodology on aof Newton's universal gravitation law, where the forceacting on two bodies with the masswith distance between objectsis. We implemented theusing two agent types: stationary gravitationaland movingthat were drawn towards thewith variable strength. Instead of using multiple parameters of the Newton's equation, theimplements a single parameterfor thethat acts as a gravitational force on the moving agents. No interactions among moving agents were assumed. We chose thisbecause varying singleparameter results in random, complex, and stable oscillatingregimes shown in Figures, and, respectively (see Appendix B for additionalsnapshots). model's toroidal simulation area was tessellated into an orthogonal 50 × 50 grid with the cell size of 50 space units. The agents had constant velocity of 1 space unit per time-step. The model was tested with the attraction strength in the range of [ 0 , 7 ] . Agent trajectories were recorded after the initial 200 time-steps to ensure a stable system behavior, and the model was executed with 150 agents for 5500 time-steps. The analysis used the agents' log of ≈ 500 000 trajectory data-points and the past and future history vectors had length range = 15 with the number of cardinal direction being set to 4 – {N, W, S, E}. The size of co-occurrence matrix C was 816 × 816 past and future histories for a given range (see χ 2 clustered co-occurrence matrix rows into the ϵ behavioral primitives for α ϵ > 0.9 . This is rather strict criteria to reject the null hypothesis, but it very well illustrates the clustering semantics and the emergence of behavioral primitives. The state-space network is constructed from matrix T and simplified by using the edge weight > 0.0004, the nodes with degree > 0, the node's size is proportional to the sum of its weighted in-degree and its weighted out-degree, and the network nodes are grouped by common color into the clusters with high mutual node connectivity. 2,16 2. V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre, “ Fast unfolding of communities in large networks ,” J. Stat. Mech.: Theory Exp. P10008, 1– 12 (2008). Laplacian dynamics and multiscale modular structure in networks ,” {arXiv preprint 16. R. Lambiotte, J.-C. Delvenne, and M. Barahona, “,” {arXiv preprint arXiv:0812.1770 } (2009). Thetoroidal simulation area was tessellated into an orthogonal 50 × 50 grid with the cell size of 50 space units. The agents had constant velocity of 1 space unit per time-step. Thewas tested with the attraction strength in the range of. Agent trajectories were recorded after the initial 200 time-steps to ensure a stablebehavior, and thewas executed with 150 agents for 5500 time-steps. Theused the agents' log oftrajectory data-points and the past and future history vectors had length= 15 with the number of cardinal direction being set to 4 – {N, W, S, E}. The size of co-occurrence matrixwas 816 × 816 past and future histories for a given range (see Appendix A for details). The hypothesis testclustered co-occurrence matrix rows into thebehavioral primitives for. This is rather strict criteria to reject the null hypothesis, but it very well illustrates thesemantics and the emergence of behavioral primitives. The state-spaceis constructed from matrixand simplified by using the edge> 0.0004, the nodes with> 0, the node's size is proportional to the sum of its weighted in-degree and its weighted out-degree, and thenodes are grouped by common color into the clusters with high mutual node connectivity. A. Analysis a t t r a c t i o n = [ 0 , 1.0 ] with agent behaviors in pseudo-straight trajectories as their behavior is not affected by the attractors' gravitational pull. Figure 2(aa) attraction = 0.0, and the resulting trajectories marked as Features i that can be informally described as “ ← , → , ↑ , ↓ , ↗ , ↘ , ↙ , ↖ ” lines. The corresponding state-space network (a) shows the results of how the analysis generalized agent behaviors using recorded trajectories. The behaviors highlighted in the model's snapshots (aa) are marked as corresponding group of network nodes in the resulting state-space network (a). Each group of nodes represents closely related agent behaviors with the recorded trajectories deviating only a few degrees. For example, agent trajectory i ↖ corresponds to a sub-graph with behavioral primitives with high values in the first two positions corresponding to moves in the N and W directions: 〈 10 , 5 , 0 , 0 〉 , 〈 11 , 4 , 0 , 0 〉 , 〈 9 , 6 , 0 , 0 〉 , 〈 8 , 7 , 0 , 0 〉 , and 〈 7 , 8 , 0 , 0 〉 connected with low weight edges. Recall that each network node represents one ϵ-cluster (or behavioral primitive), and the weighted edges are the likelihood of an agent with a given behavior transitioning to a different behavioral primitive in the next times-steps. The agents' behavior in the random regime is unlikely to change which is shown as low-weight edges connecting the network nodes and mutually disconnected node clusters. We can also say, that the behavior likely stay unchanged as seen by node's high-weight self-leading edge. The decoupled nature between the attractors and agents resulted in the mutual dis-connectivity among groups of nodes representing similar behaviors in the state-space network. The random regime occurs forwith agent behaviors in pseudo-straight trajectories as their behavior is not affected by thegravitational pull. Figureshows the agent behavior for= 0.0, and the resulting trajectories marked as Featuresthat can be informally described as “” lines. The corresponding state-space(a) shows the results of how thegeneralized agent behaviors using recorded trajectories. The behaviors highlighted in thesnapshots (aa) are marked as corresponding group ofnodes in the resulting state-space(a). Each group of nodes represents closely related agent behaviors with the recorded trajectories deviating only a few degrees. For example, agent trajectorycorresponds to a sub-graph with behavioral primitives with high values in the first two positions corresponding to moves in the N and W directions:, andconnected with low weight edges. Recall that eachnode represents one-cluster (or behavioral primitive), and the weighted edges are the likelihood of an agent with a given behavior transitioning to a different behavioral primitive in the next times-steps. The agents' behavior in the random regime is unlikely to change which is shown as low-weight edges connecting thenodes and mutually disconnected node clusters. We can also say, that the behavior likely stay unchanged as seen by node's high-weight self-leading edge. The decoupled nature between theand agents resulted in the mutual dis-connectivity among groups of nodes representing similar behaviors in the state-space attractors. Figures 2(b) 2(bb) model configuration with attraction = 1.5 that represents a borderline between the random and complex regimes. The resulting agent trajectories are a mix of pseudo-linear and various “L, ⊂ , ⊃ , ∩ , ∪ ” shapes. Unlike the random regime, all identified behaviors are connected showing that the agents can and did transition from one stable behavior to another. For example, mutually connected groups of nodes that correspond to Features i ↖ and i ↑ are now connected forming a group of nodes marked as a feature ii that corresponds to a “L-shape” trajectory. A path through a state-space network describes a non-trivial agent behavior as well as defined the behavioral precursors that led to agent's current behavior. As the value of the attraction parameter increases, the agent behavior is increasingly altered by theFiguresandshow theconfiguration with= 1.5 that represents a borderline between the random and complex regimes. The resulting agent trajectories are a mix of pseudo-linear and various “” shapes. Unlike the random regime, all identified behaviors are connected showing that the agents can and did transition from one stable behavior to another. For example, mutually connected groups of nodes that correspond to Featuresand↑ are now connected forming a group of nodes marked as a featurethat corresponds to a “L-shape” trajectory. A path through a state-spacedescribes a non-trivial agent behavior as well as defined the behavioral precursors that led to agent's current behavior. a t t r a c t i o n = [ 2.0 , 5.0 ] as the coupled dynamics between the attractors and the agents increases the structured behavior. Figures 2(c) 2(cc) model configuration with attraction = 3.0 with the model behavior of all behavioral Features i-vi. The agents are intermittently captured and then released by the attractors with their trajectories significantly altered from their original movement, resulting in linear, pseudo-linear, hyperbolic, partially orbiting, S, U, and ∞ behavior shapes. The agents likely transition between similar behaviors, which resulted in a 63% increase in the number of edges with the weight greater than 200 compared to the random regime networks. The prevalent hyperbolic trajectories, marked as Features iv, are now represented by a single network node that represents trajectories 〈 5 , 9 , 1 , 0 〉 , 〈 0 , 2 , 9 , 4 〉 , 〈 0 , 3 , 7 , 5 〉 , 〈 0 , 7 , 8 , 0 〉 , 〈 2 , 8 , 5 , 0 〉 , 〈 3 , 0 , 5 , 7 〉 , … , 〈 7 , 3 , 0 , 5 〉 , 〈 8 , 3 , 0 , 4 〉 , 〈 8 , 4 , 0 , 3 〉 , 〈 8 , 5 , 0 , 2 〉 , 〈 9 , 5 , 0 , 1 〉 , 〈 10 , 2 , 0 , 3 〉 (6) as one behavioral primitive. For example, a hyperbolic trajectory ⊂ is abstracted as transition from 〈 0 , 7 , 8 , 0 〉 → 〈 0 , 3 , 7 , 5 〉 → 〈 3 , 0 , 5 , 7 〉 of agent behaviors that are now clustered into the same behavioral primitive. The clustering of the past behaviors with statistically similar future behavior is symmetric, so the same sequence of three transitions also represents ⊃ trajectory. The complex regime occurs inas the coupledbetween theand the agents increases the structured behavior. Figuresandshow theconfiguration with= 3.0 with thebehavior of all behavioral Features. The agents are intermittently captured and then released by thewith their trajectories significantly altered from their original movement, resulting in linear, pseudo-linear, hyperbolic, partially orbiting, S, U, andbehavior shapes. The agents likely transition between similar behaviors, which resulted in a 63% increase in the number of edges with the weight greater than 200 compared to the random regimeThe prevalent hyperbolic trajectories, marked as Features, are now represented by a singlenode that represents trajectoriesas one behavioral primitive. For example, a hyperbolic trajectoryis abstracted as transition fromof agent behaviors that are now clustered into the same behavioral primitive. Theof the past behaviors with statistically similar future behavior is symmetric, so the same sequence of three transitions also representstrajectory. Features v. and vi. are the other frequent, non-trivial agent trajectories that are represented as single network nodes but abstract multiple behaviors. The agents are frequently attracted into oscillating, ellipsoid behaviors that they can revert out. The analysis abstracted these ellipsoid trajectories as one behavioral primitive. The node representing an ellipsoid behavior is marked as feature vi and contains 〈 4 , 5 , 5 , 1 〉 , 〈 3 , 4 , 4 , 4 〉 , 〈 3 , 5 , 5 , 2 〉 . Feature v marks other infrequent, transitional, complex behaviors found in this regime. a t t r a c t i o n ≥ 5.0 . Figure 2(dd) model configuration for attraction = 7.0 with agents captured by an attractor and unable to escape. The resulting agent trajectories were various “ ◯ ” shapes (Feature vi) which created a dynamic but stable model behavior. Feature vii shows agents transitioning from one stable oscillating behavior to another with a similar shape. Although the resulting transition network (d) appears to be complicated, the overall system behavior can be described as nine stable ellipsoid behaviors showed in nine distinct node colors. Each color cluster of ellipsoid behavioral primitives represents one general orientation of ellipsoid trajectories. The stable oscillating regime is found for the. Figureshows theconfiguration for= 7.0 with agents captured by anand unable to escape. The resulting agent trajectories were various “” shapes (Feature) which created abut stablebehavior. Featureshows agents transitioning from one stable oscillating behavior to another with a similar shape. Although the resulting transition(d) appears to be complicated, the overallbehavior can be described as nine stable ellipsoid behaviors showed in nine distinct node colors. Each color cluster of ellipsoid behavioral primitives represents one general orientation of ellipsoid trajectories. B. Transitional map similarity between ABMs model behaviors as transition networks of behavioral primitives and weighted transition edges for stochastic agent based models, we re-initialized the model with random seed of each attraction configuration and recorded the agent trajectories for 20 independent runs of the model for attraction in the range of 0 to 7 in half steps, resulting in 300 independent data sets to analyze. Figure 3(a) analysis for each attraction value. The error bars for each data point show the variance in the number of behavioral primitives found with a minimum and maximum variance of 0% (at attraction = 0.0) and 3.5% (at attraction = 0.5), respectively. The maximum variance is when the system dynamics start transitioning between two regimes such as the random and complex behavior, or from the complex behavior to the stable oscillating. The number of behavioral primitives for each attraction value in Figure 3(a) α and after certain execution time no new behavioral primitives were found. To assess the consistency of the above methodology to abstract commonbehaviors as transitionof behavioral primitives and weighted transition edges forwe re-initialized thewith random seed of each attraction configuration and recorded the agent trajectories for 20 independent runs of thefor attraction in the range of 0 to 7 in half steps, resulting in 300 independent data sets toFigureshows the number of behavioral primitives generalized by thefor each attraction value. The error bars for each data point show the variance in the number of behavioral primitives found with a minimum and maximum variance of 0% (at= 0.0) and 3.5% (at= 0.5), respectively. The maximum variance is when thestart transitioning between two regimes such as the random and complex behavior, or from the complex behavior to the stable oscillating. The number of behavioral primitives for each attraction value in Figureasymptotically converged to the reported (y-axis) value as a function of time. This attests to the methodology's stability to generalize behavior, since for a chosenand after certain execution time no new behavioral primitives were found. T corresponding to two independent model executions. The results in Figure 3(b) ABM executions. If the transition matrix T captures the shape of the behavioral landscape of system's dynamics, then the scalar MSE value is an indicator on how close the geometries are in the two behavioral landscapes. The geometry of the behavioral spaces is close to identical for the attraction = 0.0 (random regime), self-similarity for each attraction configuration (the diagonal MSE scores), and for the configurations with a t t r a c t i o n ≥ 5.5 when the system dynamics are described as a network of stable oscillating behavioral primitives. The high self-similarity values at these attraction values are indicative of the behavioral similarity as well as the methodology's consistency to generalize the system's tell-tale behaviors. We also use the Mean Square Error (MSE) test to measure the distance between the result matricescorresponding to two independentexecutions. The results in Figureshow the pairwise comparisons of the self-similarity for theexecutions. If the transition matrixcaptures the shape of the behavioral landscape of system'sthen the scalar MSE value is an indicator on how close the geometries are in the two behavioral landscapes. The geometry of the behavioral spaces is close to identical for the= 0.0 (random regime), self-similarity for each attraction configuration (the diagonal MSE scores), and for the configurations withwhen theare described as aof stable oscillating behavioral primitives. The high self-similarity values at these attraction values are indicative of the behavioral similarity as well as the methodology's consistency to generalize the system's tell-tale behaviors. The increasing structure in the random regime can be seen as the dark gray self-similarity for a t t r a c t i o n = [ 0.5 , 1.5 ] . Model configurations for a t t r a c t i o n = [ 1.5 , 5.5 ] have emergent, complex system dynamics with diverse geometry of the corresponding state-space transition networks with low MSE values. Finally, the convergence of the system's dynamics into a stable oscillating behaviors is separated by the dark gray a gradient for the MSE values of a t t r a c t i o n = [ 5.0 , 6.0 ] with the stable oscillating behaviors of matrices for a t t r a c t i o n = [ 6.0 , 7.0 ] . The MSE measure not only identifies if a system behavior has changed with changes to the model's initial conditions but also gives a qualitative comparison of the overall systems dynamics. The MSE heat map plot supports the classification of a system's dynamics into random, complex, and stable oscillating regimes. The regime changes are also visible as a darker band of dis-similar MSE values at the attraction = 1.5 and convergence to stable oscillating behavior as the gray-scale gradient at a t t r a c t i o n = [ 3.0 , 5.0 ] .

V. CONCLUSION, DISCUSSION, AND FUTURE WORK Section: Choose Top of page ABSTRACT I. INTRODUCTION II. BACKGROUND III. METHODOLOGY IV. RESULTS AND ANALYSIS V. CONCLUSION, DISCUSSION... << CITING ARTICLES We present a computational framework to analyze the system behavior for any stochastic ABM. The Geometry of Behavioral Spaces generalizes common not dependent on in-depth domain knowledge of the analyzed system or the nature of agents' interactions. The Geometry of Behavioral Spaces is not dependent on in-depth domain knowledge of the analyzed system or the nature of agents' interactions. Our methodology generalizes common behavioral patterns of complex systems dynamics and creates a state-space transition network of the agent's behaviors. This network view represents the behavioral landscape of an ABM and can be used as a meta-heuristic for a fitness function of stochastic search algorithms, an analytical tool for optimizing the model's parameter configuration space, or identifying the behavioral building blocks in the ABM's system dynamics. Our proposed methodology is empirical in its nature and should not be mistaken for a theoretical model or a proof of computation. network offers a quantitative and information rich view of spatially explicit ABM behaviors using the discrete temporal and spatial analyses presented in Sec. network topology can be used to measure the distance between two behavioral primitives, calculate the diversity of the exhibited behaviors in the ABM, or analyze the structural holes in the behavioral landscape. The state-spaceoffers a quantitative and information rich view of spatially explicitbehaviors using the discrete temporal and spatialpresented in Sec. IV . Thetopology can be used to measure the distance between two behavioral primitives, calculate the diversity of the exhibited behaviors in theorthe structural holes in the behavioral landscape. To date, the common approaches to evaluate ABM dynamics in search of emergent properties have either been through a binary decision or subjective human analysis of system dynamics. Our proposed methodology will explain analytically how the final model configuration was reached, assess the model's parameter sensitivity, and allow analysis to be scaled to a large number of model executions. The explanation of a model's system dynamics is performed by constructing a transition network of generalized behavioral primitives found in the agent's behaviors. Emergent system dynamics can be represented as network of structures that correspond to the introduction of structured spatio-temporal patterns of behavior. Our methodology's unique utility is to identify different system behaviors, outline the transitions between different system regimes, and consistently generalize stochastic behavior were illustrated on an ABM model with well understood yet complex emergent behaviors. Although we tested the proposed methodology on one well-understood ABM, its system dynamics cover a wide range of behaviors which gives us confidence that our proposed methodology will scale to a wide range of potential applications. Our future work will focus on testing our proposed methodology on models with asynchronous execution, agents with non-uniform velocities, and simulations with strong and weak emergence. We will investigate the trade-off between the lossy compression of agent trajectories and the scalability and the framework's ability to detect complex attractors. As an engineering application, we propose using the proposed methodology in the evolutionary algorithms to evolve high diversity model configurations that satisfy a desired simulation goal, and classifying behavior of bacterial colonies in system biology.

ACKNOWLEDGMENTS This work would not be possible without generous funding from Alaska NSF EPSCoR Grant No. 1208927. The authors would like to thank all reviewers and project collaborators for their comments and suggestions.