Undergraduate Tutorial for Simulating Flocking with the Vicsek Model

The biological world is full of fascinating collective behaviors across an enormous range of length scales. At the nanometer-micrometer length scale, cytoskeletal proteins organize and exhibit long-range order while coordinating motion leading to cell division (1). In vitro mixtures of cytoskeletal proteins also develop fascinating long-range patterns (2–4).

Ordering and collective motion also occur within expanding bacterial colonies on agar plates (5, 6). Group motion of bacteria within bacterial suspensions is capable of affecting mechanical properties, such as the viscosity of a fluid (7). Furthermore, bacterial motility is coupled to biofilm growth (8, 9).

Eukaryotic cells coordinate with their neighbors in both their motion and the forces they generate when closing wounds (10–12). Within epithelial sheets, cells exhibit interesting and complex patterns of motion while maintaining a mechanically robust cell layer (13–15). Although those coordinate behaviors contribute positively to human health, coordinated behavior in the form of cancer metastasis does not (16, 17).

Collective behaviors reach significantly larger length scales as well. Most notably, collective behaviors are observed in the flocking of birds and the schooling of fish. Order and complex dynamics of humans within dense crowds (18) and animals within migrating mammal herds (19) are also observed.

The emergence of group motion from individual motion, which is present from micrometer length scales to hundreds of meters, is common to the biological world. Although the specific details associated with collective motion likely differ at different length scales, a comprehensive understanding of these behaviors would be helpful in defining the effective interactions between individuals that lead to these interesting behaviors.

The Vicsek model (20), proposed in 1995, is a simplified computational model that captures the emergence of collective behavior. In this model, objects move through space and tend to orient themselves with their neighbors. Thus, the details of the interaction between individuals change as an object accumulates more neighbors—a hallmark of an emergent system. The Vicsek model is a valuable tool to study collective behaviors at any length scale because it makes no assumptions about the origin of the interaction between individuals.

The utility of this model is also to help identify new characteristics of collective behavior. For example, short-range interactions between individuals that lead to network-level coordination must transfer information. As such, the Vicsek model is a tangible method to study information flow (21, 22). Additionally, comparisons between group-level properties in the model and experimental systems can elucidate the true interactions that likely deviate from the simplified Vicsek interaction (23, 24).

The Vicsek model also predicts a phase transition from a disordered state (individuals acting independently) to an ordered state (collective behavior). Most phase transitions that students learn about in a typical undergraduate curriculum are thermodynamic in origin, such as crossing of solid/liquid or liquid/gas phase boundaries by changing the temperature or pressure of water. However, this example is a different type of transition—a kinetic phase transition.

In this manuscript, we present a worked-through tutorial for building the Vicsek model from scratch. This example is the product of a research project in the lab of A. Pasha Tabatabai (APT), where APT led undergraduate physics majors, coauthors MacQuarrie Thomson (MT) and Reece Keller (RK), to build their own active particle simulations, which included the Vicsek model. A unique feature of this tutorial is that MT and RK played an integral role in identifying particular skills that they needed to learn to build their models independently (i.e., skills that may not be reinforced in a standard undergraduate curriculum). In response, we highlight and provide support in these areas by generating specific problems within the tutorial that address these skills. It is our hope that the discrete nature of the questions within this tutorial with clear actionable items makes it accessible to all students regardless of their background in computational physics or collective motion. Additionally, we provide detailed discussion in other sections for the topics of periodic boundary conditions (IV.B.3), troubleshooting as a skill (IV.B.4), and defining steady state (IV.B.5). A conservative estimate is that this tutorial can be completed within an academic quarter or semester.

We conclude the manuscript with a discussion of some current work with flocking models and how an undergraduate reader would implement these alterations to perform novel research in emergent behaviors.

The Vicsek model is a valuable pedagogical tool for many reasons.

• (a) An interdisciplinary approach to science encourages students to build a well-rounded skill set making them better prepared for scientific careers. Although the basic Vicsek model obscures biologically relevant details in the interactions between objects, it is possible to adapt this model to better reflect true biological systems.

• (b) This model provides an approachable strategy to connect an extremely complicated observation in the real world (i.e., flocking) with basic mathematical and physical rules. This model is computationally simpler and easier to integrate than other active simulations where a true equation-of-motion with interparticle interactions is required (25, 26).

• (c) This model provides an example of a phase transition. It would be beneficial for students to see this model after taking a chemistry or thermodynamics course. In this way, comparisons can be made to a thermodynamic phase diagram. Within a physics course, this could be compared with the phase transition in the ferromagnetic two-dimensional (2D) Ising model. Introducing the Vicsek model (a non-thermodynamic phase transition) provides an opportunity to explore the transition through simulation. Implementing this within a noncomputational course could be accomplished by providing code to students and encouraging them to change variables.

• (d) Computational skills are increasingly important for students interested in a career in science. Little computational background is needed to build this model, particularly with the assistance of the tutorial. As such, exploring the Vicsek model can be implemented as part of a thermodynamics course, statistical physics course, or a biological physics course to learn new physics. Alternatively, this example could be implemented as a module within a computational physics course to emphasize computational skills in the context of a physics problem. Depending on the course goal, students can be provided with sections of prebuilt code.

• (e) The ability to troubleshoot issues in research is a skill that translates across all disciplines and careers. In this manuscript and associated tutorial, we discuss troubleshooting strategies and break down the problem into small components that build on each other. Providing pre-built code to students, as mentioned above, may limit the troubleshooting growth of students.

Often times students learning challenges are defined by the instructor, which may or may not be aligned with actual student challenges. As such, MT and RK identified skills that they needed to learn to implement a computational model. We hope that this style of tutorial is helpful for future undergraduate students to work independently through the manuscript.

Other pedagogical tools based on the Vicsek model exist. For example, the authors are aware of a video describing the original paper (20) in detail (27), as well as fully worked example code accompanied by blog-style commentary (28). Although resources like these are helpful, we believe that students will more deeply understand the Vicsek model by building their own code from scratch as opposed to referencing the prebuilt code of others.

Although it is possible to write these simulations in any computer language, we chose to write in Python 3 because it is free, open source, and taught in the introductory computer science course at our institution. We downloaded the Anaconda distribution and worked within a Jupyter notebook.

We provide a tutorial in the Supplemental Material that describes building the Vicsek model from scratch by using questions that develop and test the model step-by-step. We provide solutions to each question so that students can be sure that they are moving in the right direction at each step.

The angular noise Δθ varies for each object at each time step and is chosen from a uniformly distributed random angle within the interval [–η/2, η/2]. Therefore, η sets the possible values of the random noise Δθ. The effect of this angular noise is to add some uncertainty in the trajectory of a given particle, and its effects can be seen by looking at an isolated particle where becomes θ_i(t) (Fig 3).

In a system with many objects, the term in Equation 2 acts to align neighboring objects, whereas Δθ acts to randomize orientations. When η is small, the ordering effect dominates, and long-range order is present (Fig 4).

Thus we see that by incorporating these rules, we capture an essential result of the Vicsek model: that the noise η controls a dynamic first-order phase transition from ordered flocking (high v_a) to disorder (v_a → 0), even when the number of objects per unit area (i.e., the density) is kept constant (Fig 5).

When fully describing a phase transition, it is common to determine exactly how an order parameter changes at this transition point. In the case of the Vicsek model, a noise level η_c separates the regime of collective motion from the disordered regime. In fact, this noise value is also a function of system size η_c(L), as seen in Figure 5. In the original publication, Vicsek et al. (20) treated this transition as continuous (i.e., a second-order phase transition like the 2D Ising model). However, later works with larger systems argue that the canonical Vicsek model is a discontinuous (i.e., a first-order) phase transition (29) and the transition occurs at a noise level that is also a function of system density η_c(ρ) (30).

One reason that the Vicsek model is interesting to study is that the long-range order associated with large values of the order parameter v_a are not possible in 2D for equilibrium systems. The Mermin–Wagner theorem states that a system with short-range interactions at thermal equilibrium cannot lead to a long-range order (31, 32). In the Vicsek model, the interactions are short range (r₀ ≪ L), yet the steady state is capable of achieving long-range order (v_a ∼ 1) because self propulsion (v ≠ 0) is a nonequilibrium feature. An exercise in the tutorial demonstrates that measurements of v_a for systems with and without swimming are very different despite all other parameters being the same.

Computational abilities vary widely across an undergraduate cohort partially because computational courses are introduced at different points in the undergraduate curriculum, and different departments have different major requirements for computational courses. As such, many of the challenges in developing this model are rooted in the learning curve associated with Python.

The tutorial assumes students know the following Python skills before starting the tutorial:

• (a) How to download and install the Anaconda distribution (www.anaconda.com) of Python 3 and work within a Python environment. We recommend Spyder or Jupyter notebooks.

• (b) How to import libraries such as matplotlib and NumPy.

• (c) How to use PyPlot from matplotlib to display and save plots.

• (d) How to write a function.

Students who are unfamiliar with these items can find myriad online resources to help them prepare for this tutorial.

The tutorial is designed to approach particular computational problems step-by-step. The tutorial covers the following computational topics in the context of problem solving:

• (a) Importing and using built-in Python functions.

• (b) Random number generators.

• (c) Using loops and if statements.

• (d) Generating and displaying images and plots.

• (e) Calling multiple functions and using return statements within functions.

The use of periodic boundary conditions in computational simulations is common. When simulating something, you want to avoid the effects of a boundary. If a boundary exists, you have to consider new rules for interaction with that boundary. In fact, what you really want is an infinitely large simulation so that the effects of these boundaries are minimal, but that would require infinite computing resources (and take infinitely long to run).

We use periodic boundary conditions such that objects leaving the bounding area reappear on the opposite side of the simulation area. As such, objects simply pass through the boundary and do not interact. Implementing periodic boundary conditions is an explicit step in the tutorial.

A primary challenge in building these models, and research in general, is finding troubleshooting methods for suspected errors. Within the typical structure of an undergraduate physics course, students tend to determine whether their solutions to problems are correct on the basis of a comparison with a known solution (either given or from a peer). However, in research scenarios, a known solution may not exist. As such, students must find ways to convince themselves and their advisor that their results are correct, determine how to correct their results, or both.

This process deviates significantly from the learning process within a traditional physics course. In these courses, students tend to submit assignments and receive a grade—end of transaction. However research is iterative, and progress relies on reflecting on the work done and implementing improvements.

As a consequence of this difference between research and coursework, troubleshooting is not an explicit skill learned in courses. Undergraduate research is important for this reason. We emphasize here that troubleshooting is an active part of the research process. Typical methods used to confirm that the model is implemented correctly (e.g., visual checks, building toy models, and taking limits of results) are demonstrated in the tutorial.

In a typical simulation, particles are distributed in a region of space with random positions and orientations. The system is allowed to evolve, and in this evolution process the positions and orientations may become more correlated across objects depending on the experimental conditions. The goal is to determine the effect of parameters on the ultimate behavior of the system that is independent of the initial conditions and configuration.

As such, we separate the behavior of our system into two regimes: the “transient” and the “steady state.” In the steady state behavior, fluctuations in a measured quantity such as the order parameter occur, but the mean value is constant. In the preceding transient regime, a measured order parameter fluctuations and the mean are changing in a way that depended on the initial configuration. Here, we only make measurements in steady state. A strategy to determine whether the simulation has reached the steady state regime is to compare averages of the order parameter over different time windows. If the average is unchanged, then each time window is in steady state. Note that reaching steady state is not necessarily the same as reaching thermodynamic equilibrium.

Although the simplicity of the Vicsek model is sufficient to capture this order–disorder transition, the model is still being adapted in current research. We highlight a few recent works.

The primary assumption in the standard Vicsek model is that objects interact with all other objects within a given radius. However, the visual information collected by a bird, for example, is not isotropic (i.e., the bird cannot see behind itself). To this end, nonisotropic interactions with limited fields of view, as shown in Figure 6, have recently been integrated into the Vicsek model (33–35).

It is important to consider that true interactions between animals may not be governed by metric interactions (i.e., how far the animals are apart). For example, starlings in a flock have been shown to interact topologically (individuals interact with some number of nearest neighbors, regardless of distance) (23), and fish in a school use line-of-sight interactions (24).

Another principle assumption of the standard Vicsek model is that objects move with a constant speed at all times. This model has been adapted to better mimic the variability of fish speed, which has been determined to depend on the local alignment of neighbors (36).

Additionally, most real systems exist in three dimensions, whereas this manuscript and the tutorial are explicitly done for the 2D Vicsek model. This can be extended to three dimensions.

Furthermore, self-propelled particle models such as the Vicsek model are used as a reference to other active models like the vertex model when investigating the order–disorder transition (37).

We have motivated the reasons to study the Vicsek model in the undergraduate curriculum both from the perspective of a student researcher and as an educator. We have defined the details of the Vicsek model and implemented this model to show a canonical result—a kinetic phase transition between order and disorder that stems from local interactions between active objects. We provide discussion about concepts that are necessary to build this model that may be new intellectual territory for undergraduate students, and we accompany a tutorial to help students build the model step-by-step. Finally, we have discussed a few ways in which current researchers are adapting the Vicsek model for their own research.

It is our hope that undergraduate students feel comfortable approaching this model and that our provided resources, in addition to the original manuscript (20), are sufficient to assist in an undergraduate student building their own model from scratch.

A supplemental tutorial that builds the Vicsek model is available at https://doi.org/10.35459/tbp.2022.000227.s1.

APT designed and conceived the project and wrote the manuscript. MT, RK, and APT performed simulations and analysis. MT and RK identified and defined learned skills and concepts.

APT was supported by an M. J. Murdock Charitable Trust Award (FSU-201913717) and a National Science Foundation Award (2137509). The authors declare no conflicts of interests.