Table of Contents Author Guidelines Submit a Manuscript

A corrigendum for this article has been published. To view the corrigendum, please click here.

Volume 2017 (2017), Article ID 3548591, 16 pages
Research Article

Self-Organized Societies: On the Sakoda Model of Social Interactions

1Departamento de Ingeniería Industrial, Universidad de los Andes, Bogotá, Colombia
2CeiBA Complex Research Center, Bogotá, Colombia
3Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Avda. Diagonal las Torres 2640, Peñalolén, Santiago, Chile
4UAI Physics Center, Universidad Adolfo Ibáñez, Santiago, Chile

Correspondence should be addressed to Sergio Rica

Received 3 October 2016; Revised 7 December 2016; Accepted 14 December 2016; Published 23 January 2017

Academic Editor: Mattia Frasca

Copyright © 2017 Pablo Medina et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


We characterize the behavior and the social structures appearing from a model of general social interaction proposed by Sakoda. The model consists of two interacting populations in a two-dimensional periodic lattice with empty sites. It contemplates a set of simple rules that combine attitudes, ranges of interactions, and movement decisions. We analyze the evolution of the 45 different interaction rules via a Potts-like energy function which drives the system irreversibly to an equilibrium or a steady state. We discuss the robustness of the social structures, dynamical behaviors, and the existence of spatial long range order in terms of the social interactions and the equilibrium energy.

1. Introduction

Throughout the history, the society has self-organized in a myriad of social structures and behaviors which appears as a response to the attitudes, decisions, and expectations among individuals, that is to say, from local (usually simple) rules to global behavior [1]. This capacity appears in animals [2], insects [3, 4], among others, and it flourishes in human beings since the early life [5]. The basic underlying mechanism of this self-organization lies on the evaluation of the expectations resulting from the attitudes among individuals inducing, or not, a mobility decision.

Interestingly, despite the complex nature of social interactions, the social behavior shares many common features with a variety of physical systems. Segregation, inclusion, and aggregation are examples of a collective order arising from simple and local individual rules based on attitudes, decisions, and expectations among individuals [69].

In the late forties, Sakoda introduced in his Ph.D. thesis [10] a first general discrete dynamical model for social interactions and published later [11], at the same time of the well-known Schelling’s social segregation model [12, 13], which is, as we will show, only a special class already included in the original Sakoda model. The particular case of the Schelling segregation model and its variations have been studied widely through agent based models [9, 1423]. However, the Sakoda dynamics and its underlying richness, which explain other different social phenomena far from segregation, are mostly unknown; thus they have not been, to the best of our knowledge, studied previously in a deeper manner.

The original Sakoda model roughly consists of social interactions among two groups of individuals evolving in a network according to specific attitudes of attraction, repulsion, and neutrality. An individual evaluates its social expectative at all possible available locations, preferring those near individuals associated with attractive (positive) attitudes and avoiding locations near individuals associated with repulsive (negative) ones. This procedure is repeated randomly among all possible individuals; henceforth Sakoda’s algorithm is iterated recursively driving the system, under some conditions, to a well-organized spatial pattern.

The dynamics is quite rich because of three aspects: () the large number of combinations of the possible attitudes, () the effect of the separation distance among individuals (i.e., the interaction could be of short or long range), and () the mobility of the individuals (the individuals may move into its own neighborhood or they may migrate far away).

This produces a myriad of possible patterns of the individuals’ arrangements in the space (social structures) to be analyzed. Some of them were already recognized by Sakoda in his early work [11].

In the present paper we discuss and complement the Sakoda original classification of social structures. More precisely, we characterize all independent individual interactions (45 different interactions) in the case of short range interaction and long range movement. We introduce a Potts energy [24, 25], which has been the natural -state generalization of the Ising model. As the former, the Potts model also displays a phase transition in statistical physics. In the context of the Sakoda model, the Potts energy happens to be a Lyapunov function in the case of symmetric interactions that governs the future dynamics and quantifies different attractors in terms of the parameters of the problem.

This Potts-like energy turns out to be a common principle for all possible interactions and describes the evolution of different social interacting systems. Finally, we discuss the role of empty sites and the existence of long range order in the social structures.

The present paper is organized as follows. Section 2 presents the Sakoda model and the numerical scheme. Section 3.1 shows the energy principle that mandates the dynamical evolution of the system. Next, Section 4 characterizes the various social structures arising from different interactions possibilities. A quantitative energy-based study is done in Section 5.1. The role of vacancies is studied in Section 5.2 and the existence of long range order is discussed in Section 5.3. Finally, we conclude in Section 6 and complement the analytics with a series of appendices.

2. The Basic Model

The Sakoda model consists in a regular and periodic lattice with nodes in a two-dimensional space. Each node is associated with a discrete variable and may take the values . States denote an occupancy by an individual of one of the two types (+ or −); the state denotes an empty place or a vacancy. Only when a node is empty may an individual (+ or −) occupy the free spot and reassign its value with its current one.

The social interaction is mainly characterized by the possible attitudes among individuals which are summarized in a matrix. We call hereafter the “-matrix,” which for the sake of space we will symbolize in a vectorized form: . The entries take three possible values: . These indicate, respectively, an attractive , a neutral (0), or a repulsive attitude from members of the type to the individuals of the type . There exist different possible cases which are reduced, because of symmetry, to 45 independent social structures. In Appendix A, we discuss the symmetries of the “-matrices” and its properties in a deeper manner.

Next, we model the preferences of one individual over a particular place in the lattice based on the spatial locations and the attitudes of other individuals. Following the original ideas of Sakoda, we propose a function that quantifies the social expectations of the individual :where denotes a symmetric () interaction strength. We avoid self-interactions taking . On the other hand, the sign of the interaction is given by (do not confuse the variable that takes the values with the entries of the matrix that take the same values ):

In general, this expression is not symmetric; that is, . We underline that though is symmetric, the full interaction is not necessarily symmetric. The coefficients may include short and long range interactions. The mobility of an individual could be also of long or short movement. An individual located at the node would move towards an empty node , if . Here represents the potential value (1) evaluated at the empty node, , occupied by the individual . As a general rule, the individual chooses the place everywhere in the lattice that makes the social expectation function the highest. In case of degeneracy, that is, if there are nodes, , such that , one of these is selected randomly with the same probability. In a short movement case, the mobility of the individual is restricted to its Moore neighborhood (the eight closest neighbors).

To summarize, the movement of individuals occurs in the following manner: during an iteration step, an individual is selected randomly. Then, it evaluates the potential function at every available empty site in its range of movement for a potential new location. Next, the individual selects the site with the highest value of the potential expectation (1) and then it moves to this site. If there is not a site that improves the potential function, then the individual stays in place. The algorithm is iterated until the system reaches a fixed point; otherwise, it runs forever.

We point out that the rule preserves the initial number of individuals of all types. That is, if , , and are the number of individuals of the types +, −, and the number of vacancies, respectively, they keep their initial value during the evolution. We define the fractions by and (naturally ). Hereafter, as in Sakoda original paper [11], we limit our study to the case of .

2.1. Numerics

For the following, we restrict our central results only to the case of short range interactions and long range movements: an individual evaluates its preferences considering the Moore neighborhood; that is, , if is in the 8 closest neighbors of the site , and elsewhere. Then, the individual may look forward the highest expectation at any empty node in the lattice. We discuss the cases of short and long range interactions/mobility at the end of the paper and in Appendix F.

The numerical simulations for all 45 possible -matrices were run in a two-dimensional periodic lattice of . The evolution usually runs for time steps, which happens to be sufficient to reach equilibrium or a steady state. We have also simulated the case of without significant dependence on the system size.

The initial condition consists of a random distribution of states characterized by a fraction of vacancy states. For all of our simulations, we have checked the robustness of the social structures considering different random initial configurations. We focus our attention on five levels () of different initial concentration of vacancies as a good representative sample (see Section 5.2). Other initial distributions may be studied using the same frame and tools developed in the present paper.

Figure 1 summarizes some social structures resulting from the afore-mentioned algorithm for four characteristic -matrices and for three different values of the vacancy fraction, . Over all the manuscript, the light gray cells (red cells online) represent individuals at the state and the dark gray cells (blue cells online) represent individuals, while the white ones represent the empty spaces, .

Figure 1: (Color online) snapshots of the few basic social structures: (a) , (b) , (c) , and (d) for various density of vacancies: (I) , (II) , and (III) . The simulations consider the case of short range of influence and long range of movement and were run in a lattice of sites and for steps.

From the simulations, we notice the similarities of our results with the physics of surface tension [26]. Indeed all possible situations among the three phases follow after an interface equilibria argument (something that will be clear in Section 3). For instance, Figure 1(a) shows a situation of a miscible and phases, but both are nonmiscible with respect to white. Similarly, Figure 1(b) presents a miscible phase between and 0 (white), but both are nonmiscible with phase. Figures 1(c) and 1(d) present cases with the three nonmiscible phases , but the interface equilibria are different. While in Figure 1(c) there are interfaces between /, , and ; in the case of Figure 1(d), / interface does not exist and it is replaced by a monolayer of white vacancies (excepting the case of a very rare number of vacancies, see Figure 1(d)(I) for ). For a given interaction, a robust behavior is clearly noticed for different vacancy fractions, but the patterns may differ substantially from one case to another because of the role of the empty sites. We shall discuss more in detail this point later on in Section 5.2.

Finally, we have checked the robustness of the phenomena for the case of a nondeterministic Sakoda instance, through considering the rule: An individual located at node would move towards an empty node , with a probability . Here . As expected, in the large limit one recovers the deterministic rule; on the contrary, in the limit the individuals move randomly to any available empty site. In this case the system behavior is highly fluctuating and it does not self-organize into any spatial structure. Therefore, it is expected that, for a given -matrix, a critical parameter exists, which controls the existence or not of a social structure.

3. Potts Energy for the Sakoda Model

3.1. An Energy Theorem for the Dynamics of the Sakoda Model

The evolution of the Sakoda model is ruled by an energy principle that we formulate in this section. After [27], the following theorem is shown.

Theorem 1. If the -matrix is symmetric, then the Potts-like energy function [24, 25], of a given configuration of the network ,is a nonincreasing function after a movement resulting from the Sakoda algorithm; that is, . In other words, is a nonstrict Lyapunov function.

Proof. Suppose that is a symmetric matrix. Thus, for two arbitrary sites and in the lattice. Consider the exchange of the sites (occupied) and (empty), that is, and . Hence, the initial and the subsequent configuration after a swap between the and sites are and If the -matrix is symmetric, energy (3) difference between the initial and the actual configuration becomesHere we have used that the self-interaction, , vanishes and for all . This occurs because the -function vanishes exactly if it is evaluated at a vacancy state (for details see Appendix B). The intrinsic mechanism of Sakoda’s rule imposes that the social expectation (1), , of the individual , is smaller than the social expectation at the site , that is, :Rewriting the r.h.s. of (4) in terms of the expressions and of (1) and (5), respectively, we obtain We emphasize that the equality holds in the case of ; in this case no exchange is carried out. In conclusion, if is a symmetric matrix, then ; that is, energy (3) is a nonincreasing quantity after a movement, concluding the proof. We emphasize that the proof is completely general; hence, it does not depend neither on the range of the interactions and the range of motion.

If the -matrix is not symmetric, energy (3) is not necessarily a Lyapunov function; thus the energy may increase or decrease indistinctly. However, as we shall see next, there is an energy decreasing principle for both cases, symmetric and nonsymmetric interactions. This principle is understood in a statistical sense: for a random configuration, in average, it is more probable to decrease than to increase the energy after a movement [23].

This principle is understood in a statistical sense: for a random configuration, in average, it is more probable to decrease than to increase the energy after a movement [23]. The observation of an energy decreasing principle for nonsymmetric interactions may indicate behind this principle a “free energy.” In [21] we concluded that this free energy is , where and are the Lagrange multipliers associated with the constant number of individuals of each type and the constant number of vacancies. Hence, the macrobehavior minimizes energy with the corresponding restrictions (see next section).

In Appendix C, we show that the energy (3) corresponds, in the present case, exactly to the spin 1 Ising model [20]. Moreover, if the -matrix is , one readily gets , which is exactly the expression of the Ising energy associated with Schelling’s model [21].

Finally, notice that the energy (3) may be written in terms of geometrical parameters [16, 20, 21]. For instance, in the case considered here of short range interaction, the energy readswhere is an effective interface length which is directly related to the number of edges, , and corners, , among the interface created by the states and (see Appendix D for a derivation of (6)). Finally, the term represents the bulk energy contribution. However, because are constant during the evolution, the bulk energy term does not change; therefore variations of the energy are produced only by interface length variations, if these exist. As we will see later, the parameter plays an important role in the social structure characterization.

3.2. Energy Behavior in the Sakoda Model

To illustrate the energy principle, Figure 2(i) plots the temporal evolution of the energy (3) for the same cases of Figure 1. The selection includes both symmetric and nonsymmetric interactions. It is easy to check that this energy function is a Lyapunov function for the former cases, while for the later it is not (see Figure 2(i)). Nevertheless, for all nonsymmetric interactions, the energy tends to a certain value indicating the existence of an equilibrium state (see curves (a) and (d) in Figure 2). This case suggests the existence of some kind of “statistical attractor,” something already noticed in Schelling’s model [23]. Despite the complexity involved in the Sakoda model, the long-time dynamics reaches a statistical attractor with a well-defined average energy and social structure (see Section 5.1).

Figure 2: (Color online) temporal evolution of the energy (i) and the effective perimeters (ii) for the same cases as in Figure 1 row (II), that is, for . The curves are labeled with an (a) for the snapshot (a II) of Figure 1, with a (b) for (b II) and so forth. (i) In this plot one notices clearly energy decreasing tendency in time. Here the curves (b) and (c) correspond to a nonsymmetric -matrix; hence its energy may decrease or increase indistinctly. Some increments of energy are labeled with . (ii) The perimeters are explicitly labeled in the plot with , , or . It is observed that the qualitative evolution agrees with the energy behavior. These special cases ran for time steps.

Quantitatively, the energy expansion (6) allows us to distinguish between the segregation and aggregation patterns. Indeed, the energy contributions , , and correspond, respectively, to the interface energies among the two populations + and −, +, and 0 state, and − and the vacancies. In complete analogy with the physics of surface tension [26], the coefficients , , and would be a kind of “social tension” among the interfaces. From a physical point of view, an interface is nonmiscible if the respective social tension, namely, , , or , is positive. Hence, because the system minimizes the energy, it minimizes the corresponding perimeter; otherwise, the interface would be miscible. Thereby the interface increases its length mixing the respective states.

The previous arguments allow us to give an explanation about why segregation is more efficient in the case of nonsymmetric interactions; for example, the segregation is stronger in Figure 1(c)(II), than in the symmetric case Figure 1(d)(II). In the nonsymmetric interaction case, the energy is not mandated strictly to decrease (see Figure 2(a)). Therefore, the system may evolve forever, searching in average the lowest energy as possible [23]. On the other hand, in the symmetric interaction case, the system reaches a frustrated state with a larger energy.

As a consequence of this, the whole patterns of Figure 1 may be understood according to (6). In Figure 1(a), that is, and , the perimeter energy reads . In this case, the energy decreases if the interface perimeter between and populations increases (see blue curve in Figure 2(b)), thus creating a labyrinth pattern in analogy with patterns of antiferromagnetic materials.

In the case of Figure 1(b), that is, and , the energy term is simplified to , being this a minima if, simultaneously, is maximized (see the upper red curve in Figure 2(b)) and is minimized (see the red curve in Figure 2(b)). In the former case, individuals arrange in stripe patterns maximizing the perimeter with the white vacancies, while in the latter case, individuals make clusters surrounded by a white layer of . Cases of Figure 1(c), , that is, , and Figure 1(d), (), produce segregation patterns. For these cases, the perimeter energy terms are similar: and for Figures 1(c) and 1(d), respectively. Both cases tend to minimize all perimeters among interfaces (see the green and orange curves in Figure 2(b)), but a notorious difference appears: in Figure 1(d) there is no / interface. As it can be seen in Figure 2 (curve (d)), the interface perimeter, , vanishes as the time goes by. This happens because forming this interface “costs” twice compared to forming the same interface of Figure 1(d); this produces that and regions are surrounded by a white layer which has no cost.

4. The Basic Social Structures

Under the same conditions of Figure 1, we have run numerical simulations for all 45 possible -matrices. The social structures may be classified in terms of the values of the -matrix based on three different aspects: () the particular social structure, that is, the pattern, and () the pattern power spectrum which provides information on the spatial correlations (see Section 5.3). The best classification scheme appears to be in terms of the parameters and . The first parameter gives the information of the attitudes towards their own group, while the second provides information about the attitudes of individuals towards the opposite group. Considering and , it makes classifying 25 different cases possible as shown in Figure 3 and in Table 1.

Table 1: Summary information on the patterns of Figure 3. The table is indexed by and . Each box contains the precise S-matrix and the exponent obtained from the corresponding Fourier spectrum (see Section 5.3). The cases labeled “ = none” mean that the slope may not be defined.
Figure 3: (Color online) classification of the attractors for different -matrices according to the and plane. The precise -matrices are given in Table 1. The 25 numerical simulations were done in a lattice with random initial configurations and . As expected from the energy expansion (6), the diagonal splits roughly the segregation-like and the aggregation-like patterns: / segregation behavior is observed below this line, while a heterogeneous attraction behavior is mostly observed above this line. We highlight the following social structures: segregation (S1), crossroads (S2), mutual suspicion (S3), social workers (S5), husband and wife (S6), boys and girls (S7), inclusion (S8), individualism (S9), and exclusion (S10), which are representative of the collective behavior.

The phase diagram of Figure 3 allows us to classify the basic social behaviors characterized by well-defined social structures. In his original work, Sakoda [11] summarized this rich dynamics by eight relevant behaviors with a particular interpretation of the real life situation (we keep original Sakoda’s terminology): crossroads, mutual suspicion, segregation, social workers, boys and girls, couples, and husband and wife. Based on the symmetries of the diagram, we suggest in the present work three new social behaviors: Inclusion, Repulsion, and Exclusion, which complete the whole picture. Naturally, the Sakoda model behavior is only simplified view of real social behaviors that may be absent in the current model.

We underline that all patterns as well as its evolution may be understood quantitatively in terms of the energy minimization principle introduced in the present paper. In particular, as expected, the parameter splits roughly the segregation-like and the aggregation-like pattern.

The description of the social behaviors presented in Figure 3 and its respective social structures are as follows:(S1) Segregation typifies attraction among individuals of the same group and repulsion among members of the other group (). In this case, the system minimizes the interface perimeter inducing segregation (see Figure 1(d)).(S2) Crossroads manifests also segregation; however in this case, the individuals of other groups do not interact. The -matrix reads (S3) Mutual suspicion is a characteristic typified by neutrality among members of the same group and repulsion towards the members of the other group. In this case the -matrix reads .(S4) Social climbers represents the interaction between a lower class and an upper class this is characterized by a -matrix .(S5) Social workers is represented in Figure 1(b). In this case, one group feels attraction by everybody while the other group feels repulsion for everybody. Its -matrix reads .(S6) Husband and wife is the case where the members of one group feel attraction by everybody, while individuals of the opposite group feel attraction for the former group and repulsion for individuals of the same group. In this case .(S7) Boys and girls is exactly the opposite of segregation, . This is an example of antiferromagnetic physics. In this case (), the perimeter is maximized [22].(S8) Inclusion manifests the social situations where all individuals want to be always together. In this case .(S9) Individualism is a case where the individuals want to be as far as possible from others. Here .

The two previous cases correspond to , which characterizes a “mixed phase”; indeed, though the former, Inclusion, maximizes the interface perimeter , the later, Individualism, minimizes it.(S10) Exclusion is a case where the individuals of one group want to be accepted by the opposite group, but the later do not accept the former group. Its -matrix is .

5. Quantitative Description of the Social Structures

5.1. Phase Diagram in the - Plane

All patterns may be understood quantitatively in terms of the energy minimization principle presented throughout this paper. Therefore, the energy seems to be a pertinent quantity for a characterization of the resulting social structures, at least concerning the social segregation [28]. For all cases, the energy reaches a steady state allowing us to characterize the social structures as statistical attractors.

In the special case , the energy expression (6) reads () From the numerical simulations, we conclude that the behavior of the corners and the interface length scale proportional also to , that is, the total energy scales, as , where , is a function of . The relationships among perimeters, corners, and the energy are more exhaustively developed in Appendix E.

Figure 4 shows a set of 40 of the 45 attractors (the missing 5 social structures are overlapped) classified in the - plane for . This figure shows that the energies are below an upper envelope, labeled by in Figure 4. It may be noticed the existence is an absolute energy maximum for all configurations for that correspond to and and an absolute minimum (, ). The social structures are spanned all over the - plane below this envelope and above the line for and the horizontal line for . The existence of two different behaviors for (aggregation) and (segregation) is a clear consequence of the existence of these two phases.

Figure 4: (Color online) phase diagram of 40 attractors in plane (normalized energy per occupied site, ) for (notice that the attractors , , , , and are overlapped by other attractors in the figure). These simulations correspond to a two-dimensional lattice running up to simulation steps in the case of short influence/long range movement evolution rule. Attractors with social interpretations are labeled (from (S1) to (S11)). Notice the pattern social climber (S4) missing in Figure 3. The missing five attractors are overlapped by others because of an approximate degeneracy of the energy for the same parameter . These are , , , , and . The curves and are sketches of the upper and lower bounds, respectively (see text).

Figure 4 presents a more complete landscape of the social structures than Figure 3, showing that energy appears to be a good characterization of social structures. However, the parametrization includes the energy that itself is a function on the interaction parameters. The energy spreads successfully the large number of degeneracies of the restricted classification in planes and of Figure 3. However, the energy itself depends on the parameters of the problem (see (A.3)); therefore, the characterization is visually attractive, but the characterization reveals out the real symmetries of the interactions and its consequences.

Lastly, we notice the existence of the 11th social structure, labeled by S11 and named as Neutrals. The social structures belong to Neutrals class where the individuals of different groups do not express neither rejection nor attraction among others. This behavior is represented by the interaction . Similarly, we identify the 9th group of social behavior (G9) named by social neutrality.

5.2. The Role of the Vacancy Fraction on the Social Structures

As it can be noticed in Figure 1, the final patterns depend on the empty spaces in the lattice, . In particular, we observe that the number of representative social structures diminishes as both and also in the limit . In the case of large number of vacancies, the individuals move almost freely through the space; as the vacancies diminish, finite space effects play a role. This is explained because of the smaller number of options that individuals have to move in the lattice. This effect modifies the basic interaction dictated, a priori, by (1). For instance, an individual that feel repulsion in free space would not necessarily follow a repulsive movement in a restricted environment (see Figures 1(c)(I) and 1(d)(I)). We investigated systematically the effect of on the final patterns.

Figure 5 summarizes the phase diagram of some attractors in the space. The reader may note that different vacancy concentrations, , do affect dramatically the final attractors. Nevertheless, one notices that the lines and characterize two critical lines which separate clearly three domains. For , we observe a heterogeneous phase, without long range order. For , a neutral phase appears with coexistence of a disordered phase. Lastly, for , a segregation phase with quasi-long range order appears. In particular, we notice that, in the special case of Schelling’s model, that is, , one observes that if we increase the number of total vacancies we recover the phenomenology already found in [20].

Figure 5: (Color online) examples of different attractors classified by parameter for different values of as indicated in the plot. Numerical simulations run in a two-dimensional periodic lattice of dimensions for simulation steps and considering short influence/long range movement rule.
5.3. Existence of Long Range Order on Social Structures

Generally speaking, it is observed as a common rule that all fractions present phases of segregation (ferromagnetic phase), heterogeneous attraction (antiferromagnetic phase), and a neutral or disordered phase (paramagnetic phase). More deeply, in a close analogy with phase transitions, we observe that the segregation-like patterns do manifest long range order (more precisely quasi-long range order), but the cases of heterogeneous attraction and disordered phases do not.

The quasi-long range order is readily quantified by the spatial correlation function which reads . It may be computed directly in terms of the Fourier transform of the spatial pattern: Here is the spatial pattern, is the linear size of a square lattice, and is an integer ranging within . A straight-forward calculation relates the correlation function with the spectrum. In the limit , we have thus yielding In practice, we compute directly the fast Fourier transform (FFT) of the spatial patterns and after performing an angular average, one gets the isotropic spectrum (in ): .

The spectrum allows us to distinguish quantitatively the existence of different social structures; moreover, a sudden change of these correlations may denote a phase transition from one regime to another. For instance, in the segregation regime that corresponds to a ferromagnetic case with well-defined domains, one expects that , where is a critical exponent. This power law behavior is a characteristic of quasi-long range order spatial correlation. Figure 6 presents the spectra of the four patterns of Figure 1 for . One notices that in the cases of (b), (c), and (d) one gets quasi-long range order with the following critical exponents: (b) ; (c) ; and (d) . Notice that (c) possesses an exponent close to which corresponds to Porod’s law, indicating that the observed phenomena is related to phase separation [29]. Besides the slight differences among the spectrum exponents, cases (b), (c), and (d) are systems that possess quasi-long range order, implying that faraway individuals are correlated among them. On the other hand, in the case of Figure 1(a), the spectrum is almost flat in the long length scales (short wave number ) but it shows a peak at the scale of the discrete lattice, revealing a short wavelength modulation (the stripe pattern). Table 1 summarizes the exponents for all cases corresponding to Figure 3.

Figure 6: (Color online) spectra of the social structure for the same patterns of Figure 1 row (II) (). In (a) and (b) the existence of a short wavelength structure is observed (in (b) this peak is labeled by an arrow).

6. Discussion

We have reported the dynamics and the social structures of the, almost forgotten, social interaction model proposed by Sakoda in the 1940s. The model is based on two different groups of individuals moving through a background of empty sites. It includes different attitudes that individuals assume (attractive, neutral, or repulsive) towards their own group and towards the other group, resulting in 45 fundamental interactions. All of them were studied in the present paper.

We observe that different attractive/repulsive/neutral social behavior was mainly characterized depending on the parameters and and . These parameters allow us to classify the pattern by its key features. The entire spectra of possible attractors may be studied in the frame of well-known physical analogies, like ferromagnetic (segregation) and antiferromagnetic (heterogeneous attraction) patterns, kinetics ordering, phase transitions, existence of quasi-long range order, and surface tension phenomena. In this sense, we propose a Potts-like energy that turns out to be a Lyapunov function in the case of symmetric interactions. This energy characterizes different attractors in terms of the parameters of the problem, namely, the initial fraction of vacancy places, the range of interaction, and the range of motion.

Surprisingly, the Sakoda model amalgamates its underlying behavior under an energy principle which, under the condition of symmetric interactions, mandates the system to equilibrium. This energy appears to be a common principle for all possible interactions because it decreases statistically in time, beyond the range of validity of symmetric interactions.

In conclusion, typical social structures and behaviors arise as a consequence of the interactions among individuals and the migration to a more attractive place. Depending on the intrinsic interactions we observe a variety of patterns (see Figures 3 and 4); some of them may be a manifestation of our societies. Despite the complexity of the interactions, expectative, and motions, an energy minimization principle surges as common law that rules the evolution of social structures. We suggest that the pertinence of the energy principle may deserve more attention from the social science community. We underline that the present model is only an approximative view of social behaviors. As an example, all individuals are considered of the same type, they are not distinguishable, and all of them are simply characterized by a single -matrix. Real societies are certainly more complex and harder to understand. The aim of this paper is to stimulate further quantitative description based on Sakoda’s approach.

We end this paper commenting on possible generalization of the Sakoda model studied here:(1)There are many variations of this model that may be considered. First, our results may be generalized to any complex network of interactions, more general graphs, and other space dimensions. A second possibility is to attribute a nonzero entry in the -matrix for the interactions among states and vacancies. In this case, the -matrix becomes a matrix, and the total number of possible interactions would be quite large: , but only 3903 of them are independent. Similarly, in the case of -states of the variable , the -matrix becomes interaction matrix and the total possible number of different cases to be considered would become huge: . In practice, in these cases, only a small sample of social structures could be studied. Therefore, an energy unified picture, as the one presented here, would be helpful.(2)Though we focused our work on the local range of influence of individuals and the long range of movement, the dynamical behavior for different ranges of influence and different movement of the individuals should be considered. Both the interaction range and the range of movements influence macroscopically the formation of the social structure. As a general rule, the typical length scale of the pattern depends on the range of interaction and/or the scale of the movement. The largest pattern arises for the case of long range movement/long range interaction. In this sense, agglomeration sizes are produced by the available information that individuals have access: individuals may explore more empty sizes if they consider a long range movement; on the contrary, short movements allow individuals to explore just few options of empty sizes; hence small clusters are expected. All possible cases of long/short interaction and long/short movement will be discussed in a separate publication.

Figure 7 summarizes the attractors for the same -matrices of Figure 1 but for the cases of long movement/long range interactions and short range movement/short range interactions and short range movement/long range interactions. As a general rule, the typical length scale of the pattern depends on the range of interaction and/or the scale of the movement. In those attractors that exhibit patterns of aggregation or segregation, long movements produce larger clusters of individuals than short movements. The range of interaction also plays a role on the scale of the patterns in the social structures. The largest pattern arises for the case of long range movement long range interaction (Figure 7(III)); the next intermediate size would be for the case of short range interaction long range movement (Figure 1(II)). For the cases of short range movement, the pattern’s scale is usually smaller. Figure 7(I) shows that the case of short range interaction-short range movement displays the smallest pattern’s size. Finally, Figure 7(II) corresponds to the case of long range interaction-short range movement with a slightly larger scale pattern. In this sense, agglomeration sizes are produced by the available information that individuals have access: individuals may explore more empty sizes if they consider a long range movement; on the contrary, short movements allow individuals to explore just few options of empty sizes; hence small clusters are expected.

Figure 7: Snapshots for different interaction and movement (rest) rules for the same cases as in Figure 1. (a) , (b) , (c) , and (d) . The first line (I) corresponds to the case of short range interaction and short range movement. The second row (II) corresponds to long range interaction/short range movement, and, finally the third line (III) corresponds to the case of long range interaction and long range movement. Simulations were run in a lattice with and for steps for (I) and (II) and steps for (III).


A. The Social Interaction Symmetries

Because may take three different values (), the total amount of different –matrices is 81. Therefore the set of all -matrices may be arranged in a matrix, namely, :The indexes of the matrix will be labeled by , and . Nevertheless, Sakoda’s interaction possesses symmetry under the transformation together with the exchange . This symmetry transformation reduces the 81 possible cases to 45. It is equivalent to take the lower or the upper triangular part of the matrix including the diagonal.

The Sakoda matrices can be classified by the parameters:The reader may note that any kind of linear combination among these allows classifying the 45 basic elements. However, our choice of the above parameters considers special invariances under the -transformation. Indeed, and do not change under , while changes its sign under this transformation.

B. Proof of (4)

Consider the exchange of the sites (occupied) and (empty); that is, and . Therefore, the energy for the initial configuration, , reads

We have split the terms with an index and/or . Because the terms and vanish and imposing the symmetry conditions, the initial energy is simplified to

Similarly the energy of the updated configuration readsAs before, all quantities involving vanish. The first term is identical in both expressions (B.2) and (B.3), because those terms do not involve the sites and . The resulting energy difference becomes which is exactly (4) of the paper.

C. Remark on the Potts Energy (3)

The -function equation (2) may be written as (we acknowledge J. P. Nadal for this remark):

Notice that in the second equality the first three terms are symmetric under the exchange , but the last one is not. After (C.1), the energy equation (3) takes the form of a spin-1 Ising model [20]:Because this energy is a symmetric sum on indexes and , the antisymmetric terms of (C.1) cancel out.

D. Derivation of the Energy-Perimeter Formula (6)

The terms in expression (C.2) contribute if and only if and are both different from zero. Hence, the sums run over producing, approximately, an effective sum only over sites. In the case of Moore neighborhood, the edges and corners formula stated in (6) may be derived after the following considerations: let us define regions , , and (which are not necessarily simply connected); the regions composed by the sates are , , and , respectively. These domains possess , , and individuals. Further, we define , , and as the interfaces among /, /, and / ; these interfaces consider the presence of different states in the same neighbors. The interfaces are geometrically expressed in terms of edges and corners of the sates and . (See Figure 8 as a sketch of the interfaces.)

Figure 8: (Color online) scheme for different interaction energy contributions in the case of the Moore neighborhood in a two-dimensional regular lattice. bond is represented in green, −− bond is drawn in orange, is drawn in purple, is drawn in sky-blue, is drawn in pink, and finally 00 bond is not drawn. (a) It shows a vertical interface among two populations and . The interface energy appears as a consequence of the asymmetry among the individual in the bulk, 8 individuals of the same type, and the individuals at the interface, 5 individuals of one type and 3 from the other. (b) It shows a general interface with both edges and corners. The reader may check the validity of (D.4) and (4).

The energy is nonzero only if and are in or . For the individuals in the bulk of (the individuals and whose neighbors are all ) the energy contribution becomes where differs from because of the boundary or interface terms. Rewrite this expression in terms of ; that is, Similarly, for the individuals and whose neighbors are all , the energy becomes

At the interface, the contribution on a site, say , to the energy comes from the eventual , , and interaction. Assume that, in its neighborhood, has individuals +, individuals of sign −, and individuals at 0 (naturally ), Then, after (C.2), the energy at an interface becomes In the case of a linear interface , the number of neighbors is , , and . To summarize, in the case of a linear interfaces one gets

Therefore the total energy becomes

Finally, one notices that as an effect of the discrete nature of the lattice, an interface may be composed of vertical and horizontal edges (joined by a corner). Hence, the total perimeter is formed by the contribution of edges and corners . A careful examination provides the relation for the perimeter in terms of edges and corners. One concludes that ; hence, we readily get (6).

E. The Edge-Corners Relationships

The numerical simulations provide evidence of correlations among the number of edges and corners of an interface for all concentration vacancies together. Figure 9(a) plots the edges versus the respective corners for the 45 Sakoda cases. This linear relation among the perimeter and the corners has only a statistical sense which follows from the obvious equivalence among the number of edges and the number of corners in a diagonal interface. Naturally there are patterns with a large/small perimeter but with a small/large number of corners which are out of this description; for example, the heterogeneous attraction, , and the repulsion are far away from the trend line.

Figure 9: (a) Plots of the rescaled interface edges versus the rescaled number of corners for all the types: are symbolized by and are symbolized by ■ and for case. The colors represent the five different vacancy concentrations: the cyan markers stand for , the red ones stand for , the green ones are for , blue ones stand for , and orange ones stand for . One notices that the correlation is linear: . The numerical simulation considers periodic lattice and includes the 45 different cases and 5 different initial fractions of vacancies. (b) plots the total rescaled number of edges versus ; (c) the total rescaled number of edges versus ; and finally (d) plots the total rescaled energy versus . In plots (b), (c), and (d) the surfaces were obtained after an interpolation with a quadratic fit of the form , where a, b, c, d, e, and f are fitting constants. The simulations were developed in a two-dimensional periodic lattice of size , considering short influence/long movement range evolution rule and with a running time of steps.

Figures 9(b) and 9(c) show the rescaled perimeter interface, , and number of corners, , as a function of and , respectively.

Figure 9(d) also presents the relation of the attractor energies as a function of the parameters and , for the 45 -matrices. A dependence of the rescaled energy as a function of for various values of is observed. From the data of Figure 9(d), it is possible to conclude that the behavior of the corners and the interface length scale are proportional to ; that is, , where is a function of . Similar relations hold for the number of edges and corners, hence the interface length . The scaling for the energy follows directly after these relations: because the “bulk” energy becomes in the special case of , the perimeters scale as , hence yielding the total energy scales as .

F. The Long Range Interaction Case

In this case, let us consider the interaction exchange coefficient a function that depends explicitly on the spatial distance among sites and , being , the Euclidean distance. Therefore, (2) of the main text may be written as

To avoid any problem with a divergent self-interaction, more precisely, case in sum (F.1), assume that in our model. Though could be any kind of function, as a general rule, consider that is a decreasing function in its argument; that is, the shorter the distance is, the stronger the influence is.

In the present chapter, for the long range interaction, the potential function is proposed as follows:

Originally, Sakoda reported his results considering , because he argued that the weight contribution of faraway individuals should be quite strong; that is, the social interactions are of “extremely” long range. In general, because strong long range interactions do not evolve to an equilibrium (at least in short time) these cases are avoided and limited to “weak” long range interactions, say in two-dimensional space. In this case it is expected that the robust behavior does not depend explicitly on the detail of the individual interaction, something already noticed by Schelling [1].

Figure 10 presents a comparison between the case on long range interaction/long range mobility in the particular case of and the Sakoda strongly long range interaction . The simulations show that the case does not depend strongly on the size of the system that is confirmed comparing Figure 7 (long influence/short movement case) with Figure 10.

Figure 10: Snapshots for the same cases as in Figure 1 of the paper considering the extremely long range interaction : (a) , (b) , (c) , and (d) . These snapshots correspond to a lattice size and . These simulations run up simulation steps.

Competing Interests

The authors declare no competing interests or any conflict of interests regarding the publication of this manuscript.


The authors thank Cyrille Piatecki for bringing Sakoda’s original reference to their attention and A. Valdivia who kindly offers them computational power at the Universidad de Chile. This work is supported by Núcleo Milenio Modelos de Crisis NS130017 CONICYT (Chile) and Eric Goles and Sergio Rica also acknowledge the CONICYT Grant BASAL-CMM and the FONDECYT Grants 1140090 and 1130709, respectively. Pablo Medina acknowledges the support from Colciencias (Colombia) through their doctoral funding program.


  1. T. C. Schelling, Micromotives and Macrobehavior, Norton, New York, NY, USA, 2006.
  2. E. Wilson, Sociobiology: The New Synthesis, Harvard University Press, Cambridge, Mass, USA, 1975.
  3. E. Wilson, The Insect Societies, Harvard University Press, Cambridge, Mass, USA, 1971.
  4. L. Quevillon, E. Hanks, S. Bansal, and D. Hughes, “Social, spatial, and temporal organization in a complex insect society,” Scientific Reports, vol. 5, article 13393, 2015. View at Publisher · View at Google Scholar
  5. J. K. Hamlin, K. Wynn, and P. Bloom, “Social evaluation by preverbal infants,” Nature, vol. 450, no. 7169, pp. 557–559, 2007. View at Publisher · View at Google Scholar · View at Scopus
  6. R. Hegselmann, “Modeling social dynamics by cellular automata,” in Computer Modeling of Social Processes, W. B. G. Liebrand, A. Nowak, and R. Hegselmann, Eds., p. 37, Academic Press, London, UK, 1998. View at Google Scholar
  7. C. Castellano, S. Fortunato, and V. Loreto, “Statistical physics of social dynamics,” Reviews of Modern Physics, vol. 81, no. 2, article 591, 2009. View at Publisher · View at Google Scholar
  8. L. Guo, Y. Cheng, and Z. Luo, “Opinion dynamics with the contrarian deterministic effect and human mobility on lattice,” Complexity, vol. 20, no. 5, pp. 43–49, 2015. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  9. L. Gauvin, J. Vannimenus, and J.-P. Nadal, “Phase diagram of a Schelling segregation model,” European Physical Journal B, vol. 70, no. 2, pp. 293–304, 2009. View at Publisher · View at Google Scholar · View at Scopus
  10. J. Sakoda, Minidoka: an analysis of changing patterns of social interactions [Ph.D. thesis], University of California, Berkeley, Calif, USA, 1949.
  11. J. M. Sakoda, “The checkerboard model of social interaction,” The Journal of Mathematical Sociology, vol. 1, no. 1, pp. 119–132, 1971. View at Publisher · View at Google Scholar
  12. T. C. Schelling, “Models of segregation,” American Economic Review, vol. 59, no. 2, pp. 488–493, 1969. View at Google Scholar
  13. T. C. Schelling, “Dynamic models of segregation,” Journal of the Mathematical Society, vol. 1, no. 2, pp. 143–186, 1971. View at Publisher · View at Google Scholar
  14. M. Pollicott and H. Weiss, “The dynamics of Schelling-type segregation models and a nonlinear graph Laplacian variational problem,” Advances in Applied Mathematics, vol. 27, no. 1, pp. 17–40, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  15. D. Stauffer and S. Solomon, “Ising, Schelling and self-organising segregation,” The European Physical Journal B, vol. 57, no. 4, pp. 473–479, 2007. View at Publisher · View at Google Scholar
  16. D. Vinković and A. Kirman, “A physical analogue of the Schelling model,” Proceedings of the National Academy of Sciences of the United States of America, vol. 103, no. 51, pp. 19261–19265, 2007. View at Publisher · View at Google Scholar
  17. R. Pancs and N. Vriend, “Schelling's spatial proximity model of segregation revisited,” Journal of Public Economics, vol. 91, pp. 1–24, 2007. View at Google Scholar
  18. A. Singh, D. Vainchtein, and H. Weiss, “Schelling's segregation model: parameters, scaling, and aggregation,” Demographic Research, vol. 21, pp. 341–366, 2009. View at Publisher · View at Google Scholar · View at Scopus
  19. S. Grauwin, E. Bertin, R. Lemoy, and P. Jensen, “Competition between collective and individual dynamics,” Proceedings of the National Academy of Sciences of the United States of America, vol. 106, no. 49, pp. 20622–20626, 2009. View at Publisher · View at Google Scholar · View at Scopus
  20. L. Gauvin, J. Nadal, and J. Vannimenus, “Schelling segregation in an open city: a kinetically constrained Blume-Emery-Griffiths spin-1 system,” Physical Review E, vol. 81, no. 6, Article ID 066120, 2010. View at Publisher · View at Google Scholar
  21. N. Goles Domic, E. Goles, and S. Rica, “Dynamics and complexity of the Schelling segregation model,” Physical Review E, vol. 83, no. 5, Article ID 056111, 13 pages, 2011. View at Publisher · View at Google Scholar
  22. P. Collard, “Beyond the schelling's segregation: is it equivalent to be repulsed by dissimilar rather to be attracted by similar?” in Advances in Artificial Life, ECAL 12, pp. 480–487, 2013. View at Google Scholar
  23. V. Cortez, P. Medina, E. Goles, R. Zarama, and S. Rica, “Attractors, statistics and fluctuations of the dynamics of the Schelling's model for social segregation,” The European Physical Journal B, vol. 88, no. 1, article 25, 2015. View at Publisher · View at Google Scholar · View at Scopus
  24. R. B. Potts, “Some generalized order-disorder transformations,” Mathematical Proceedings of the Cambridge Philosophical Society, vol. 48, no. 1, pp. 106–109, 1952. View at Google Scholar · View at MathSciNet
  25. F. Y. Wu, “The Potts model,” Reviews of Modern Physics, vol. 54, no. 1, article 235, 1982. View at Publisher · View at Google Scholar · View at MathSciNet
  26. P. G. De Gennes, “Wetting: statics and dynamics,” Reviews of Modern Physics, vol. 57, no. 3, pp. 827–863, 1985. View at Publisher · View at Google Scholar · View at Scopus
  27. E. Goles and S. Martínez, Neural and Automata Networks: Dynamical Behavior and Applications, vol. 58 of Mathematics and Its Applications, Kluwer Academic, Dordrecht, 1990. View at Publisher · View at Google Scholar · View at MathSciNet
  28. V. Cortez and S. Rica, “Dynamics of the Schelling social segregation model in networks,” Procedia Computer Science, vol. 61, pp. 60–65, 2015. View at Publisher · View at Google Scholar
  29. A. Bray, “Theory of phase-ordering kinetics,” Advances in Physics, vol. 43, no. 3, pp. 357–459, 1994. View at Publisher · View at Google Scholar