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

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.


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 r i r j = t t 0 t t 0 v i (τ)v j (ξ ) dξ dτ which in LGCA notation are written as whereê x andê y are the two orthonormal unit vectors in Cartesian coordinates. In 2D the diagonal elements are given by: and where θ k = arg c i k . Adding them up gives the MSD which is just a Taylor-Green-Kubo formula 47,54 . When i = j we have cos (θ i − θ j ) = 1, so by taking these terms out of the sum we get On the one hand v 2 τ 2 = ε 2 and on the other ε 2 = 2dDτ, so using these relations on the first term on the right hand side yields Because the cosine is an even function cos (θ i − θ j ) = cos (θ j − θ i ), which means that we are adding two identical terms for every i = j (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,

2/17
If the orientations of the particle are completely uncorrelated from one another, then it is possible to write the sums on the right as If the reorientation probabilities are even functions of the angle θ k , then we have that sin (θ i ) = 0, and so the second sum on the right hand side disappears leaving We can choose our coordinate system such that θ 0 = 0 without loss of generality. With this choice of the coordinate system, we can rewrite the MSD as Using the definition of the VACF this becomes Expanding the difference on the right hand side we get The last sum on the right hand side can be simplified, so we get Finally, using the relation between the intantaneous particle velocity v and the diffusion constant in the random walk limit D rw , and reordering terms, we obtain

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 v 0 and v t are the velocities of the particle at time 0 and t, respectively. Using this definition, we can calculate the VACF of a stochastically moving particle as where P( v,t) is the probability of the particle having a velocity v at time t.
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 3/17 that at time step k time has elapsed by kτ. We can then rewrite the definition of the velocity autocorrelation in the case of an LGCA in the following way 47 : where c i 0 is the orientation of the particle at time step k = 0 and c i k is the orientation of the particle at time step k. To calculate the VACF, we start by defining a function as: where c i k is the particle orientation at time step k. Having defined this function, we can rewrite Eq. (10) as follows: where the partition function is defined as The expected value of the function is given by 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 Using the last two equations, we get an expression for the single step correlation: In this single-particle model, the partition function can be easily calculated. For a 2D square lattice the partition function reads The particle orientations c i are normalized vectors. This allows us to rewrite the single step correlation as and the VACF as where θ k = arg [ c k ]. Supplemenary Eq. (S25) can be rewritten by adding zeroes in the following way:

4/17
which after rearranging terms, has the form Using trigonometric identities and the linearity of the expected value operator we can expand this expression to where f is a sum of expected values of products of sines and cosines. Because the model is Markovian, the i-th orientation is only correlated with the (i − 1)-th orientation. This allows us to write Now, because the reorientation probabilities Eq. (10) are even functions with respect to Using these relations we find that the VACF is given by Using Supplementary Eqs. (S23) and (S24) in Supplementary Eq. (S29) yields which can be written as if we define the exponent α as The exponent α depends on the lattice dimension and geometry, as follows: • In 1D the exponent is given by: • In 2D with a triangular lattice the exponent is: • In 2D with a square lattice the exponent is given by: • In 2D with an hexagonal lattice it takes the form: • With a cubic 3D lattice the exponent reads:

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 The expected value on the right hand side is the (m − n)-step correlation. From the previous VACF calculation we know that in the Markovian model which can be substituted in the expression of the MSD to obtain Because the sums on the right hand side only depend on the interval length | n − m | and not on the specific values of the indices n and m we can replace both sums by a sum over all posible interval lengths. There are k − j ways to divide an interval of k time steps (because the sums start from n = 1) into intervals of size j. Taking all into account, the MSD becomes which can also be written as where α is given by Supplementary Eq. (S32). We now distribute the two multiplying time steps τ on the second term on the right hand side, and multiply by one the exponent of the exponential function thus leaving it unchanged: We now use the definitions of the diffusion coefficient and the particle speed to obtain the following expression for the time step length: τ = 2dD v 2 , and use it to substitute for τ on the denominator of the exponent

VACF in homogeneous Markovian models
We will now consider a general Markovian model for a single moving particle. The model is then a Markov chain of particle orientations, i.e. the particle can transition between different orientations at each time step.
Definition 1. The state space of the Markov chain is E = { c 0 , c ±1 , · · · , c ±n , · · · , c N }, where the 2N different states are given by c n = cos π n , sin π n , n = 1, · · · , N − 1, Definition 2. The state space subset E 0 is defined as If the space is isotropic, then it is reasonable to require that the probability of the particle turning left or right be identical. Furthermore, we assume that the probability of turning does not depend on the specific time step, i.e. that the Markov process is homogeneous.
Definition 5. The velocity autocorrelation function (VACF) is given by Theorem 2. The velocity autocorrelation function of a particle whose orientations are given by a homogeneous, symmetric Markov chain is either delta-correlated, i.e. g k = δ 0,k , where δ is the Kronecker delta; alternating, i.e. g k = (−1) k a k , a ∈ R + ; or exponentially decaying, i.e. g k = e αk , α ≤ 0.
Proof. The proof is by induction.
Using the Chapman-Kolmogorov equation, we can calculate the VACF at the time step k + 1

7/17
We now expand the second sum on the right hand side of the equation Inserting this expression back into the VACF yields We can rewrite a as a = ∑ θ p(θ ) cos θ , where θ = arg( c 0 , c n ), ∀ c n ∈ E . Using the fact that 0 ≤ p(θ ) ≤ 1 and ∑ θ p(θ ) = 1, we have We have three cases: • −1 ≤ a < 0, then a = −1 | a | and g k = a k = (−1) k | a | k .

One dimension
We will now sketch our method for obtaining the reorientation probabilities P i k ,k in 1D. We start by expanding g(k) for the first two time steps after kτ = t ≥ ∆ (see Eq. (16)) for a 1D lattice. We will denote by the subscript f the lattice direction parallel to the original orientation of the particle. Similarly, the subscript r denotes the direction opposite to the original orientation of the particle. Numerical subscripts denote the time step at which the reorientation probability is evaluated.
Time step k = 1 Only two trajectories are possible after one time step. Their probabilities are given by P f ,1 and P r,1 . The normalization condition for these probabilties reads We now expand the VACF: We can substitute P r,1 from Supplementary Eq. (S42) into Supplementary Eq. (S43) Rearranging terms we obtain the probability for having the same orientation as originally to

8/17
Substituting Supplementary Eq. (S44) into Supplementary Eq. (S42) we obtain the probability for the particle to turn around after the first time step, After inspection of Supplementary Eqs. (S44) and (S45), and recalling that in the 1D lattice c 1 = 1 and c 2 = −1, both probabilities can be written as a single expression, where i is a placeholder variable for either f or r.
Time step k = 2 After two time steps, we have four different possible paths for the particle, with four different probabilities. If we assume the probabilities at each time can be written as P f f = P f ,1 P f ,2 , we can expand the VACF to obtain: P f ,1 P f ,2 − P f ,1 P r,2 + P r,1 P f ,2 − P r,1 P r,2 = P f ,2 − P r,2 P f ,1 + P r,1 = g (2), which by employing Supplementary Eq. (S42) can simplified to: Given that the probabilities in the previous time step were normalized, it is sufficient to require that the probabilities in the current time step be normalized: Inspecting Supplementary Eqs. (S47) and (S48) and comparing them with Supplementary Eqs. (S42) and (S43) we can see that they are identical except for the evaluation of g(k). Therefore, for the second time step it holds that Any k It is easy to see that for further times we can always assume that the probabilities are uncorrelated so that only the last orientation in the particle's orientation history is relevant for the calculation. If we do, Supplementary Eqs. (S46) and (S49) can be generalized for any time step k in the following way:

Two dimensions: Triangular lattice
We will repeat the calculation we did in 1D now in 2D for two different lattice geometries to identify possible dependencies on the lattice dimension and/or geometry.
Time step k = 1 We have three possible lattice directions with lattice vectors given by either on alternating nodes. In the first time step there are three possible paths given by P r,1 , P a,1 , and P u,1 , where P r,1 is the probability to reverse orientation. The normalization condition is, in this case, given by P r,1 + P u,1 + P a,1 = 1, while the VACF is given by 1 2 (P u,1 + P a,1 ) − P r,1 = g(1).
We need to make an assumption to continue, as there are more variables than equations. We assume that the probability of turning left or right is identical, P u,1 = P a,1 := P f ,1 . (S53)

(S55)
Substituting P f ,1 from Supplementary Eq. (S55) into Supplementary Eq. (S54) we obtain P r,1 + 2(P r,1 + g(1)) = 3P r,1 + 2g(1) = 1, which, after rearranging, gives the expression for the probability of the particle to go back: Using Supplementary Eq. (S56) in Supplementary Eq. (S54) we obtain the probability Examining Supplementary Eqs. (S53), (S56) and (S57) we can summarize them as Time step k = 2 In this case there are 9 different possible orientation histories with 9 different probabilities. If we now denote by f and r the lattice directions parallel and antiparallel to the original particle orientation, respectively, and by u and a the remaining lattice directions and assume that the probabilities are uncorrelated, we require that probabilites at the present time step are normalized: while the VACF has the form P u,1 P f ,2 − 1 2 P u,1 P u,2 − 1 2 P u,1 P a,2 + P a,1 P f ,2 − 1 2 P a,1 P u,2 − 1 2 P a,1 P a,2 + P r,1 P f ,2 − 1 2 P r,1 P u,2 − 1 2 P r,1 P a,2 = P f ,2 − 1 2 (P u,2 + P a,2 ) (P u,1 + P a,1 + P r,1 ) = g(2), which, by Supplementary Eq. (S51), is simplified to: To continue, we impose the isotropy condition Supplementary Eq. (S53) denoting by P r,2 the probabilities P u,2 and P a,2 . With these assumptions the normalization condition reads while the VACF is now Inserting P f ,2 from Supplementary Eq. (S61) into Supplementary Eq. (S62) we obtain the probability P r,2 : and, using the normalization condition Supplementary Eq. (S61) we obtain the probability P f ,2 : These probabilities can be written in the general form

10/17
Any k From Supplementary Eqs. (S58) and (S65) we can see that for any further time step k and making the same assumptions as before the probabilities are given by Two dimensions: Square lattice Time step k = 1 There are four possible lattice directions with lattice vectors c 1 = (1, 0), c 2 = (0, 1), c 3 = (−1, 0), and c 4 = (0, −1). Therefore there are four possible probabilities, so the normalization condition reads P r,1 + P f ,1 + P u,1 + P a,1 = 1, where P u,1 and P a,1 are the probabilities of going in the two directions orthogonal to the original orientation of the particle. We now expand the VACF to obtain P f ,1 − P r,1 = g(1).
Right from the start we have more variables than equations, so we need to make one more assumption in order to continue with the derivation. To simplify we assume the following: With this assumption the normalization condition becomes Rearranging terms we obtain the probability Inserting Supplementary Eq. (S71) into Supplementary Eq. (S70) we obtain the remaining probability Supplementary Equations (S69), (S71) and (S72) can then be summarized as Time step k = 2 There are now 16 different possible histories for the traveling particle. As before we assume that the probabilities are uncorrelated which, together with Supplementary Eq. (S67), allows us to write the normalization condition as P r,2 + P f ,2 + P u,2 + P a,2 = 1, and the VACF now is P f ,1 P f ,2 + P u,1 P f ,2 + P r,1 P f ,2 + P a,1 P f ,2 − (P f ,1 P r,2 + P u,1 P r,2 + P r,1 P r,2 + P a 1 P r,2 ) = P f ,2 − P r,2 P f ,1 + P u,1 + P r,1 + P a,1 = g(2), which, by using Supplementary Eq. (S67), is simplified to P f ,2 − P r,2 = g(2).
We see that Supplementary Eqs. (S67) and (S74), and (S68) and (S75) are practically identical. Therefore, by making the same assumptions, we arrive at the following expression for the probabilities at k = 2:
The VACF is given by Similarly as in the case of the square lattice, we impose the following condition which allows deriving the reorientation probabilities: which enables us to simplify the normalization condition in the following way: Using Supplementary Eq. (S81) to substitute P r,1 into Supplementary Eq. (S79) we obtain which, after rearranging terms yields the probability Now, using Supplementary Eq. (S81) we can obtain the remaining probability Examining Supplementary Eqs. (S80), (S82) and (S83) we arrive at the general expression Any k As done before we can continue the process for further times and, making the same assumptions, we arrive at an equation as Supplementary Eq. (S84) for any time k: Any dimension, any lattice geometry, any time Now that probabilities were derived for several dimensions, geometries, and times, we can see from Supplementary Eqs. (S50), (S66), (S77) and (S85) that the general form of the probabilities is given by where d is the spatial dimension and b is the number of nearest neighbors.

MSD of the piecewise process
We will now calculate the MSD using probabilities such that the VACF is a power-law decaying piecewise function defined as where t is such that C 0 ∆ t φ = 1. It is straightforward to see that the probabilities which define such a VACF obey where i 0 is the index of the velocity channel the particle started in and ω is such that t = ωτ. It is easy to see that in the first ω time steps the MSD is defined by r 2 (k) = k 2 ε 2 , or, using the definition of the particle speed and taking the limit τ → 0: We will now calculate the MSD for time steps greater than ω. The calculation will be made for a 1D lattice, but the results are identical for any dimension and lattice geometry. To ease notation, we will omit any subindices refering to time steps k ≤ ω, as we know that, given Supplementary Eq. (S88), only those trajectories where the first ω orientations of the particle are identical to the original orientation of the particle have non-zero probabilities.

13/17
Any k We can proceed for further k and will arrive at the following expression for any k > ω: which by using the definition of the diffusion coefficient and particle speed can be converted to Combining Supplementary Eqs. (S89) and (S93) we obtain the MSD of a particle with a piecewise power-law decaying VACF:

Generalized time-correlated random walk: rule derivation
We maximize the caliber subject to observing a certain VACF, which translates into the Lagrange multiplier problem This yields the trajectory probabilites where Z = exp (1 − λ ) is called the dynamical partition function which, by optimizing the functional with respect to λ (i.e., ∂C ∂ λ = 0), is given by Z = ∑ Γ exp ∑ k i=1 β (i) c n 0 · c n i . Optimizing with respect to β (i) yields our original constraint g(k) = ∑ Γ P Γ c n 0 · c n i .
Solving for β (i) using Supplementary Eqs. (S97) and (S98) can be quite challenging, so we expand P Γ in a Taylor series around β (i) = 0, which reduces to g(k) ≈ ∑ Γ 1 Z 1 + ∑ k i=0 β (i) c n 0 · c n i c n 0 · c n k , where the dynamical partition function is 14/17 simplified as where b is the number of lattice directions and θ i is the angle between the original particle orientation and the particle orientation at time step i. ∑ Γ cos θ i = 0 because the lattice directions and hence the possible particle orientations are symmetrically and homogeneously distributed. Substituting Z, using the same notation as previously, employing trigonometric identities, and denoting the spatial dimension by d, we proceed with the calculation: cos θ 1 cos θ k + β (2) ∑ Γ cos θ 2 cos θ k + · · · + β (k − 1) ∑ Γ cos θ k−1 cos θ k + β (k) ∑ Γ cos 2 θ k which determines the Lagrange multiplier β (k) = dg(k).
So the generalized probabilities are finally which is the probability for the whole trajectory. Due to the exponential form of this probability, we can decompose the trajectory probability into reorientation probabilities: