Cellular automaton models for time-correlated random walks: derivation and analysis

Cellular automaton models for time-correlated random walks: derivation and analysis

J. M. Nava-Sedeño Technische Universität Dresden, Center for Information Services and High Performance Computing, Nöthnitzer Straße 46, 01062 Dresden, Germany H. Hatzikirou R. Klages School of Mathematical Sciences, Queen Mary University of London, Mile End Road, London E1 4NS, United Kingdom A. Deutsch Technische Universität Dresden, Center for Information Services and High Performance Computing, Nöthnitzer Straße 46, 01062 Dresden, Germany
Abstract

Many diffusion processes in nature and society were found to be anomalous, in the sense of being fundamentally different from conventional Brownian motion. An important example is the migration of biological cells, which exhibits non-trivial temporal decay of velocity autocorrelation functions. This means that the corresponding dynamics is characterized by memory effects that slowly decay in time. Motivated by this we construct non-Markovian lattice-gas cellular automata models for moving agents with memory. For this purpose the reorientation probabilities are derived from velocity autocorrelation functions that are given a priori; in that respect our approach is “data-driven”. Particular examples we consider are velocity correlations that decay exponentially or as power laws, where the latter functions generate anomalous diffusion. The computational efficiency of cellular automata combined with our analytical results paves the way to explore the relevance of memory and anomalous diffusion for the dynamics of interacting cell populations, like confluent cell monolayers and cell clustering.

Introduction

Within the past two decades transport processes in many branches of the sciences were observed to be anomalous, in the sense that they do not obey the laws of conventional statistical physics and thermodynamics [1, 2, 3, 4, 5, 6, 7, 8]. Important cases are diffusion processes where the long-time mean square displacement (MSD) does not grow linearly in time. That is, , where the angular brackets denote an ensemble average, does not increase with as expected for Brownian motion but either subdiffusively with or superdiffusively with [9, 10, 11]. After pioneering work on amorphous semiconductors [12], more recently anomalous diffusion has been detected in many other complex systems [3, 4, 6, 8]; here well-known examples of physcial systems are nanopores [13], plasmas [14] and glassy material [15].

Biological systems frequently exhibit anomalous properties as well: Prominent examples are the foraging of organisms [16], epidemic spreading [17] and the diffusion of macromolecules in biological cells [6]. Especially, it was found that many types of cells migrate anomalously: Hydra cells [18], mammary gland epithelial cells [19], MDCKF cells [20], amoeboid Dictyostelium cells [21, 22], T cells [23], breast carcinoma cells [24] and stem cells [25] were all experimentally observed to move superdiffusively, typically with non-Gaussian position and/or velocity distribution functions [18, 19, 20, 21, 22, 23, 24] accompanied by either exponential or non-exponential position [23], and exponential [26, 21, 22] or power law [18, 19, 20] velocity autocorrelation function (VACF) decay. For T cells it was argued that superdiffusion optimizes their search to kill intruding pathogens [23, 27]. While all these results are on single cell migration, currently collective cell migration is moving into the center of interest [28], where cells interact with each other, e.g., by chemical signalling [29]. Interesting phase transitions inside dense tissues of epithelial cell monolayers were reported [30] and partially traced back to particular features of single-cell migration [31]. It was also observed experimentally that superdiffusion appears to foster the formation of clusters of stem cells leading to tissue formation [25]. Other works investigate the role of interacting agents for phase transitions in active matter [32], and collective anomalous dynamics emerging from the interaction of single agents [33, 34].

On the theoretical side there are many different ways to model anomalous diffusion in terms of stochastic processes, such as continuous time random walks (CTRW) [1, 5], generalized Langevin equations [2], Lévy flights and walks [8], fractional diffusion equations [1, 4, 5], scaled Brownian motion and heterogeneous diffusion processes [7]. A subset of these models, most notably generalized Fokker-Planck equations [18, 19, 20], generalized Langevin equations [21, 22] and generalized random walks [23, 24], has been used to model anomalous movement of single biological cells. However, solving equations for anomalously moving single particles analytically or numerically is typically difficult already [1, 2, 4, 5, 7, 8]. To our knowledge there exists no systematic attempt to generalize this theory to model interacting many-particle systems; the only exception we are aware of is a line of work in plasma physics [14].

On the other hand, several models have been introduced to study the collective movement of particles and cells [35, 36, 37] . Cellular automata (CA) in particular have the advantage of being less computationally demanding than continuous models when performing simulations. A specific type of CA is the so-called lattice-gas cellular automata (LGCA) [38]. In LGCA, each lattice node can contain several particles, which at each time step are rearranged within the lattice node according to the interaction rule, and subsequently moved to a neighboring node. In a biological context particles can be regarded as cells, while the LGCA rules mimic cell migration and interaction. Furthermore, LGCA have proved to be amenable to mathematical analysis [39]. For this reason, LGCA have been introduced as mesoscopic models for single and collective cell migration [40, 41, 42, 43]. So far, none of the mentioned models has considered anomalous migration of single cells.

It thus arises the need to design simple fundamental schemes by which the collective properties of interacting agents can be studied whose individual dynamics is anomalous. Using our methods will enable to explore the relevance of microscopic single-particle dynamics for emerging collective phenomena. We thus devise a scheme by which anomalous dynamics of many interacting agents can be simulated efficiently, which is based on capturing the non-trivial decay of VACFs. This approach generates superdiffusion if the correlation decay is of power law-type [18, 19, 20]. We emphasize that our data-driven approach can be applied to any moving entity that exhibits dynamics with non-trivial correlation decay, a feature that may be expected to hold more generally for the movement of biological organisms [16].

We use the LGCA modeling framework and construct various time-correlated random walk models. After briefly introducing the LGCA concept, we define an LGCA model for unbiased random walk. Next, motivated by the biophysical mechanism of single cell crawling we construct a persistent random walk LGCA model wherein angular (orientation) correlations give rise to temporal correlations. Subsequently, we construct a first LGCA model for time-correlated random walk which is data-driven, as the model’s reorientation probabilities are derived by assuming that the exact temporal dependence of the VACF is known a priori. Finally, we develop a generalized time-correlated random walk LGCA model for cell movement at short and medium time regimes by curing a deficiency of our first time-correlated random walk model for short times. Figure 1 shows single cell tracks with the corresponding VACFs and MSDs for our main two classes of LGCA models we are dealing with, which are Markovian and non-Markovian random walks, exemplified by showing their basic features.

Figure 1: Basic types of diffusive movement in two dimensions characterized by two key quantities (VACF and MSD). Shown in the first column are the tracks of a single particle starting at that exhibits either long-time normal diffusion (top row) or superdiffusion (bottom row). The color gradient changes from blue to yellow with elapsed time. The second column displays the particle’s velocity autocorrelation function (VACF) Eq. (3), the third column its mean square displacement (MSD) Eq. (4).

Lattice-gas cellular automata

Cellular automata are mathematical models where the states of discrete lattice nodes are updated at discrete time steps. If the states of the lattice sites are Boolean, such states can be interpreted as presence/abscence of a particle at a particular node. The lattice-gas cellular automaton is a specific CA type, which has two important characteristics: first, particle reorientation and migration are separated into a probabilistic and a deterministic step, respectively. Secondly, to each node, velocity channels are associated which can be occupied by at most one particle (exclusion principle). The set of velocity channels is given by , (see Fig 2). Particles move in discrete time steps of duration to neighboring nodes located a distance away in the lattice. At each time step particles adopt the orientation with a probability , called the reorientation probability, where , is the elapsed simulation time. Subsequently, cells will be deterministically translocated to the nearest neighbor located in the direction of ; see Fig. 2 for a sketch of LGCA dynamics.

Figure 2: LGCA dynamics. At each time step a particle is assigned an orientation with a probability . Subsequently, the particle is translocated to the nearest neighbor in the direction of its orientation.

If we were to have different particles in a single node, then the probability of the particles adopting the orientations would be given by

(1)

Classical random walk

Here we define an LGCA model for unbiased random walk[44]. In this model, at all timesteps, all orientations are chosen with equal probability [45, 46]. This means that the reorientation probability is given by

(2)

In order to characterize the movement of a particle in this model we calculate time-dependent expressions for the VACF and the MSD, which measure the persistence, in terms of memory decay in time, and the spatial exploratory power of a moving particle, respectively. In LGCA models space, time and particle velocities are discrete so that the VACF is given by [47]

(3)

where is the orientation of the particle at the -th time step. The MSD is calculated by

(4)

where is the norm of the particle displacement at time step defined as , where is the position of the particle at time step . The probability can be calculated from Eq. (2) by noticing that . In this simple random walk model Eq. (3) reduces to a sum of cosines over homogeneously distributed angles. Hence the VACF is given by , where is the Kronecker delta. In the limit the VACF tends to

(5)

where is the Dirac delta function. Eq. (5) means that the movement of the particle is uncorrelated as soon as it starts moving, i.e. the orientation of the particle at any time step is completely independent from its previous orientation, which is the Markov property. On the other hand, simple combinatorics can be used to calculate the particle’s MSD yielding . We can rewrite this expression by using the general definition of the diffusion coefficient

(6)

where is the dimension of space. For a memoryless random walk this equation boils down to with an MSD of . Given that is the time step length and that is the elapsed time, the MSD is

(7)

Eq. (7) shows that the classical random walk model trivially yields a normal diffusion process, where the MSD increases linearly in time [48].

Persistent random walk

The assumption of the classical random walk in Eq. (2) that all the directions of movement are equally probable is generally not true. We will now construct an LGCA model motivated by a simple biophysical model for a single, persistently moving cell.

Rule derivation

Figure 3: Reorientation of a biological cell moving persistently. A cell which moves in a direction feels an intracellular force and reorients towards the direction of the force due to the torque .

Biological cells move by exerting forces to propel themselves. In the case of eukaryotic cells such as fibroblasts, movement is achieved by crawling over the substrate. Crawling is performed by polymerization of the actin cytoskeleton at the leading edge propelling the cell in this direction. We can identify the direction of the actin concentration gradient with the direction of movement of a cell . Furthermore, intracellular forces due to actin activity point towards the direction of new actin polymerization. The cell then rotates/reorients/repolarizes due to the torque defined by , whose norm is given by

(8)

where , , , and and are the angular components of and , respectively; see Fig. 3. In the overdamped regime, characteristic of the cellular environment, the intracellular force will not cause the cell to rotate indefinitely but rather will cause the cell to rotate until and are parallel. Taking this into account, it is possible to rewrite Eq. (8) as , where is the direction of motion of the cell after it has finished its reorientation.

The torque is given in terms of an energy of rotation as . Using this equation, the energy of rotation is then given by

(9)

where is the amplitude of the torque generated inside the cell. Having defined the energy of rotation, Eq. (9), we can describe the cell’s reorientation by a Langevin equation [49] as with relaxation constant and a zero-mean, delta correlated noise term such that , where is the rotational diffusion coefficient. Based on this we can immediately derive the LGCA reorientation probabilities [50].

Figure 4: Comparison between random walk and persistent random walk models. Shown are theoretical results for VACF (top) and MSD (bottom) in the random walk (dashed maroon line) and persistent (solid blue line) models. The parameter values are , , and .

These probabilities for a single cell then read

(10)

where is the orientation of the cell at the previous time step, is the normalization constant (also known as the partition function) and is the sensitivity, where is the effective relaxation constant.

Model analysis and results

Using Eq. (10) we can calculate the VACF and MSD for this model. By using the properties of the partition function in Eq. (10) the VACF at every time step is (see Sec. B in the Supplementary Information)

(11)

where the exponent depends heavily on the lattice dimension and geometry. In particular, in a 2D square lattice we have . In all geometries the exponent is over its domain (see again Sec. B in the Supplementary Information).

Equation (11) can be generalized to continuous time and space by employing the relations between the time and space scalings, namely the diffusion coefficient in the random walk limit and the instantaneous cell speed , where is the lattice spacing and the time step length. Taking the limit yields the VACF in continuous time and space

(12)

When time is discrete the MSD is given by [51, 52] (see Sec. A in the Supplementary Information) . Calculating the expected value on the right hand side we obtain (see again Sec. A in the Supplementary Information)

(13)

Using Eq. (11) and taking the limit of small time step length we get

(14)

which can be easily integrated to obtain the MSD for continuous time

(15)
Figure 5: Comparison between persistent, time-correlated and generalized time-correlated random walk models. Shown are simulation results (circles, mean standard error of the mean) and theoretical prediction (solid line) for VACF (top row) and MSD (bottom row). Parameters are and in all cases. i) Sensitivity values: (left), (right). ii) and iii) , , and exponents: (left), (right).

Equation (14) agrees with the formal solution of the MSD for an overdamped Langevin equation with colored noise [53]. Correspondingly, Eq. (10) coincides with a Langevin process where the noise is not white but colored whose correlation is given by Eq. (12). Using Eq. (15) we find that when rather quickly.

Comparing Eq. (5) to Eq. (12) in Fig. 4 we see that in this second model the velocities are no longer delta correlated but now decay exponentially in time. On the other hand, for long times both Eqs. (7) and (15) yield normal diffusion. However, for short times cells performing persistent random walks move superdiffusively, contrary to cells performing classical random walks; see again Fig. 4. In Fig. 5i, Eqs. (12) and (15) are compared with results from LGCA simulations, where we see that the theory adequately predicts the observed simulation results. Details on the computational implementation are found in Sec. H of the Supplementary Information.

Time-correlated random walk

Due to the exponential decay of correlations the previous model did not show superdiffusion at long time scales. It turns out that finding a homogeneous, isotropic Markovian model that shows superdiffusion and power-law decaying correlations is not possible (for a proof see Sec. C of the Supplementary Information).

Theorem 1.

The velocity autocorrelation function of a particle whose orientations are given by a homogeneous, symmetric Markov chain is either delta-correlated, i.e. , where is the Kronecker delta; alternating, i.e. , ; or exponentially decaying, i.e. , .

To reproduce superdiffusion and power law decaying autocorrelations, we will construct a non-homogeneous model by assuming that the time dependency of the VACF is a known power law.

Rule derivation

We now assume that the VACF is known. In particular, if the movement is power law-correlated the VACF has the form [53]

(16)

where and and assume that , to disregard the movement at short times, where Eq. 16 diverges. The rate of decay of the VACF is proportional to the exponent . The crossover time specifies the time a which . The walk is positively correlated if , and anti-correlated if . Because the process is non-homogeneous, in Eq. (3) explicitly depends on the velocity channels and the current time step . Combining Eqs. (3) and (16) we obtain the following relation [47]:

(17)

It is possible to derive the reorientation probabilities by expanding Eq. 17 for every time step. Additionally, in order to reduce the number of equations, we make the following assumptions:

  • The reorientation probabilities are independent, that is, the probability of following a certain trajectory is .

  • There is symmetry around the initial orientation, i.e. if then , .

Figure 6: Comparison between persistent and time-correlated random walk models. Shown are theoretical results for VACF (top) and MSD (bottom) in the persistent (solid blue line) and time-correlated (dotted green line) models. The parameter values are , , , , , and .

Using these assumptions we can derive the general expression for the reorientation probabilities determining a certain VACF (see Sec. D in the Supplementary Information)

(18)

where is the dimension of space and is the number of lattice directions given by the lattice geometry. If the VACF follows a power law, then Eq. (18) is always valid if the crossover time is smaller than the time step length, as the divergence of Eq. (16) is avoided. If , we assume that the movement at short times is completely correlated, i.e. the VACF is given by

(19)

where is such that . We can then define a piecewise reorientation probability

(20)

where is the duration of ballistic motion.

Model analysis and results

For this model the VACF is obviously known as the reorientation probabilities were calculated specifically to reproduce it. We now calculate the MSD for this model. The reorientation probabilities given by Eq. (18) only depend on the initial cell orientation and on the time step , so they are independent of other orientations at other times. Because of this independence of orientations the MSD is given by (see Sect. A in the Supplementary Information)

(21)

yielding a Taylor-Green-Kubo formula [47, 54]. In the limit of small time step lengths the MSD is given by

(22)

Equation (22) shows that there is a correction to the random walk diffusion coefficient as well as a new term depending on the particle speed. If is given by Eq. (16), Eq. (22) can be integrated yielding

(23a)
(23b)
(23c)

These expressions for the MSD are valid when . For long crossover times when reorientation probabilities are given by Eq. (20), the MSD is (see Sec. E in the Supplementary Information):

(24)

Figure 5ii shows a comparison of Eqs. (16) and (22) with LGCA simulations. Details on the computational implementation are found in Sec. H of the Supplementary Information. From Eqs. (23) we see that in general , which defines three regimes:

  1. : superdiffusive regime, arising from the term

  2. : subdiffusive regime, as the term dominates at short times only

  3. : normal diffusive regime, as the linear term is the dominating term

Figure 6 shows that, while the VACF decays rapidly in the persistent model, in the time-correlated model the VACF decays much more slowly. Additionally, the movement is superdiffusive in both models are short times, however this behavior is long-lasting in the time-correlated model.

Generalized time-correlated random walk

The reorientation probabilities derived for the time-correlated random walk are only valid for certain time ranges, due to the divergence of the VACF when . We will now derive a generalized model which is valid on both short and long time scales. For this purpose we use what is called the maximum caliber formalism [55], which we introduce briefly.

Rule derivation

The maximum caliber formalism has been proven successful to derive models for dynamic systems from data. The procedure consists in maximizing the entropy over a path of system evolutions, with the constraint of reproducing certain observables. The procedure of entropy maximization does not only ensure that the resulting model contains as few assumptions as possible, but is also considered the only method of correctly obtaining unknown probability distributions from known data [55]. The procedure is as follows:

Let the path entropy, or caliber, be defined as where is a possible path followed by the system during its evolution. The probability of following such a path is given by . In the case of a single random walker, the path is the entire history of particle velocities up to the last time step . Furthermore, we constrain the unknown probabilities by a normalization constant and the observed VACF (in this case Eq. (17)). Then the problem translates into optimizing the functional [56] , where and are Lagrange multipliers to be determined. The Lagrange multiplier is given by (see Sec. F in the Supplementary Information), and determines the normalization constant.

Using the expression for we obtain the reorientation probability

(25)

where is the normalization constant for the reorientation probability.

If one required not only that the VACF was observed but also that the autocorrelation function would decay similarly independently of the start and end points, i.e. , if , the problem would be similar with the exception that we would now have constraints if the trajectory consists of time steps. Analogously, the probabilites would then be given by

(26)

Model analysis and results

The VACF can be easily calculated by using the properties of the partition function . Given a distribution , the expected value of the function is

(27)
Figure 7: Comparison between generalized and time-correlated random walk models. Shown are theoretical results for VACF (top) and MSD (bottom) in the time-correlated (dotted green line) and generalized time-correlated (solid red line) models. The parameter values are , , , , .

Combining Eqs. (3), (25) and (27) the VACF is given by , where in the case of of probabilities given by Eq. (26). In the case of a 2D square lattice the partition function can be easily calculated as . Taking the logarithm and differentiating we obtain . If we consider power law correlations and take the limit we obtain

(28)

Eqs. (16) and (28), corresponding to the time-correlated random walk and generalized time-correlated random walk, respectively, are visually compared in Fig. 7. A Taylor expansion of around (i.e. valid for or ) shows that up to second order terms , as expected. For the VACF decays as a power law as well (see Sec. G in the Supplementary Information), , and the difference between both functions behaves as . When the probabilites are given by Eq. (25) particle orientations are independent of one another, as it was the case in the time-correlated model, and the MSD is given by Eq. (21). Therefore, in the limit the MSD is

(29)

If we consider power law correlations and expand in Taylor series around up to second order, we recover Eq. (22) when is given by Eq. 16. Therefore the MSD in this model follows similar regimes as those of the time-correlated model Eq. (23). Eqs. (24) and (29) are visually compared in Fig. 7, and see Fig. 5iii for a comparison of Eqs. (28) and (29) with LGCA simulations. Details on the computational implementation are found in Sec. H of the Supplementary Information.

When the probabilities are given by Eq. (26) instead the MSD is given by Eq. (13), as all orientation pairs and are correlated, where the correlation is given by thus not depending specifically on the values of and but on their difference only. Therefore, for a 2D square lattice and a power-law decaying VACF in the limit the MSD is given by

(30)

Again we see that this expression agrees with the formal solution of the MSD for an overdamped Langevin equation with colored noise [53]. In this case, however, the noise correlation is not decaying exponentially. If we expand the hyperbolic tangent on the right hand side around (i.e.  or ) up to the second term and integrate we obtain the MSD as

(31a)
(31b)
(31c)

We conclude from Eqs. (31) that

When we have and the process is superdiffusive. When we have , at short times the term dominates, and the process is subdiffusive while at long times the linear term dominates yielding normal diffusion. Finally, when we have and the process is completely normal diffusive.

Summary and discussion

The goal of our study was to design a simple model for a single particle moving with memory in abscence of any environmental cue. We chose a cellular automaton, specifically, an LGCA because it is a flexible and computationally efficient framework and has the potential to analyze collective behavior in populations of moving particles or cells. After having introduced an LGCA model for unbiased random walk, we have derived three different novel time-correlated LGCA models for single particle migration.

The subsequent persistent random walk LGCA model was derived from a biophysically-motivated Langevin equation for particle reorientation in an overdamping environment. We showed that in this model the VACF decays exponentially. Furthermore, we proved that a particle in this model moves superdiffusively only at short times while it diffuses normally in the long time limit. This behavior as well as the expression we found for the MSD agree with that found by Othmer [57] using the telegrapher’s equation and also to the one by Chechkin et al. [53] for exponential noise.

The time-correlated random walk model was derived by assuming that the specific form of the VACF is known. We also assumed that reorientation probabilities were completely independent, and that particles have no preference in turning left or right. We considered the specific case of a power law-decaying VACF and showed that the MSD exhibits two transitions when and . For small exponents the particle moves superdiffusively on every time scale. At intermediate exponents there are non-linear contributions dominating at short times. For large exponents, all non-linear contributions vanish in the long time limit resulting in normal diffusion.

Finally, we derived a generalized LGCA model by maximizing the diffusing particle’s path entropy while retaining the constraint of reproducing a certain VACF. In this model the reorientation probability Eq. (25) is similar to the reorientation probability of the persistent random walk Eq. (10), with some differences: in Eq. (25) the particle’s orientation is compared to its initial orientation while in Eq. (10) the particle’s orientation is compared to the particle’s orientation at the previous time step. Furthermore, in Eq. (10) we have a constant parameter while in Eq. (25) this parameter decays with time, i.e. . We recall that Eq. (10) results from considering a Langevin equation for the particle’s reorientation. The parameter depends on the magnitude of the reorienting force, the relaxation constant (related to friction) and the angular diffusion constant . A time-dependent parameter would be obtained when considering a generalized Langevin equation resulting from either a time-dependent reorientation force or friction if these values changed much more slowly than the time needed for the particle to be displaced. Taking this into account, Eq. (10) describes the movement of a particle when reorientations can be performed almost instantaneously compared to the time required for the particle to move in space. Eq. (25) on the other hand describes movement when the particle keeps moving but needs considerably more time to change its initial orientation. When in addition the VACF is required to be invariant under time translations we showed that the corresponding MSD time regimes match to those found by Chechkin et al. [53] for power law-correlated noise.

We have verified our analytical results of all constructed models by comparing them to LGCA computer simulations. In order to derive the analytical VACF and MSD expressions, we have considered the limit . In this limit, the macroscopic time remains small even after several time steps. It stands to reason that, for , , the difference between our analytical expressions and simulations becomes negligible.

In their present form our LGCA models assume (i) the particle has constant instantaneous speed ; (ii) the particle moves to a neighboring site at every time step; and (iii) the particle moves on a regular lattice, which impacts the specific expression of the VACF. All these models could be extended by considering different instantaneous speeds, as well as waiting times between subsequent displacements, by using multispeed LGCA and adding rest (zero velocity) channels, respectively. Effects of the lattice regularity on the single particle movement can be compensated by choosing the sensitivity appropriately in the persistent random walk model as well as the crossover time in the generalized time-correlated model. In the time-correlated model the VACF does not depend on the lattice geometry.

Our new models could also be extended to account for external forces acting on the particle, independent from its intrinsic anomalous movement. For extending the models, we can consider that particle reorientations are caused by internal correlations of individual cells and by particle interactions. The probability of having particles in a node with a certain orientation due to internal particle orientations has already been introduced in Eq. (1), while individual particle reorientation probabilities would be given according to one of the models introduced in this work. On the other hand, the probability of having particles in a node with a certain orientation due to particle interactions would be a function of other particles’ positions and orientations (see [50], for examples of such probabilities). If we assume that both probabilities are independent, then the reorientation probability for all particles would simply be . Such an extension could be useful for studying physical systems such as plasma gases [14]. Furthermore, it would be interesting to construct LGCA models generating Lévy walks exhibiting non-Gaussian probability density functions [5, 8]. More importantly, and tracing back to the original motivation of our work, due to the computational efficiency of LGCAs, our schemes could be applied to model large groups of interacting cells to study the impact of persistence and time correlations in single-cell dynamics on collective phenomena. Highly promising examples are coordination and swarming in bacteria [34], pluripotent cells during early development [25], and the emergence of phase transitions in collective cell migration [28, 29, 30, 31]. Moreover, the non-cellular microenvironment is crucial for cell migration phenomena. Recently, the impact of complex environments on cell dissemination has been studied with a cellular automaton model [58]. It would be interesting to extend the models introduced here to analyze the impact of anomalous dynamics and complex microenvironments on cell dissemination and cancer invasion.
The LGCA modeling framework followed in this work is characterized by simplifying the concept of a moving particle to movement in discrete time steps between discrete nodes on a regular lattice, possessing only a finite, discrete set of velocities. On one hand, this “discrete approach” is decidedly more simplified and abstract than “continuous approaches” such as continuous time random walks or fractional diffusion equations. On the other hand, as we have shown in the present work, the LGCA offers an advantage not only in computational efficiency and straightforward multiparticle extension, but also in ease of model analysis. This gets rapidly complex in the aforementioned continuous approaches [1, 2, 3, 4, 5, 6, 7, 8] but remains feasible in the LGCA, even when dealing with systems of interacting particles [35, 38, 42].
In the era of “Big Data”, there is an abundance of biological data. Single or collective cell motility can be measured in vitro or in vivo via various experimental methods such as in vivo two-photon imaging [59] or cell cytometry [60], respectively. In this regard, there is a need for “data-driven” modeling frameworks. Our work comes timely to fulfill this scope by proposing the “data-driven” modeling of single particle superdiffusive behavior without prior knowledge of the mechanisms at work. Such an approach is vital for the study of phenomena whose driving mechanisms are currently unknown or challenging to model [61, 62, 63].

Acknowledgements

The authors thank the Centre for Information Services and High Performance Computing (ZIH) at TU Dresden for providing an excellent infrastructure. The authors acknowledge support by the German Research Foundation and the Open Access Publication Funds of the TU Dresden.The authors would like to thank Anja Voß-Böhme, Lutz Brusch, Fabian Rost, Osvaldo Chara, Simon Syga, and Oleksandr Ostrenko for their helpful comments and fruitful discussions. Andreas Deutsch is grateful to the Deutsche Krebshilfe for support. Andreas Deutsch is supported by the German Research Foundation (Deutsche Forschungsgemeinschaft) within the projects SFB-TR 79 “Materials for tissue regeneration within systemically altered bones” and Research Cluster of Excellence “Center for Advancing Electronics Dresden” (cfaed). Haralampos Hatzikirou would like to acknowledge the SYSMIFTA ERACoSysMed grant (031L0085B) for the financial support of this work and the German Federal Ministry of Education and Research within the Measures for the Establishment of Systems Medicine, project SYSIMIT (BMBF eMed project SYSIMIT, FKZ: 01ZX1308D). Josué Manik Nava-Sedeño is supported by the joint scolarship program DAAD-CONACYT-Regierungsstipendien (50017046) by the German Academic Exchange Service and the National Council on Science and Technology of Mexico.

Author contribution statement

All authors formulated the mathematical model. J.M.N.S. performed the analysis and simulations. All authors interpreted results. J.M.N.S. wrote the manuscript with contributions from all authors. All authors read and approved the final manuscript.

Additional information

The lattice-gas cellular automata code used in this study is available from the corresponding author on reasonable request.
Competing financial interests The authors declare no competing financial interests.

References

  • 1. Metzler, R. & Klafter, J. The random walk’s guide to anomalous diffusion: A fractional dynamics approach. Phys. Rep. 339, 1–77 (2000).
  • 2. Coffey, W., Kalmykov, Y. & Waldron, J. The Langevin Equation (World Scientific, Singapore, 2004).
  • 3. Metzler, R. & Klafter, J. The restaurant at the end of the random walk: Recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A: Math. Gen. 37, R161–R208 (2004).
  • 4. Klages, R., Radons, G. & Sokolov, I. (eds.) Anomalous transport (Wiley-VCH, Berlin, 2008).
  • 5. Klafter, J. & Sokolov, I. First Steps in Random Walks: From Tools to Applications (Oxford University Press, Oxford, 2011).
  • 6. Höfling, F. & Franosch, T. Anomalous transport in the crowded world of biological cells. Rep. Prog. Phys. 76, 046602/1–50 (2013).
  • 7. Metzler, R., Jeon, J.-H., Cherstvy, A. G. & Barkai, E. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. Phys. Chem. Chem. Phys. 16, 24128–24164 (2014).
  • 8. Zaburdaev, V., Denisov, S. & Klafter, J. Lévy walks. Rev. Mod. Phys. 87, 483–529 (2015).
  • 9. Shlesinger, M., Zaslavsky, G. & Klafter, J. Strange kinetics. Nature 363, 31–37 (1993).
  • 10. Klafter, J., Shlesinger, M. F. & Zumofen, G. Beyond Brownian motion. Phys. Today 49, 33–39 (1996).
  • 11. Sokolov, I., Klafter, J. & Blumen, A. Fractional kinetics. Phys. Today 55, 48–54 (2002).
  • 12. Scher, H. & Montroll, E. Anomalous transit-time dispersion in amorphous solids. Phys. Rev. B 12, 2455–2477 (1975).
  • 13. Kukla, V. et al. NMR studies of single-file diffusion in unidimensional channel zeolites. Science 272, 702–704 (1996).
  • 14. Balescu, R. Aspects of anomalous transport in plasmas, vol. 18 of Series in Plasma Physics (CRC Press, London, 2005).
  • 15. Barthelemy, P., Bertolotti, J. & Wiersma, D. A Lévy flight for light. Nature 453, 495–498 (2008).
  • 16. Viswanathan, G., da Luz, M., Raposo, E. & Stanley, H. The Physics of Foraging (Cambridge University Press, Cambridge, 2011).
  • 17. Brockmann, D., Hufnagel, L. & Geisel, T. The scaling laws of human travel. Nature 439, 462–465 (2006).
  • 18. Upadhyaya, A., Rieu, J., Glazier, J. & Sawada, Y. Anomalous diffusion and non-Gaussian velocity distributions of hydra cells in cellular aggregates. Physica A 293, 549–558 (2001).
  • 19. Thurner, S., Wick, N., Hanel, R., Sedivy, R. & Huber, L. Anomalous diffusion on dynamical networks: a model for interacting epithelial cell migration. Physica A: Statistical Mechanics and its Applications 320, 475 – 484 (2003).
  • 20. Dieterich, P., Klages, R., Preuss, R. & Schwab, A. Anomalous dynamics of cells migration P. Natl. Acad. Sci. USA 105, 459– 463 (2008).
  • 21. Takagi, H., Sato, M., Yanagida, T. & Ueda, M. Functional analysis of spontaneous cell movement under different physiological conditions. PLoS ONE 3, e2648/1–8 (2008).
  • 22. Bödeker, H., Beta, C., Frank, T. & Bodenschatz, E. Quantitative analysis of random ameboid motion. Europhys. Lett. 90, 28005/1–5 (2010).
  • 23. Harris, T. et al. Generalized Lévy walks and the role of chemokines in migration of effector cd8+ T cells. Nature 486, 545–548 (2012).
  • 24. Metzner, C. et al. Superstatistical analysis and modelling of heterogeneous random walks. Nat. Comm. 6, 7516–7523 (2015).
  • 25. Barbaric, I. et al. Time-lapse analysis of human embryonic stem cells reveals multiple bottlenecks restricting colony formation and their relief upon culture adaptation. Stem Cell Reports 3, 142–155 (2014).
  • 26. Rieu, J., Upadhyaya, A., Glazier, J., Ouchi, N. & Sawada, Y. Diffusion and deformations of single hydra cells in cellular aggregates. Biophys. J. 79, 1903–1914 (2000).
  • 27. Krummel, M., Bartumeus, F. & Gérard, A. T cell migration, search strategies and mechanisms. Nat. Rev. Immun. 16, 193––201 (2016).
  • 28. Méhes, E. & Vicsek, T. Collective motion of cells: from experiments to models. Integr. Biol. 6, 831–854 (2014).
  • 29. McCann, C., Kriebel, P., Parent, C. & Losert, W. Cell speed, persistence and information transmission during signal relay and collective migration. J. Cell Sci. 123, 1724–1731 (2010).
  • 30. Malinverno, C. et al. Endocytic reawakening of motility in jammed epithelia. Nat. Mat. 16, 587–597 (2017).
  • 31. Bi, D., Yang, X., Marchetti, M. & Manning, M. Motility-driven glass and jamming transitions in biological tissues. Phys. Rev. X 6, 021011/1–13 (2016).
  • 32. Fodor, É. et al. How far from equilibrium is active matter? Phys. Rev. Lett. 117, 038103/1–6 (2016).
  • 33. Fedotov, S. & Korabel, N. Emergence of Lévy walks in systems of interacting individuals. Phys. Rev. E 95, 030107/1–6 (2017).
  • 34. Ariel, G. et al. Swarming bacteria migrate by Lévy walk. Nat. Commun. 6, 8396/1–6 (2015).
  • 35. Mente, C., Voß-Böhme, A., & Deutsch, A. Analysis of individual cell trajectories in lattice-gas cellular automaton models for migrating cell populations. B. Math. Biol. 77, 660–697 (2015).
  • 36. Chen, J. C., & Kim, A. S. Brownian dynamics, molecular dynamics, and Monte Carlo modeling of colloidal systems. Adv. Colloid Interfac. 112, 159–173 (2004).
  • 37. Weitz, S., Deutsch, A., & Peruani, F. Self-propelled rods exhibit a phase-separated state characterized by the presence of active stresses and the ejection of polar clusters. Phys. Rev. E 92, 012322 (2015).
  • 38. Deutsch, A., & Dormann, S. Cellular Automaton Modeling of Biological Pattern Formation:
    Characterization, Applications, and Analysis
    (Birkhauser, Boston, 2017).
  • 39. Hatzikirou, H., Brusch, L., & Deutsch, A. From cellular automaton rules to a macroscopic mean-field description. Acta Phys. Pol. B Proc. Suppl. 3, 399–416 (2010).
  • 40. Reher, D., Klink, B., Deutsch, A., & Voß-Böhme A. Cell adhesion heterogeneity reinforces tumour cell dissemination: novel insights from a mathemtical model. Biol. Direct 12, 18 (2017).
  • 41. Böttger, K., Hatzikirou, H., Voß-Böhme, A., Cavalcanti-Adam, E. A., Herrero, M. A., & Deutsch, A. An emerging Allee effect is critical for tumor initiation and persistence. PLoS Comput. Biol. 11, e1004366 (2015).
  • 42. Böttger, K., Hatzikirou, H., Chauvière, A., & Deutsch, A. Investigation of the migration/proliferation dichotomy and its impact on avascular glioma invasion. Math. Model. Nat. Phenom. 7, 105–135 (2012).
  • 43. Hatzikirou, H., Basanta, D., Simon, M., Schaller, K., & Deutsch, A. ‘Go or Grow’: the key to the emergence of invasion in tumour progression? Math. Med. Biol. 29, 49–65 (2012).
  • 44. Pearson, K. Mathematical contributions to the theory of evolution - a mathematical theory of random migration. Biometric ser. 3, 54 (1906).
  • 45. Barber, M. N., & Ninham, B. W. Random and restricted walks - Theory and applications. (Gordon and Breach, New York, 1970).
  • 46. Lawniczak, A.T. Lattice gas automata for diffusive-convective transport dynamics. Center for Nonlinear Studies, Newsletter No. 136 (1997).
  • 47. Klages, R. & Korabel, N. Understanding deterministic diffusion by correlated random walks. J. Phys. A: Math. Gen. 35, 4823–4836 (2002).
  • 48. Reif, F. Fundamentals of statistical and thermal physics (McGraw-Hill, Auckland, 1965).
  • 49. Peruani, F., Deutsch, A., & Bär, M. A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms. Eur. Phys. J.-Spec. Top. 157, 111–122 ( 2008).
  • 50. Nava-Sedeño, J. M., Hatzikirou, H., Peruani, F., & Deutsch, A. Extracting cellular automaton rules from physical Langevin equation models for single and collective cell migration. J. Math. Biol., 1–26 (2017).
  • 51. Schüring, A. Molekulardynamik-Simulationen und Sprungmodelle zur Diffusion in Zeolithen. (Technische Universität Leipzig, 2003).
  • 52. Shalchi, A. Applicability of the Taylor-Green-Kubo formula in particle diffusion theory. Phys. Rev. E 83, 046402 (2011).
  • 53. Chechkin, A., Lenz, F. & Klages, R. Normal and anomalous fluctuation relations for Gaussian stochastic dynamics. J. Stat. Mech.: Theor. Exp. 2012, L11001/1–13 (2012).
  • 54. Korabel, N., Klages, R., Chechkin, A., Sokolov, I. & Gonchar, V. Fractal properties of anomalous diffusion in intermittent maps. Phys. Rev. E 75, 036213–1–14 (2007).
  • 55. Pressé, S., Ghosh, K., Lee, J., & Dill, K. A. Principles of maximum entropy and maximum caliber in statistical physics. Rev. Mod. Phys. 85, 1115–1141 (2013).
  • 56. Hazoglou, M. J., Walther, V., Dixit, P. D., & Dill, K. A. Communication: Maximum caliber is a general variational principle for nonequilibrium statistical mechanics. J. Chem. Phys. 143, 051104 (2015).
  • 57. Othmer, H. G., Dunbar, S. R. & Alt, W. J. Models of dispersal in biological systems. J. Math. Biol. 26, 263–298 (1988).
  • 58. Talkenberger, K., Cavalcanti-Adam, E. A., Deutsch, A., & Voß-Böhme, A. Amoeboid-mesenchymal migration plasticity promotes invasion only in complex heterogeneous microenvironments. Sci. Rep. (2017).
  • 59. Cahalan, M. D., Parker, I., Wei, S. H., & Miller, M. J. Real-time imaging of lymphocytes in vivo. Curr. Opin. Immunol. 15, 372–377 (2003).
  • 60. del Álamo, J. C., Meili, R., Alonso-Latorre, B., Rodríguez-Rodríguez, J., Aliseda, A., Firtel, R. A., & Lasheras, J. C. Spatio-temporal analysis of eukaryotic cell motility by improved force cytometry. P. Natl. Acad. Sci. USA 104, 13343–13348 (2007).
  • 61. Zienkiewicz, A., Barton, D. A. W., Porfiri, M., & di Bernardo, M. Data-driven stochastic modelling of zebrafish locomotion. J. Math. Biol. 71, 1081–1105 (2015).
  • 62. Christopher, R., Dhiman, A., Fox, J., Gendelman, R., Haberitcher, T., Kagle, D., Spizz, G., Khalil, I. G., & Hill, C. Data-driven computer simulation of human cancer cell. Ann. NY Acad. Sci. 1020, 132–153 (2004).
  • 63. Sisan, D. R., Halter, M., Hubbard, J. B., & Plant, A. L. Predicting rates of cell state change caused by stochastic fluctuations using a data-driven landscape model. P. Natl. Acad. Sci. USA 109, 19262–19267 (2012).

Supplementary information

Appendix A MSD in correlated systems

Here we derive the general expression for the MSD when there are temporal correlations. We use the derivation by Schüring [51]. To start with, we consider the “MSD” as defined by Shalchi [52]

which in LGCA notation are written as

(S1)

where and are the two orthonormal unit vectors in Cartesian coordinates. In 2D the diagonal elements are given by:

(S2)

and

(S3)

where . Adding them up gives the MSD

(S4)

which is just a Taylor-Green-Kubo formula [47, 54]. When we have , so by taking these terms out of the sum we get

(S5)

On the one hand and on the other , so using these relations on the first term on the right hand side yields

(S6)

Because the cosine is an even function , which means that we are adding two identical terms for every (this condition is already satisfied due to the Kronecker delta). This situation allows us to rewrite the limits of the interior sum if we take care of counting each term twice. Furthermore, using trigonometric identities it is possible to expand the cosine on the right hand side,

(S7)

If the orientations of the particle are completely uncorrelated from one another, then it is possible to write the sums on the right as

(S8)

If the reorientation probabilities are even functions of the angle , then we have that , and so the second sum on the right hand side disappears leaving

(S9)

We can choose our coordinate system such that without loss of generality. With this choice of the coordinate system, we can rewrite the MSD as

(S10)

Using the definition of the VACF this becomes

(S11)

Expanding the difference on the right hand side we get

(S12)

The last sum on the right hand side can be simplified, so we get

(S13)

Finally, using the relation between the intantaneous particle velocity and the diffusion constant in the random walk limit , and reordering terms, we obtain

(S14)

Appendix B VACF and MSD derivation in the persistent random walk

First we analytically derive the expected form of the VACF for a single particle in an LGCA where the lattice is a 2D square lattice.

Vacf

As mentioned before, the orientation probability is given by Eq. (10). The VACF is formally defined as , where and are the velocities of the particle at time 0 and , respectively. Using this definition, we can calculate the VACF of a stochastically moving particle as

where is the probability of the particle having a velocity at time .

In an LGCA particle velocities are given by the velocity channels they are located in, which belong to a finite set of unit vectors depending on the lattice dimension and geometry. Furthermore, time is also discrete with time steps of length such that at time step time has elapsed by . We can then rewrite the definition of the velocity autocorrelation in the case of an LGCA in the following way [47]:

(S15)

where is the orientation of the particle at time step and is the orientation of the particle at time step . To calculate the VACF, we start by defining a function as:

(S16)

where is the particle orientation at time step . Having defined this function, we can rewrite Eq. (10) as follows:

(S17)

where the partition function is defined as

(S18)

The expected value of the function is given by

(S19)

that is, the energy of the system is the single-step correlation. Due to the distribution of the reorientation probabilites the total energy can be calculated by the well-known relation

(S20)

Using the last two equations, we get an expression for the single step correlation:

(S21)

In this single-particle model, the partition function can be easily calculated. For a 2D square lattice the partition function reads

(S22)

Substituting Supplementary Eq. (S22) into Supplementary Eq. (S21), we have

(S23)

The particle orientations are normalized vectors. This allows us to rewrite the single step correlation as

(S24)

and the VACF as

(S25)

where . Supplemenary Eq. (S25) can be rewritten by adding zeroes in the following way:

which after rearranging terms, has the form

(S26)

Using trigonometric identities and the linearity of the expected value operator we can expand this expression to

(S27)

where is a sum of expected values of products of sines and cosines. Because the model is Markovian, the -th orientation is only correlated with the -th orientation. This allows us to write

(S28)

Now, because the reorientation probabilities Eq. (10) are even functions with respect to , the expected values become

which in turn implies

Using these relations we find that the VACF is given by

(S29)

Using Supplementary Eqs. (S23) and (S24) in Supplementary Eq. (S29) yields

(S30)

which can be written as

(S31)

if we define the exponent as

(S32)

The exponent depends on the lattice dimension and geometry, as follows:

  • In 1D the exponent is given by:

    (S33)
  • In 2D with a triangular lattice the exponent is:

    (S34)
  • In 2D with a square lattice the exponent is given by:

    (S35)
  • In 2D with an hexagonal lattice it takes the form:

    (S36)
  • With a cubic 3D lattice the exponent reads:

    (S37)

Mean square displacement

To calculate the MSD of particles performing persistent random walks, we start with Supplementary Eqs. (S6) and (S24) and rewrite the sum limits taking into account that the cosine is an even function to obtain