Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Theor Biol. Author manuscript; available in PMC 2013 September 21.
Published in final edited form as:
PMCID: PMC3409318

Modeling Local Interactions during the Motion of Cyanobacteria


Synechocystis sp., a common unicellular freshwater cyanobacterium, has been used as a model organism to study phototaxis, an ability to move in the direction of a light source. This microorganism displays a number of additional characteristics such as delayed motion, surface dependence, and a quasi-random motion, where cells move in a seemingly disordered fashion instead of in the direction of the light source, a global force on the system. These unexplained motions are thought to be modulated by local interactions between cells such as intercellular communication. In this paper, we consider only local interactions of these phototactic cells in order to mathematically model this quasi-random motion. We analyze an experimental data set to illustrate the presence of quasi-random motion and then derive a stochastic dynamic particle system modeling interacting phototactic cells. The simulations of our model are consistent with experimentally observed phototactic motion.

Keywords: Cell-motility, Agent-based models, Cyanobacteria, Group dynamics

1 Introduction

Unicellular microorganisms have evolved to live in variable and extreme environments. Some are capable of intercellular signaling and appear to utilize group dynamics to achieve desired actions, such as moving toward a food source [21] or towards light [4]. These group dynamics often result in emergent patterns which can be modeled and analyzed using mathematics. Synechocystis species PCC6803 (hereafter Synechocystis sp.) is a commonly studied unicellular freshwater Pasteur Culture Collection cyanobacterium that displays the ability to move toward light, a phenomenon referred to as phototaxis, forming finger-like projections in the direction of a light source. Synechocystis sp. is a model organism for studying phototaxis in a laboratory setting and extensive genetic and microscopic analyses have been carried out to characterize the molecular bases for motility [3]. It has been demonstrated that this surface dependent motility requires Type IV pili, photoreceptors [20] and a number of other proteins necessary for phototaxis [3]. When wild type Synechocystis sp. is exposed to light, cells begin to move, although not necessarily in the direction of light. They instead form small aggregations of cells and eventually, with a time delay, cells may move toward light. At the front of a spot of plated cells, cells align along the boundary of the spot before forming the characteristic finger-like projections or swarms of cells [8]. Yet the motion of individual cells is not as directed toward the light source as is the observed group behavior. Individual cells instead display a quasi-random motion, that is, they move in seemingly random directions. In this paper, we mathematically model the local interactions of Synechocystis sp. and observed quasi-random motion in order to address a series of questions:

  1. If motile cells are not moving exclusively toward the light, are they moving in random directions or do they move following other non-random rules of motion?
  2. Is there a characteristic distance, beyond which cells can no longer sense the presence and behavior of neighbors?
  3. If such a distance exists, how do the patterns of motion vary with respect to this distance?

An extensive amount of mathematical research has been conducted in related fields. Specifically, we would like to mention some works on animal flocking and on chemotaxis. The Couzin-Vicsek model of flocking (and its many extensions) allows animals or individual agents to be repelled by near neighbors, align with the average directional heading of not-so-near neighbors, and be attracted to far neighbors [9, 36]. This model has lent itself to many applications as well as thorough mathematical analysis, for example see [12]. Cucker and Smale offer a dynamical system which models the flocking of opinions in human networks [10, 11]; this model has also been subjected to significant analysis, for example see [15, 16]. Similar models of flocks and schools have been developed for a variety of self-propelling agents such as birds and fish, e.g. [2, 19, 27, 28, 32].

Chemotaxis, i.e., motion towards a chemical attractant, is also a field that has been extensively studied in recent decades, starting from the celebrated work of Patlak, Keller and Segel, [23, 33]. For completeness, we refer the interested reader to the following papers and to the references therein [1, 17, 18, 31, 34]. Before migrating toward the light source, Synechocystis sp. form small aggregations whose location is unrelated to a spatial concentration of a chemoattractant as is the case with chemotactic motion. Consequently, most of the mathematical modeling and analysis on the topic is irrelevant in the present context.

Furthermore, phototaxis has not been extensively subjected to mathematical modeling. Only a few models of phototaxis have been developed, for example see [7, 29]; however, these models do not consider the intercellular group dynamics. A recent agent-based model of phototaxis considers cell interactions by transmission of light by individual cells [13]. In a series of papers we have developed several families of mathematical models for describing phototaxis. [6, 8, 24, 25, 26]. In all these models, the primary focus was on the phase of the initiation of the motion towards a light source (including the associated time delay) and the resulting overall migration of the colony of cells towards the light source (including the modeling of the finger formation). The emphasis of these analyses was on the role of the group dynamics, as opposed to what can be associated with the behavior of the individual cell. Missing from these analyses was the description of what happens in regions of low to medium cell density. The purpose of this paper is to develop, present, and study mathematical models for such regions.

Our mathematical models follow the time-discrete dynamics of a finite set of particles that are interacting in a two-dimensional domain according to rules that involve certain random terms. The rules for the local interactions between particles, are based on our experimental observations given by the analysis of time-lapse movies of the bacteria under a plethora of controlled conditions.

The structure of the paper is as follows. In section 2 we present the biological background, describe the experimental setup, demonstrate some of the experimental results, and formulate a list of observations that lead us in developing the mathematical models. The mathematical models are then presented in section 3. Two mathematical models are described. Our first model is a discrete-time model of local interactions that allows cells to keep moving in their previous direction, stop moving, or move in the direction of one of their neighbors. The second model increases the randomness in the motion. This time, cells can elect to change their direction of motion, but only if there are other cells nearby. Simulations of these models are then presented in section 4. These simulations show that the models, in particular the local-interactions model, provide results that replicate the experimental observations. Concluding results are given in section 5.

2 Biological background

In our study, we consider a phototactic microorganism Synechocystis sp. This unicellular cyanobacterium has been studied extensively by time-lapse video microscopy (TLVM). The experimental observations that serve as the basis to the mathematical models that follow, were generated in the Bhaya lab. Time-lapse movies of moving bacteria were acquired under an optical microscope and then analyzed as described below.

In Figure 1, we show several examples of some of the experimentally observed characteristics of motion. Each image contains wildtype Synechocystis sp., exposed to a directional light source coming from the top left of the domain. In Fig. 1(a), the beginnings of finger formation are observed. In Fig. 1(c), the same part of the spot is observed 24 hours after (a) with significantly more prominent fingers. The cells observed in Fig. 1(b) and 1(d) are taken from the same experiment at different time points and locations on the plate and highlight the wide range of observable cell behavior. Fig. 1(b) illustrates the effect of the surface on movement – cells preferentially move on a pre-wetted surface which contains extracellular polysaccharides produced by the cells. In (d), cells form small aggregations and move toward the light-source. Such aggregations of cells can be seen in all four images. We also observe the tendency of cells to align along the boundary between the pre-wetted and agarose-only surface in Fig. 1(a), (c), and (d). In each figure, we note the tendency of some bacteria to move in pairs, although it is important to note that in some cases, these pairs of cells are cells which have divided but have not separated.

Figure 1
Experimental results illustrating (a) beginnings of finger-like projection formation, (b) surface dependence, and (c/d) aggregation and finger formation. The bacteria within a matrix of extracellular polysaccharides appear as white dots or circles. The ...

Not illustrated in Fig. 1, the wavelength of light and density of cells also impact finger formation and general cellular movement. Additionally, directional cell movement is typically delayed, i.e. after cells are plated, they do not appear to start moving in the direction of light until a few hours have passed. Furthermore, cells have been observed to sense the presence of and move toward other cells over a large distance, as is shown in [25]. This seems to imply the existence of at least one communication mechanism between cells.

To analyze the experimental data, we used ImageJ, a public domain image processing software produced by the National Institutes of Health, and a plug-in called Particle Tracker developed by Sbalzarini and Koumoutsakos [35]. Particle Tracker determines the position of the cells in each picture, as demonstrated in Fig. 2 and then links the cell positions together to determine trajectories for each cell. Similar algorithms to [35] are utilized by other groups (such as, e.g., the algorithm used to track fish in [22]). In Fig. 2, we demonstrate the cell-position-identifying capabilities of Particle Tracker. We observe that most of the cells in the image are correctly identified.

Figure 2
An example of using the Particle Tracker plugin for ImageJ. Cells that were identified in this frame have red circles around them. The light source imposes a force on the cells in the direction indicated on the picture. This same directional force is ...

After the cells have been identified in a series of images, we further use the plug-in to track cells as they move. In Fig. 3, the different subfigures are shown at consecutive times, taken 50 frames, or 100 seconds, apart. These cells, observed under an optical microscope, are exposed to a directional light source, the force pointing toward the upper left corner of each image, as illustrated in Fig. 2 and Fig. 3(g). Particle Tracker is used to identify trajectories over 453 images, taken at 5 second intervals. In Figure 3, we show the 44 trajectories which are at least 100 frames in length.

Figure 3
Example of trajectories produced using Particle Tracker. (a) Indicates cells which were stationary for the entire 2000 s. (c) Illustrates a pair of cells which move together for frames (a)–(e). (e) Highlights a cell which turns to move toward ...

In Figure 3, a few patterns of motion can be identified using the trajectory analysis. In Fig. 3a, a couple of cells are identified as stationary. The small dot-like trajectories of these cells illustrate that these cells do not move significantly from their original positions. Cells moving as a pair are marked in Fig. 3c. In Fig. 3e, we see a cell that abruptly changes its direction to move toward a neighboring cell. Additionally, in comparing the series of figures, it is possible to identify at least one cell moving in the direction of light, although this is clearly not the case for the majority of cells.

It is clear from Fig. 3 that the particles do not obey a simple rule of motion. They do not all move towards the light source. Yet, the motion is very different than a random walk. The particles can maintain their direction of motion for several minutes, after which they may change their direction of motion, remain stationary, or adhere to a neighboring particle. In any event, this motion does not correspond to any of the classical variations of a random walk.

Is there any bias in the direction of their motion? To study this question we use numerical data generated by Particle Tracker to compute a histogram illustrating the distribution of angular direction of cellular movement. The data set we considered this time is shown in Fig. 4(a). We cropped this picture to the small rectangle shown in the northeast corner of the picture and used Particle Tracker to detect and track particles. Using a small section of experimental video allows us to isolate the influence of light and focus on the local cell-cell interactions. The light source in this picture is to the northwest, or 34π radians (about 2.4), and directly above and to the left of the small selection, in the region with a visibly higher cellular density, is a finger moving toward the light. A scatter plot of the distances and directions travelled by the particles from the selected rectangle is shown in Fig. 4(b). Fig. 4(c) is a histogram for angular directional movement of cells between time frames. One might expect that the cells in the selected area would also be traveling in the direction of the finger; however, as can be seen in the histogram in Fig. 4(c) this is not the case. The cellular movement is almost uniformly distributed with no apparent bias toward light and possibly a slight bias toward moving to the right during this 10 minute interval.

Figure 4
An example of using Particle Tracker data to determine the direction of the cells. In (a), we see a large mass of cells with a smaller region selected. This figure shows one image taken from a sequence of 301 with 2 seconds between each frame and spanning ...

The focus of this work is on regions of low to medium density, where we observe a quasi-random motion. Cells do not necessarily move only towards light, but also in a different, seemingly random motion as demonstrated by the particle trajectories in Fig. 3. Even at the base of a directed finger as seen in Figure 4, the overall movement of cells is not directed toward the light source. This quasi-random motion has also been observed within fingers on short time scales.

The basic characteristics of the observed dynamics of the cells in the low to medium density area can be described follows:

  1. At any time, cells that are moving may stop moving. Conversely, at any time, stationary cells may start moving. This is not clear from the tracking example in Fig. 3, but has been observed in movies.
  2. There is a clear persistence in the direction of the movement. Cells that move in a given direction may continue moving in the same direction. This is apparent in Fig. 3 with the relatively straight tracks that the cells follow.
  3. The duration of the time in which a cell persists moving in a given direction is not fixed. It can vary over a large range of values. This is also evident from Fig. 3 as the cells travel different distances before turning.
  4. The direction of the motion of a cell may change at any time. While this is illustrated in Fig. 3, it is more obvious in the TLVM movies.
  5. When the direction of the motion changes, there is some tendency of moving in the direction of a neighboring cell. One example of this was highlighted in Fig. 3.

While analysis of time-lapse movies and snapshots of Synechocystis sp. is useful to characterize general movement of cells, it is still quite difficult to understand the fundamental interactions that are taking place between cells that allows for the observed movement and patterns.

3 A mathematical model for local interactions of Synechocystis sp

In this section we present mathematical models of the quasi-random motion of cells. Two models are discussed: a local-interactions model and a neighbor-dependent random motion model.

In exploring local interactions between cells, it is important to note that simple random motion does not elucidate the patterns we observe experimentally. This is evident, e.g., by the nature of the particle trajectories shown in Fig. 3.

3.1 Model assumptions

We begin by stating our formative assumptions.

  1. Neighbor detection. We assume that cells not only can detect the presence of neighbors within a given radius but also can detect the directions to neighbors within that radius. The radius of a neighbor detection will be considered to be a model parameter. This assumption is motivated by a few biological observations. First, the pili which grow on the cell are known to get longer with time. Cells have been observed to make connections through these pili. These connections are not easily observed with an optical microscope and are, hence, not shown in our images illustrating typical cell behavior in Fig 1. The length of the pili could be directly related to the neighbor detection distance. Second, a signaling molecule, which quickly degrades, could have a similar effect. A larger neighbor detection distance could be observed over longer time scales or for differing agarose densities. Lower densities may allows for faster diffusion of a signaling molecule. Such a signaling molecule has not yet been experimentally identified for this bacterium, but it is reasonable to assume one may exist. Overall, not that much is known about either mechanism biologically, and therefore we assume the existence of a finite neighbor detection distance without making any further assumptions about the underlying biological mechanism, leading to our next assumption.
  2. Communication mechanism unknown. As mentioned, we do not make any assumptions about the specifics of the communication mechanisms which allows cells to locate neighbors; though various mechanisms are feasible, including chemical signaling, local force sensing, and physical interaction. While these mechanisms are of interest, they will not play any role in our model, at least at this stage.
  3. Cell density. We assume that cell density is relatively low and that spatial hindrance effects are negligible so that we can use a particle model.
  4. Movement. When it comes to the motion, we assume cells may follow one of the following rules:
    1. A cell may continue moving in the direction it was previously moving,
    2. A cell may stop moving, or
    3. A cell may orient itself to move in a different direction. For our local interaction model, we will assume that the new direction is set toward a neighboring cell. For the two random models discussed above, any direction is equally viable.
  5. Memory. We assume that cells remember their preferred direction of motion even when they stop moving. In formulating the mathematical model we assume that this information can be retrieved at any time. While this has not yet been shown for Synechocystis sp., the Myxococcus bacterium has been observed to leave a polarized trail of slime, detectable by the cell, which supports the idea of a cell having a long-lasting memory [37].
  6. No directional persistence with change of direction. We also observe that, unlike flagellar chemotactic cells, Synechocystis sp. have no visually apparent directionality. That is, the cells do not have an obvious head or tail which indicates the preferred directionality of the cell and whether the cell is moving forward or turning. We model this by assuming that if a cell turns, the turning angle is not biased based on the previous direction. Other mathematical models have included a bias in which agents, when turning, tend to choose new directions which are similar to the current forward-moving direction, e.g. [30]. That is, we assume our cells have no inherent directionality as we do not have experimental evidence to justify imposing such a bias, especially in the phase of cells at the base of a finger where quasi-random motion is prevalent. Note that this type of bias is different from considering global forcing on the cells due to a directed light source.

3.2 Model A: A local-interactions model

Given the three rules of motion, we formulate our model as a discrete-time dynamical particle system. Our approach resembles the method used to study persistence in chemotaxis in [30]. The state of the system at time t is {xi(t),υi(t),θi(t)}i=1N where N is the number of cells, υi(t) [set membership] {0, 1} is the velocity of cell i, θi(t) [set membership] [0, 2π) is the angular direction of cell i, and xi(t) [set membership] R2 is the position of cell i calculated using υi(t), θi(t) and xi(t − 1). Note that the velocity of each cell can be either 1 or 0. In this way, we are assuming that a cell can either move at constant velocity or stop, as discussed in the assumptions above.

In this model, we allow cells to move in the direction of a neighboring cell. We define a neighboring cell to be any cell j which is within an interaction distance D from the cell i of interest. In this stochastic model, we only allow a cell to change direction with probability 1 − a. We formulate this as:

θi(t+1)={θi(t),with probability a,ηj,with probability 1ani for js.t.xjBD(xi),

where BD(xi) is a ball of radius D centered at xi, ni is the number of cells with position xj contained in BD(xi), i.e. ‖xi(t) − xj(t)‖2 < D, and ηj is the angle of the vector pointing from particle i to particle j. In this way, when a cell changes direction, with probability 1 − a, the direction presented by each neighboring cell j has an an equal probability, (1 − a)/ni, of being selected as the new directional heading for cell i. Note that if ni = 0, a cell does not have any neighbors. In this case the variable θi remains unchanged, i.e., θi(t + 1) = θi(t).

The velocity in our model (which is either 0 or 1) is determined independently of the choice of direction. It remains unchanged with probability b. It switches to the other state with probability 1 − b. This is formulated as

υi(t+1)={υi(t),with probability b,1υi(t),with probability 1b.

Finally, after updating the velocity and angle associated with each cell, all cell positions are updated according to


3.3 Model B: Random neighbor-dependent movement model

In the random neighbor-dependent movement model, a cell is only allowed to change direction if there is a neighbor within an interaction distance D. In this case, the direction is chosen randomly. We formulate this as:

θi(t+1)={θi(t),with probability a,η(t),with probability 1a if there exists xjBD(xi).

Here, for any t, η(t) is a uniform random variable distributed in [0, 2π). The position and velocity are calculated in the same way as in Section 3.2.

Before we consider simulations, note that these agent-based models impacted our ability to perform a rigorous analysis of the system or a thorough parameter space exploration. This is a known limitation of agent-based models compared to continuum models; however, we still are able to obtain results by looking at properties of simulations of this model. Further note that there are other simple modeling alternatives that could be considered using agent-based models; for instance, one could consider a random motion model with particles binding to the boundary. In this paper to focus our study on particle interactions, we only consider random motion model R and the neighbor interaction models A and B.

4 Simulations

To justify the complexity of the Models A and B, let us first consider a simpler model. Consider a discrete-time model of random motion where cells are allowed to change direction to a uniformly chosen random variable on [0, 2π) with probability 0.2 for each time step; we refer to this as Model R. We compare this random motion Model R with a the neighbor-dependent random model, Model B, where, for each time step with probability 0.2, cells are allowed to change direction to an angle chosen uniformly from [0, 2π), but only if there is a neighboring cell within 5 pixels distance. In Figure 5, we compare the random model, shown in Fig. 5(a)–(c), to Model B, where a directional change in movement depends on the presence of a neighbor within D = 5 pixels. Clearly, Model B, which is a rather simple neighbor-dependent model, yields results that better match the experimental data, compared with the random model. This is manifested in the cells in Fig. 5(d)–(f) that tend to align along the boundary (similar to what is experimentally observed).

Figure 5
(a)–(c) Random motion model (Model R) and (d)–(f) Random motion with neighbor dependence model (Model B) on changing direction for N = 250 cells at times t=0, 350, and 1000 time steps. The red circle around the cell in the upper left hand ...

We quantify this behavior in Figure 6. To do this, we partitioned the 100 pixel radius circle into 169 equal area sections. The sections are labeled in Fig. 6(a). For each model, we then identified the numbered partition for each cell in space and times from start to 1000 time steps. This data was then compiled into a histogram counting how many cells fall into each section for each model. In this way, each cell is counted 1001 times. In Fig. 6(b), we show a histogram of the cell locations for the random motion Model R. We observe that the cells are relatively uniformly distributed in space-time. In contrast, in Fig. 6(c), we show a histogram of the cell locations over each time step for the neighbor-dependent random motion, Model B. In this case, we observe that cells are gathering on the boundary, partitions 122 through 169, with approximately twice the probability of them falling in the interior of the circle, partitions 1 through 121. This ratio is likely parameter dependent. Qualitatively, this agrees with the experimental observations.

Figure 6
(a) Partitioning the domain into 169 equal area sections. The numbers represent partition number labels. (b) A histogram of the number of cells in each partition over time for the random motion model (Model R). (c) A histogram of the number of cells in ...

While this analysis helps to further understand the dot patterns presented in Fig. 5, we do not perform such a study on all of the simulations. The example we provide does illustrate the difference between random motion and the patterns of the motions we study. This pattern of motion has not been previously studied.

We now consider a few simulations comparing Model A, the local-interactions model, from section 3.2, to Model B, the neighbor-dependent random movement model, from section 3.3.

In Fig. 7, we see a comparison of Models A and B in which cells are constrained to a circle with a diameter of 200. We compare the effect of varying interaction distance and consider the cases of D = 5 and D = 10. The length scale of interaction distance D is shown in the upper left corner of each picture around a stationary cell. This cell, contained in a circle of radius D, is not included in any of the simulations. Each simulation is run with N = 250 cells and can be seen at t = 0, 50, and 350 increments. The initial conditions, seen in 7(a), 7(d), 7(g), and 7(j), are identical in all simulations. We observe that for the local interaction model in 7(c) and the random neighbor-dependent motion model in 7(i), the same cells tend to stay around the border, similar to what is seen in experiments. We recall that there is no mechanism that is embedded into the model that forces cells to stay on the boundary of the domain once the boundary is reached. Indeed, occasionally, individual particles return back into the domain. The cells appear to be slightly more uniformly distributed in 7(i) than in 7(c). In 7(f), Model A with D = 10 yields cells aggregating in small clumps of 5 to 15 cells. The simulated cells in Fig. 7(l) appear to be distributed uniformly with little cohesion to the boundary. In fact, the motion observed in the simulation depicted in 7(j)–(l) is very similar to the random motion model from Figure 5 because the chosen surface density and interaction distance are such that, on average, each cell has 1.5 neighbors, effectively removing the neighbor-dependence from this model.

Figure 7
(a)–(f) Model A; (g)–(l) Model B. N = 250 particles are constrained to a circle of diameter 200. Shown are simulation results at times t = 0, 50, 350. The interaction distance is D = 5 for (a)–(c), (g)–(i) and D = 10 for ...

In Figure 6, we offer a quantitative comparison of Models A and B for D = 5. We partition the domain into 169 equal-area sections and count the total number of particles that are present in each cell over 1000 time steps. We observe that for Model A, in Fig. 6(d), the boundary, partitions 122 through 169, has more cells on it with larger variability than for Model B (shown in Fig 6(c)). We suspect that this is due to the formation of aggregations on the boundary. As previously stated, cells have been experimentally observed to collect along the boundary.

The boundary conditions in Fig. 7 simulate a setup that is similar to the experimental setup focusing on short-time dynamics of cells. In experiments lasting relatively short periods of time (on the order of hours), cells are confined to a given area due to the properties of the underlying surface. However, in experiments lasting longer periods of time (on the order of days), cells overcome the underlying surface and, in the absence of a directional global light source, spread out on the surface of the agarose. In contrast, we consider large-time simulations of cells that are free to move anywhere in R2. Such simulations are shown in Fig. 8. The initial conditions are identical to the initial circular drop of cells seen in Fig. 8(a). The local-interactions model, Model A, is shown in Fig. 8(a) with D = 5 and Fig. 8(b) with D = 10. The random neighbor-dependent model, Model B, is shown in Fig. 8(c) with D = 5 and Fig. 8(d) with D = 10. All results are shown at t = 1000. Clearly, in Figs. 8(a),(c),(d), the majority of cells have dispersed. Fig. 8(b) stands different. In this case, of Model A, with the longer interaction distance, few aggregations still remain after these rather long-time simulations, in spite of the particles not being confined to the domain that is shown in the figure.

Figure 8
(a), (b) Model A; (c), (d) Model B. Both models are simulated in R2 at t = 1000 for N = 250 particles. The interaction distance is D = 5 for (a), (c) and D = 10 for (b), (d). The initial conditions and parameters a, b are identical to those used ...

To complement the previous cases, we also consider periodic boundary conditions. In Fig. 9, we show results obtained from simulations of cells on a 240 × 240 region with periodic boundary conditions. The same initial conditions as in Fig. 7(a) are used. The results are shown at time t = 3000. Model A is shown in Fig. 9(a) with D = 5 and in Fig. 9(b) with D = 10. Model B is shown in Fig. 9(c) with D = 5 and in Fig. 9(d) with D = 10. In these snapshots, aggregations are visible in Fig. 9(b) for D = 10 in Model A. There appears to be slightly more white space in Fig. 9(a) than in Fig. 9(c) and Fig. 9(d). In order to quantitatively compare this clumping effect seen in Fig. 9(a) and Fig. 9(c), we show a histogram of how many neighbors are found within a distance of 5 or 10 pixels of each cell in Fig. 10. For Model A with D = 5, we observe in Fig. 10(a) that within 5 pixels of each cell, there are more cells with two or more neighbors than in Model B shown in Fig. 10(c). In comparing the same histograms, there are also slightly more cells with no neighbors for Model B than for Model A. Within 10 pixels, the phenomenon is slightly reversed; that is, within 10 pixels, in Fig. 10(b) there are more cells with no neighbors for Model A than for Model B, in Fig. 10(d). This trend reversal implies that the aggregates seen in simulations of Model A may be tighter, as there are more isolated cells, meaning the other cells must be packed into a slightly smaller space.

Figure 9
(a), (b) Model A; (c), (d) Model B. Periodic boundary conditions; t = 1000; N = 250 particles. The interaction distance is D = 5 for (a), (c) and D = 10 for (b), (d). The initial conditions and parameters a, b are identical to those in Fig. 7.
Figure 10
Histogram illustrating how many cells have 0 to 10 neighbors within 5 or 10 pixels. Subfigures (a), (b) correspond to the simulation of Model A with D = 5 in Fig. 9(a). Subfigures (c), (d) correspond to the simulation of Model B with D = 5 in Fig. 9(c) ...

In Fig. 11, the positions of five cells, with identical initial conditions and the same model and interaction distance conditions as in Fig. 9, are illustrated for t = 1000. These trajectories are plotted on top of the final locations of the other 245 cells at time t = 1000. Note that the distance travelled by the particles is longer in both models when the local interaction distance is smaller.

Figure 11
The same five cells traced over 1000 time steps for (a), (b) Model A, and (c), (d) Model B. Periodic boundary conditions; t = 1000; N = 250 particles. The interaction distance is D = 5 for (a), (c) and D = 10 for (b), (d). The initial conditions and parameters ...

We further explore properties of Model A by considering multiple runs with the same initial conditions and setup. In Fig. 12, we consider six different runs, all with an interaction distance D = 10 pixels and N = 250 cells at t = 1000 time steps. Note that in comparing the six images the resulting distribution of cells and the size and number of aggregations are similar, but there do not appear to be any stable aggregations which routinely form in the same location. In fact, the aggregations which form move over time as well as dissociate and form new aggregations.

Figure 12
(a)–(f) Six realizations of Model A with the same initial conditions and parameters a and b, using D = 10 and N = 250 cells at times t = 1000 time steps. While precise locations of aggregations vary, the aggregations sizes and number are similar. ...

Varying the interaction distance D, we can gain further insight about the properties of this model. In Fig. 13, we consider three different interaction distances D = 10, 20, and 40 at t = 1000 time steps with the same initial conditions but different boundary conditions. In Fig. 13(a)–(c), cells are constrained to a circle, whereas in Fig 13(d)–(f), cells are simply restricted to R2 (note that in these cases not all cells are necessarily visible in the viewing window). In comparing the results for D = 10 and D = 20, we observe that the distribution and number and size of aggregations of cells is very similar. For D = 40, even the locations of the aggregations is almost identical. This is typical for simulations where D = 40; however, it is also possible for three aggregations of cells to form. The positions of the four aggregations and the distances between aggregations are generally slightly different. In both cases, all cells for D = 40 typically end up in an aggregation.

Figure 13
A study on the effect of interaction distance D. We consider (a), (d) D = 10, (b), (e) D = 20, and (c), (f) D = 40 with the same initial conditions and same parameter a, b and N = 250 cells at time t=1000. For (a)–(c), the cells are constrained ...

5 Conclusions

In this paper, we sought to answer a series of four questions regarding the seemingly random motion of and communication between cyanobacteria. In our attempt to answer these questions, we presented two mathematical models for studying the local interactions during the motion of Synechocystis sp. One model (B) assumed that a cell chooses to move in random directions and another model (A) assumed that a cell may choose to move in the direction of ones of its neighbors. Comparing these models, we can discuss question 1. While we were not able to completely isolate whether the chosen direction of cellular movement is completely random or not, in model A we do observe aggregations and interactions that appear to be more consistent with experiments. Additionally, our assumptions of model A offer a plausible mechanism by which aggregations form. To approach questions 2 and 3, for both models, we assumed that cells only move if a cell has a neighbor within some interaction distance D. We briefly compared this to a model of completely random motion as well as varied the interaction distance D to study the effect of such a distance. The movement observed with models A and B was much more consistent with experimentally observed cell movement than with the purely random model. However, it is still possible that the interaction distance D should instead be modeled as an independent property of each cell, possibly even varying with time.

Beyond these questions, there are similar characteristics of the model simulations and experimentally observed cell movement. True to the experimental observations, both models never reach a stationary steady state. The cells, unless stuck to the boundary with no neighbors within interaction distance D, are in a constant state of motion. Cells which do not permanently adhere to the boundary are capable of moving again at a later time. Both models include the persistence in the motion, yet, it is possible that the direction of motion will change at any time.

While the intercellular signaling in Synechocystis sp. is not well understood, it is thought to be present. The local interaction model we have developed does not specify how the signaling takes place. It is possible that signaling depends on the proximity or connectedness of pili of individual cells. Additionally, cells might release compounds which signal to other cells which direction to move. For instance, it has been shown that cAMP may be part of the signaling pathway required for motility [5]. The interaction distance in our model could be related to the time scale for diffusion of such a molecule, or possibly be related to the length of a pilus. If the latter is true, our model suggests that cells which have been genetically modified to have longer pili might form larger aggregates of cells, a possible biological implication to answer question 3.

Both models, the local-interactions model and the neighbor-dependent random model, produce results that better capture the experimental results than a model of a simple random motion. These models provides a possible explanation of the observed quasi-random behavior. The interaction distances for simulations considered in this paper seem to provide upper and lower bounds for reasonable behavior. Further study of this parameter, in connection with density and experimental results, may provide more insight into the relevant communication length-scale between cells. Understanding the local interactions between cells may turn out to be the key ingredient in understanding the mechanisms that control phototaxis. This might ultimately lead to understanding the source of the observed dynamical patterns (such as the finger formation, delayed motion, etc.) and we have carried out a preliminary study of the integration of local-interaction models with global forcing (due to the light source) [14].


  • [check] Phototactic cyanobacteria have quasi-random motion, not uniform motion toward light.
  • [check] Experimental data is analyzed to illustrate this unexplained motion.
  • [check] We propose two local interaction models accounting for quasi-random motion.
  • [check] Simulations suggest interaction distance between cells affects motion and patterns.


This work was supported in part by the joint NSF/NIGMS program under Grant Number DMS-0758374. The work of AG and DL was supported in part by Grant Number R01CA130817 from the National Cancer Institute. SW and DB were partly supported by the Carnegie Institution of Science. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Cancer Institute or the National Institutes of Health.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


1. Alber Mark S, Kiskowski Maria A, Glazier James A, Jiang Yi. On Cellular Au-tomaton Approaches to Modeling Biological Cells. In: Arnold DN, Santosa F, editors. Mathematical Systems Theory in Biology, Communication, and Finance. New York: Springer-Verlag; 2002. pp. 1–40. IMA 134.
2. Aoki I. A simulation study on the schooling mechanism in fish. Bulletin of the Japanese Society of Scientific Fisheries. 48(8):1081–1088.
3. Bhaya Devaki, Takahashi Akiko, Shahi Payam, Grossman Arthur R. Novel Motility Mutants of Synechocystis Strain PCC 6803 Generated by In Vitro Transposon Mutagenesis. Journal of Bacteriology. 2001;183:6140–6143. [PMC free article] [PubMed]
4. Bhaya Devaki. Light matters: phototaxis and signal transduction in unicellular cyanobacteria. Molecular Microbiology. 2004;53:745–754. [PubMed]
5. Bhaya Devaki, Naksugi Kenlee, Fazeli Fariba, Burriesci Matthew S. Phototaxis and Impaired Motility in Adenylyl Cyclase and Cyclase Receptor Protein Mutants of Synechocystis sp. Strain PCC 6803. Journal of Bacteriology. 2006;188(20):7306–7310. [PMC free article] [PubMed]
6. Bhaya Devaki, Levy Doron, Requeijo Tiago. Group dynamics of phototaxis: Interacting stochastic many-particle systems and their continuum limit. Hyperbolic Problems: Theory, Numerics, Applications. 2008:145–159.
7. Burkhart Uhland, Hader Donat-P. Phototactic attraction in light trap experiments: A mathematical model. Journal of Mathematical Biology. 1980;10:257–269.
8. Burriesci Matthew, Bhaya Devaki. Tracking phototactic responses and modeling motility of Synechnocystis sp. strain PCC6803. Journal of Photochemistry and Photobiology B: Biology. 2008;91:77–86. [PubMed]
9. Couzin Iain D, Krause Jens, James Richard, Ruxton Graeme D, Franks Nigel R. Collective memory and spatial sorting in animal groups. Journal of Theoretical Biology. 2002;218:1–11. [PubMed]
10. Cucker Felipe, Smale Steve. On the mathematics of emergence. Japanese Journal of Mathematics. 2007;2:197–227.
11. Cucker Felipe, Smale Steve. Emergent behavior in flocks. IEEE Trans. Automat. Control. 2007;52 852862.
12. Degond Pierre, Motsch Sébastien. Continuum limit of self-driven particles with orientation interaction. Mathematical Models and Methods in Applied Sciences. 2008;18:1193–1215.
13. Fatehi Hasnaa, Meyer-Hermann Michael, Figge Marc Thilo. Modelling cellular aggregation induced by chemotaxis and phototaxis. Mathematical Medicine and Biology. 2010;27:373–384. [PubMed]
14. Galante Amanda, Wisen Susanne, Bhaya Devaki, Levy Doron. Stochastic models and simulations of phototaxis. In: Sayama H, Minai AA, Braha D, Bar-Yam Y, editors. Unifying Themes in Complex Systems Volume VIII: Proceedings of the Eighth International Conference on Complex Systems; New England Complex Systems Institute Series on Complexity: NECSI Knowledge Press; 2011. pp. 105–119.
15. Ha Seung-Yeal, Tadmor Eitan. From particle to kinetic and hydrodynamic descriptions of flocking. Kinetic and Related Models. 2008;1(3):415–435.
16. Ha Seung-Yeal, Lee Kiseop, Levy Doron. Emergence of time-asymptotic flocking in a stochastic cucker-smale system. Communications in Mathematical Sciences. 2009;7:453–469.
17. Hillen Thomas, Painter Kevin J. A users guide to PDE models for chemotaxis. Journal of Mathematical Biology. 2009;57:183–217. [PubMed]
18. Horstmann D. From 1970 until present: the Keller-Segel model in chemotaxis and its consequences I. Jahresberichte DMV. 2003;105(3):103–165.
19. Huth Andreas, Wissel Christian. The simulation of movement of fish schools. Journal of Theoretical Biology. 1992;156(3):365–385.
20. Ikeuchi Masahiko, Ishizuka Takami. Cyanobacteriochromes: a new superfamily of tetrapyrrole-binding photoreceptors in cyanobacteria. Photochemical & Photobiological Sciences. 2008;7:1159–1167. [PubMed]
21. Kaiser Dale. Myxococcus—from single-cell polarity to complex multicellular patterns. Annual Review of Genetics. 2008;42:109–130. [PubMed]
22. Katz Yael, Tunstrom Kolbjorn, Ioannou Christos C, Huepe Cristian, Couzin Iain D. Inferring the structure and dynamics of interactions in schooling fish. Proceedings of the National Academy of Sciences. 2011;108(46):18720–18725. [PubMed]
23. Keller Evelyn F, Segel Lee A. Model for chemotaxis. Journal of Theoretical Biology. 1971;30(2):225–234. [PubMed]
24. Levy Doron, Requeijo Tiago. Modeling group dynamics of phototaxis: from particle systems to PDEs. Discrete and Continuous Dynamical Systems - Series B. 2008;9:103–128.
25. Levy Doron, Requeijo Tiago. Stochastic models for phototaxis. Bulletin of Mathematical Biology. 2008;70:1684–1706. [PubMed]
26. Levy Doron, Ha Seung-Yeal. Particle, kinetic and fluid models for phototaxis. Discrete and Continuous Dynamical Systems - Series B. 2009;12:77–108.
27. Li Yue-Xuan, Lukeman Ryan, Edelstein-Keshet Leah. Minimal mechanisms for school formation in self-propelle particles. Physica D. 2008;237:699–720.
28. Lukeman Ryan, Li Yue-Xian, Edelstein-Keshet Leah. Inferring individual rules from collective behavior. Proceedings of the National Academy of Sciences. 2010;107:12576–12580. [PubMed]
29. Maree Athanasius F. M, Panfilov Alexander V, Hogeweg Paulien. Phototaxis during the slug stage of Dictyostelium discoideum: a model study. Proceedings of the Royal Society B. 1999;266:1351–1360.
30. Nicolau Dan V, Armitage Judith P, Maini Philip K. Directional persistence and the optimality of run-and-tumble chemotaxis. Computational Biology and Chemistry. 2009;33:269–274. [PubMed]
31. Othmer Hans G, Hillen Thomas. The diffusion limit of transport equations II: Chemotaxis equations. SIAM Journal of Applied Mathematics. 2002;62(4):1222–1250.
32. Parrish Julia K, Viscido Steven V, Grunbaum Daniel. Self-organized fish schools: An examination of emergent properties. Biological Bulletin. 2002;202:296–305. [PubMed]
33. Patlak Clifford S. Random walk with persistence and external bias. Bulletin of Mathematical Biophysics. 1953;15:311–338.
34. Tindall Marcus J, Maini Philip K, Porter Steven L, Armitage Judith P. Overview of mathematical approaches used to model bacterial chemotaxis II: bacterial populations. Bulletin of Mathematical Biology. 2008;70:1570–1607. [PubMed]
35. Sbalzarini IF, Koumoutsakos P. Feature point tracking and trajectory analysis for video imaging in cell biology. Journal of Structural Biology. 2005;151(2):182–195. [PubMed]
36. Vicsek Tamás, Czirók András, Ben-Jacob Eshel, Cohen Inon, Shochet Ofer. Novel type of phase transition in a system of self-driven particles. Physical Review Letters. 1995;6:1226–1229. [PubMed]
37. Yu Rosa, Kaiser Dale. Gliding motility and polarized slime secretion. Molecular Microbiology. 2007;63(2):454–467. [PubMed]