title: Inference of hopping rates of anisotropic random walk on a 2D lattice via covariance-based estimators of diffusion parameters
authors: Masanori Mishima
text: 
  
  
  
  
  
  
  
   Centre for Mechanochemical Cell Biology, Division of Biomedical Sciences, Warwick Medical School, University of Warwick
   
  
  
  Traditionally, time-development of the mean square displacement has been employed to determine the diffusion coefficient from the trajectories of single particles. However, this approach is sensitive to the noise and the motion blur upon image acquisition. Recently, Vestergaard et al.  has proposed a novel method based on the covariance between the shifted displacement series.  This approach gives a more robust estimator of the diffusion coefficient of one-dimensional diffusion without bias, i.e., when mean velocity is zero. Here, we extend this approach to a potentially biased random walk on a two-dimensional lattice. First, we describe the relationship between the  hopping rates to the eight adjacent sites and the time development of the higher-order moments of the stochastic two-dimensional displacements. Then, we derive the covariance-based estimators for these higher-order moments. Numerical simulations confirmed that the procedure presented here allows inference of the stochastic hopping rates from two-dimensional trajectory data with location error and motion blur.
  
  
  
  
  
  
  
  
  
  
  
                                
  Inference of hopping rates of anisotropic random walk on a 2D lattice via covariance-based estimators of diffusion parameters
      Masanori Mishima
      December 30, 2023
  =============================================================================================================================
  
  
  
  
  
  
  
  § INTRODUCTION
  
  The forward and backward hopping rates of a random walk on a one-dimensional lattice are linked to the average velocity (v) and diffusion coefficient (D). The latter macroscopic parameters can be estimated from the trajectory of the random walk, i.e., locations on the lattice measured with a regular time interval. Fitting of the mean square displacements (MSD) vs time lag data with MSD=σ^2+2Dτ +(vτ)^2, where σ is the localization noise, has been widely used to estimate the macroscopic parameters. However, it is known that there are practical issues such as the optimization of the range of the time lags and the influence of motion blurs when this approach is applied to the single particle dynamics of a biological molecule <cit.>. Recently, a novel approach based on the covariance between the displacements between adjacent time points (covariance-based estimator, CVE) has been proposed to be a more simple, accurate and robust alternative to estimate the diffusion coefficient of a unbiased random walk <cit.>.
  
  Here, we extend the CVE-based approach to the random walk on a two-dimensional (2D) lattice, in which a particle takes stochastic hops to the eight surrounding sites with distinct rates. We start with the well-known relation between velocity and the diffusion coefficient and the 1D hopping rates and the CVE of the diffusion coefficient in this case. We then show that in a 2D random walk, x- and y-velocities and diffusion coefficients plus four (higher-order) co-moments of the observed two-dimensional displacement series are linked to the eight hopping rates. The procedure for calculating the CVEs for the macroscopic parameters, i.e., the co-moments of the 2D displacements, is provided. This allows us to infer the eight hopping rates from the trajectories of a 2D random walk even with temporal and spatial resolutions at which individual hopping events can't be captured.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § RESULTS
  
  
  
   §.§ CVE-based approach for 1D diffusion with drift
  
  Here we consider a potentially biased random walk X(t) on a 1D lattice with the grid size, a (Fig. <ref>a). Let k_+1 and k_-1 be the stochastic forward and backward hopping rates , respectively (we can assume k_+1>k_-1>0 without loss of generality). X(t)  is the sum of x_i, the displacements in a small fraction of time δ t = t/n, i.e.,
  X(t) = x_1 + x_2 + … +x_n.
  
  
  
  Since x_i are indistinguishable but independent of each other, the ensemble averages of x_i
  
      ⟨ x_1 ⟩ = ⟨ x_2 ⟩ = … = ⟨ x_n ⟩≡⟨ x ⟩
      ⟨ x_1^2 ⟩ = ⟨ x_2^2 ⟩ = … = ⟨ x_n^2 ⟩≡⟨ x^2 ⟩
  
  where  ⟨ x ⟩ and ⟨ x^2 ⟩ are the mean displacement and the mean square displacement of X(t) in δ t, respectively. In general, ⟨ x_i^2 ⟩≠⟨ x_i ⟩^2, while for i ≠ j, ⟨ x_ix_j ⟩ = ⟨ x_i ⟩⟨ x_j ⟩. The ensemble average of X(t) is 
  
      ⟨ X(t) ⟩   =   ⟨ x_1 + x_2 + … +x_n ⟩ = n⟨ x ⟩
         =    n[ (+a)× k_+1δ t + (-a)× k_-1δ t] = a(k_+1 - k_-1) · nδ t  
         =    a(k_+1 - k_-1)t
  
  
  
  
  With ⟨ X(t)^2 ⟩ =  ⟨ (x_1 + x_2 + … +x_n)^2 ⟩ = ⟨∑_i x_i^2 + ∑_j ∑_i, i ≠ j x_i x_j ⟩ = n ⟨ x^2 ⟩ + n(n-1) ⟨ x ⟩^2, the variance of X(t) is
  
      var (X(t))    =   ⟨ (X(t)-⟨ X(t) ⟩)^2 ⟩ =⟨ X(t)^2 ⟩ - ⟨ X(t) ⟩ ^2 
         =     n ⟨ x^2 ⟩ + n(n-1) ⟨ x ⟩^2 - (n⟨ x ⟩)^2 = n ⟨ x^2 ⟩ - n ⟨ x ⟩^2 
         =    n[ (+a)^2× k_+1δ t + (-a)^2× k_-1δ t] - n [ (+a)× k_+1δ t + (-a)× k_-1δ t]^2 
         =    a^2(k_+1 + k_-1) · nδ t  -   a^2(k_+1 - k_-1)^2 (nδ t)δ t 
         →    a^2(k_+1 + k_-1)t   (δ t → 0)
  
  Thus, the velocity of the constant drift and the diffusion coefficient, v and D, respectively, are linked to the hopping rates
  
      v  = a (k_+1 - k_-1) 
      
          2D = a^2 (k_+1 + k_-1)
  
  As is well known, this indicates that we can infer the stochastic hopping rates by determining the macroscopic parameters, v and D. Hereafter, for simplicity, we omit a by assuming that X(t) is a dimensionless value measured with a as the unit, i.e., the physical location on the lattice is X(t)· a. Thus,
  
      v  = k_+1 - k_-1
      
          2D = k_+1 + k_-1
  
  
  In typical single molecule/particle observations, we can't determine the true coordinate X(t) as a continuous function. We can only measure the positions at a discrete time points t=0, Δ t, 2Δ t, … , n Δ t and the measured values X_0, X_1, X_2, … X_n suffer from the motion blur due to the movement during image capturing and the error in determining the position of the molecule/particle by image analysis. Under this circumstance, X_k is related to the true X(t) as
  
      X_k = ∫_0^Δ t s(t) X(t + (k-1)Δ t) dt + ε_k
  
  where s(t) defines the state of the shutter (s(t)=0 means closed shutter and s(t)>0 means open shutter, for normalization ∫_0^Δ t s(t) dt =1) and ε_k the Gaussian error in localization by image analysis (⟨ε_i⟩ =0 and ⟨ε_i^2⟩ =σ^2 represents the precision of the measurement when the target is immobile) <cit.>. The traditional approaches to determine v and D, for example, by fitting a quadratic curve to the mean square displacement data, suffer from complications arising from non-zero s(t) and ε. However, for unbiased diffusion (i.e. v=0), it has recently been shown <cit.> that a combination of the adjacent displacements Δ X_k = X_k+1 - X_k and Δ X_k+1 = X_k+2 - X_k+1 cancels out the terms containing s(t) and ε and results in a simple relation
  
      ⟨Δ X_k^2 ⟩ + 2 ⟨Δ X_k+1·Δ X_k ⟩ = 2D Δ t.
  
  D calculated with this relation provides a more reliable estimator of the diffusion coefficients.
  
  The above relation doesn't hold in the presence of a bias in the hopping rates. However, if we consider Z(t) = X(t) - vt, the deviation from the constant drift, and its observed counterpart Z_k = X_k- v· kΔ t, we obtain ⟨ Z_k ⟩ = 0, and thus the series of its displacements Δ Z_k = Z_k+1 - Z_k = Δ X_k - vΔ t satisfy 
  
      ⟨Δ Z_k^2 ⟩ + 2 ⟨Δ Z_k+1·Δ Z_k ⟩ = 2D Δ t.
  
  Since ⟨ (X(t) - ⟨ X(t) ⟩)^2 ⟩ =  ⟨ (Z(t) + vt - ⟨ Z(t) +vt⟩)^2 ⟩ =  ⟨ (Z(t) - ⟨ Z(t) ⟩)^2 ⟩, the value of D determined with this formula gives the diffusion coefficient of the biased random walk X(t). 
  
  
  
   §.§  Biased random walk on a 2D lattice: hopping rates and co-moments of the displacements
  
  Here we discuss a biased random walk (X(t), Y(t)) on a 2D lattice measured with the x- and y-grid size a and b, respectively. Hopping can occur to the eight surrounding sites with the distinct rates k_1, k_2, k_3, k_4, k_5, k_6, k_7, k_8. Let the velocities and diffusion coefficients of the movement along the grid axes be v_x and  v_y, and D_x, and D_y, respectively. Then 
  
      ⟨ X(t) ⟩ = v_x t = (k_F-k_B) t = (k_1 + k_2 + k_3 - k_6 - k_7 - k_8) t 
      ⟨ ( X(t) - ⟨ X(t) ⟩)^2 ⟩ = 2 D_x t =  (k_F+k_B) t = (k_1 + k_2 + k_3 + k_6 + k_7 + k_8) t
      ⟨ Y(t) ⟩ = v_y t = (k_L-k_R) t= (k_1 + k_4 + k_6 - k_3 - k_5 - k_8) t
      ⟨ ( Y(t) - ⟨ Y(t) ⟩)^2 ⟩ = 2 D_y t = (k_L+k_R) t=   (k_1 + k_4 + k_6 + k_3 + k_5 + k_8) t
  
  as we discussed in the above section.
  
  Here, we consider the covariance between X(t) and Y(t).
  
      cov(X, Y) = ⟨ (X - ⟨ X ⟩)(Y - ⟨ Y ⟩) ⟩ = ⟨ XY ⟩ - ⟨ X⟩⟨ Y ⟩
         =   ⟨ (x_1+x_2+ … + x_n)(y_1+y_2+ … +y_n) ⟩ - ⟨ x_1 + x_2 + … + x_n ⟩⟨ y_1 + y_2 + … + y_n ⟩
         =   ⟨∑_i=1^n x_i y_i + ∑_i≠ jx_i y_j  ⟩ - n⟨ x ⟩· n ⟨ y ⟩
         =    n ⟨ xy ⟩ + n(n-1) ⟨ x⟩⟨ y ⟩ - n^2 ⟨ x ⟩⟨ y ⟩
         =    n ⟨ xy ⟩ -  n ⟨ x⟩⟨ y ⟩
         =    n ·[ (+1)(+1)· k_1δ t + (-1)(+1)· k_3δ t + (+1)(-1)· k_6δ t + (-1)(-1)· k_8δ t ] 
                 - n · (k_F - k_B)δ t · (k_L - k_R)δ t 
         =    (k_1-k_3 -k_6+k_8) nδ t - (k_F - k_B) (k_L - k_R) nδ t ·δ t 
         →    (k_1-k_3 -k_6+k_8)t  (δ t → 0)
  
  
  
  We move on to higher-order moments. 
  Since
  
      ⟨ X^2 Y ⟩   =   ⟨∑_k ∑_j ∑_i x_i x_j y_k ⟩
         =    n  ⟨ x^2 y ⟩ +n(n-1) ⟨ x^2⟩⟨ y ⟩ + 2n(n-1)⟨ x ⟩⟨ xy ⟩ + n(n-1)(n-2) ⟨ x ⟩ ^2 ⟨ y ⟩ ,
  
  we get
  
      ⟨ (X - ⟨ X ⟩)^2( Y - ⟨ Y ⟩) ⟩ = ⟨ X^2 Y ⟩ - ⟨ X^2 ⟩⟨ Y ⟩ -2 ⟨ X ⟩ (⟨ XY ⟩ - ⟨ X ⟩⟨ Y ⟩ )  
         =    n ⟨ x^2 y ⟩ - n ⟨ x^2 ⟩⟨ y ⟩ - 2n ⟨ x⟩⟨ xy ⟩ +2n ⟨ x ⟩^2 ⟨ y ⟩
          =    (k_1 - k_3 +k_6 - k_8) nδ t - [(k_F + k_B)(k_L - k_R) + 2(k_F - k_B)(k_1 - k_3 -k_6 + k_8)] n δ t ·δ t  
                      +2(k_F - k_B)^2 (k_L - k_R) nδ t ·δ t^2 
         →    (k_1 - k_3 +k_6 - k_8)t  (δ t → 0).
  
  
  
  Similarly, 
  
      ⟨ (X - ⟨ X ⟩)( Y - ⟨ Y ⟩)^2 ⟩→ (k_1+k_3 - k_6 -k_8) t   (δ t → 0)
  
  
  Finally, using
  
      ⟨ X^2 Y^2 ⟩ = ⟨∑_l ∑_k ∑_j ∑_i x_i x_j y_k y_l ⟩
         =    n⟨ x^2 y^2 ⟩ + 2n(n-1) ⟨ x ⟩⟨ x y^2 ⟩ + 2n(n-1) ⟨ y ⟩⟨ x^2 y ⟩ + 2n(n-1)⟨ xy⟩^2 
                + n(n-1)⟨ x^2 ⟩⟨ y^2 ⟩ + n(n-1)(n-2) ⟨ x^2 ⟩⟨ y ⟩^2 +4n(n-1)(n-2) ⟨ xy ⟩⟨ x⟩⟨ y ⟩
                      + n(n-1)(n-2)⟨ x ⟩^2 ⟨ y^2 ⟩ + n(n-1)(n-2)(n-3) ⟨ x ⟩^2 ⟨ y ⟩^2
  
  we get 
  
      ⟨ (X- ⟨ X ⟩)^2  (Y- ⟨ Y ⟩)^2 ⟩
         =   ⟨ X^2 Y^2 ⟩ - 2⟨ X^2 Y ⟩⟨ Y ⟩ -2⟨ X Y^2 ⟩⟨ X ⟩  + ⟨ X^2 ⟩⟨ Y ⟩^2  + 4⟨ XY ⟩⟨ X ⟩⟨ Y ⟩ +⟨ Y^2 ⟩⟨ X ⟩^2 -3 ⟨ X⟩^2  ⟨ Y⟩^2  
         =    n ⟨ x^2 y^2 ⟩ + 2n(n-1) ⟨ xy ⟩ ^2 + n(n-1)⟨ x^2⟩⟨ y^2 ⟩
                   - 2n [ ⟨ x^2 y ⟩⟨ y ⟩+⟨ xy^2 ⟩⟨ x ⟩] - n(n-2) [⟨ x^2⟩⟨ y ⟩^2 +4 ⟨ xy ⟩⟨ x ⟩⟨ y ⟩ + ⟨ y^2⟩⟨ x ⟩^2 -3⟨ x ⟩^2 ⟨ y ⟩^2] 
         →    (k_1+k_3+k_6+k_8) t + 2(At)^2 + 2D_x t · 2D_y t   (δ t → 0)
  
  
  In summary, we have obtained a set of formulae that relate macroscopic observations to the microscopic hopping rates:
  
      ⟨ X ⟩   =    v_x t 
      ⟨ Y ⟩   =    v_y t 
      ⟨ (X - ⟨ X ⟩)^2 ⟩   =    2D_x t 
      ⟨ (Y - ⟨ Y ⟩)^2 ⟩   =    2D_y t 
      ⟨ (X - ⟨ X ⟩)(Y - ⟨ Y ⟩) ⟩   =    A t 
      ⟨ (X - ⟨ X ⟩)^2(Y - ⟨ Y ⟩) ⟩   =    B t 
      ⟨ (X - ⟨ X ⟩)(Y - ⟨ Y ⟩)^2 ⟩   =    C t 
      ⟨ (X - ⟨ X ⟩)^2(Y - ⟨ Y ⟩)^2 ⟩   - 2 ⟨ (X - ⟨ X ⟩)(Y - ⟨ Y ⟩)  ⟩^2              
      
       - ⟨ (X - ⟨ X ⟩)^2 ⟩⟨ (Y - ⟨ Y ⟩)^2 ⟩   =    E t
  
  where
  
      [  v_x;  v_y; 2D_x; 2D_y;    A;    B;    C;    E;      ] 
      =
      [  1  1  1  0  0 -1 -1 -1;  1  0 -1  1 -1  1  0 -1;  1  1  1  0  0  1  1  1;  1  0  1  1  1  1  0  1;  1  0 -1  0  0 -1  0  1;  1  0 -1  0  0  1  0 -1;  1  0  1  0  0 -1  0 -1;  1  0  1  0  0  1  0  1;    ][ k_1; k_2; k_3; k_4; k_5; k_6; k_7; k_8;     ]
  
  
  
  
  
  
  
   §.§ Covariance-based estimator for A
  
  Based on formula <ref> and <ref>, the velocities of constant drift, v_x and v_y, can be estimated with
  
      ⟨Δ X ⟩ = v_x Δ t 
      ⟨Δ Y ⟩ = v_y Δ t
  
  As discussed in section A, drift-adjusted displacements Δ Z_k = Δ X_k - v_x Δ t and Δ W_k = Δ Y_k - v_y Δ t based on Z(t)=X(t) - v_x t and W(t)=Y(t) - v_y t give the diffusion coefficients along x- and y-axes, D_x and D_y, respectively.
  
      ⟨Δ Z_k^2 ⟩ + 2⟨Δ Z_k+1Δ Z_k ⟩ = 2D_x Δ t 
      ⟨Δ W_k^2 ⟩ + 2⟨Δ W_k+1Δ W_k ⟩ = 2D_y Δ t
  
  Here we derive analogous covariance-based estimators of A, B, C, and E. Considering the symmetry between Z(t) and W(t), we guess the formula for A as
  
      ∑_α, β=0 or 1 
      not α=β=1⟨Δ Z_k+αΔ W_k+β⟩ = ⟨Δ Z_k Δ W_k ⟩ + ⟨Δ Z_k+1Δ W_k⟩ + ⟨Δ Z_k Δ W_k+1⟩ = A Δ t
  
  We shall now prove this.
  
  With
  
      Δ Z_k+α   =   ∫_0^Δ t s(t) [
              Z(t+(k+α)Δ t) - Z(t+(k+α-1)Δ t)
          ] dt +(ε_k+α+1 - ε_k+α) 
      Δ W_k+β   =   ∫_0^Δ t s(t') [
              W(t'+(k+β)Δ t) - W(t'+(k+β-1)Δ t) ] dt' +(ε'_k+β+1 - ε'_k+β)
  
  
  
      ⟨Δ Z_k+αΔ W_k+β⟩
         =   ∫_0^Δ t∫_0^Δ t s(t)s(t') ⟨( Z(t+(k+α)Δ t) -   Z(t+(k+α-1)Δ t))  
              ×( W(t'+(k+β)Δ t) -   W(t'+(k+β-1)Δ t) )
         ⟩ dt dt' +⟨ (ε_k+α+1 - ε_k+α)(ε'_k+β+1 - ε'_k+β)  ⟩
  
  Since (Z(t), W(t)) is a random walk that fulfils ⟨ Z(t)W(t) ⟩ =At,  ⟨ Z(t)W(t') ⟩ =A·min(t, t'). Using this, the factor to be integrated in the first term of ⟨Δ Z_k+αΔ W_k+β⟩ is evaluated to be
  
      ⟨( Z(t+(k+α)Δ t) -   Z(t+(k+α-1)Δ t))   ×( W(t'+(k+β)Δ t) -   W(t'+(k+β-1)Δ t) ) ⟩
         =   ⟨ Z(t+(k+α)Δ t) W(t'+(k+β)Δ t) ⟩  - ⟨ Z(t+(k+α)Δ t)  W(t'+(k+β-1)Δ t) ⟩
             - ⟨  Z(t+(k+α-1)Δ t) W(t'+(k+β)Δ t) ⟩ + ⟨  Z(t+(k+α-1)Δ t) W(t'+(k+β-1)Δ t) ⟩
         =    A [{ kΔ t + min(t + αΔ t, t'+βΔ t)} - {kΔ t + min(t + αΔ t, t'+ (β-1) Δ t)}
                                 -  {kΔ t +min(t + (α-1) Δ t, t'+ βΔ t)}  + {(k-1)Δ t + min(t + αΔ t, t'+βΔ t) }]
         =    A [ min(t + αΔ t, t'+βΔ t) - min(t + αΔ t, t'+ (β-1) Δ t)  
                                          - min(t + (α-1) Δ t, t'+ βΔ t)  + min(t + αΔ t, t'+βΔ t)  -Δ t ]
  
  Here, considering 0 ≤ t  ≤Δ t, 0 ≤ t' ≤Δ t, the factor within the above brackets 
  
      g(α, β, t, t')    =   min(t + αΔ t, t'+βΔ t) - min(t + αΔ t, t'+ (β-1) Δ t)  
                      - min(t + (α-1) Δ t, t'+ βΔ t)  + min(t + αΔ t, t'+βΔ t)  -Δ t
  
  is calculated to be 
  
      g(0, 0, t, t^')    =   Δ t - (t + t^' - 2·min(t, t^') )
      
          g(1, 0, t, t^')    =    t^' - min(t, t^') 
      
          g(0, 1, t, t^')    =    t - min(t, t^')
  
  Thus, g(0, 0, t, t') + g(1, 0, t, t') + g(0, 1, t, t') = Δ t.
  
  The x- and y-localization errors at the same time point are not independent. With Kronecker's delta δ (α, β) = 1 (α=β) or 0 (α≠β) ,
  
      ⟨ ( ε_k+α+1 - ε_k+α)  ( ε'_k+β+1 - ε'_k+β) ⟩
         =   ⟨ε_k+α+1ε'_k+β+1⟩ - ⟨ε_k+α+1ε'_k+β⟩ - ⟨ε_k+αε'_k+β+1⟩ + ⟨ε_k+αε'_k+β⟩
         =   σ_xy^2 (δ(α, β) - δ(α+1, β) - δ(α, β+1) + δ(α, β))
         =   σ_xy^2 (2 δ(α, β) - δ(α+1, β) - δ(α, β+1)) = σ_xy^2 Ω(α, β),
  
  where we define Ω(α, β)=2 δ(α, β) - δ(α+1, β) - δ(α, β+1).
  Thus, Ω(0, 0) + Ω(0, 1)+ Ω(1, 0) = (2· 1 - 0 - 0) + ( 2· 0 - 1 - 0) + (2· 0 - 0 - 1) =0.
      
  Finally, we get
  
      ∑_α, β=0 or 1 
      not α=β=1⟨Δ Z_k+αΔ W_k+β⟩
         =   ∫_0^Δ t∫_0^Δ t s(t)s(t') A ·[g(0, 0, t, t') + g(1, 0, t, t') + g(0, 1, t, t')] dt dt' 
                                                                    +σ_xy^2 [Ω(0, 0) + Ω(1, 0)+ Ω(0, 1)] 
         =    A Δ t
  
  
  
  
   §.§ Covariance-based estimators for B and  C 
  
  The above results imply that for random walks U(s), V(t),…, if a function f(x, y,… ) exists such that f(U(s), V(t), …) = K ·min(s, t, …) for some constant K (for example, for f(x, y)=⟨ xy⟩,  f(X(s), X(t))=2D_x  ·min(s,t),  and f(X(s), Y(t))=A  ·min(s, t)), then the equation below holds
  
      ∑_α, β, … =0 or 1 
      not α=β=… = 1 f(Δ U_k+α, Δ V_k+β,…) = KΔ t
  
  where Δ U_k = U_k+1 - U_k, Δ V_k = V_k+1 - V_k… are the observed displacements (adjusted for the constant drift), providing the general recipe for the CVE for K. Below is a proof for the cases of B and C.
  
  Since ⟨ Z(t)^2 W(t) ⟩ = B t, we choose f(x, y, z) = ⟨ xyz ⟩.
  
      f(Δ Z_k+α, Δ Z_k+β, Δ W_k+γ ) = ⟨Δ Z_k+α·Δ Z_k+β·Δ W_k+γ⟩
         =   ⟨(∫_0^Δ t s(t)[Z(t + (k + α)Δ t) - Z(t + (k + α-1)Δ t)] dt +(ε_k+α+1- ε_k+α) ) 
                   ×( ∫_0^Δ t s(t')[Z(t' + (k + β)Δ t) - Z(t' + (k + β-1)Δ t)] dt' +(ε_k+β+1- ε_k+β) ) 
                             ×( ∫_0^Δ t s(t”)[W(t” + (k + γ)Δ t) - W(t” + (k + γ-1)Δ t)] dt” +(δ_k+γ+1- δ_k+γ) ) ⟩
         =   ∫_0^Δ t∫_0^Δ t∫_0^Δ t s(t)s(t')s(t”) 
                  ⟨[Z(t + (k + α)Δ t) - Z(t + (k + α-1)Δ t)] ×[Z(t' + (k + β)Δ t) - Z(t' + (k + β-1)Δ t)] 
                       ×[W(t” + (k + γ)Δ t) - W(t” + (k + γ-1)Δ t)] 
        ⟩ dt dt' dt”
                                 + ⟨ (ε_k+α+1- ε_k+α) (ε_k+β+1- ε_k+β) (δ_k+γ+1- δ_k+γ) ⟩,
   
  in which we used ⟨ Z(t+(k+α)Δ t) - Z(t+(k+α-1)Δ t) ⟩ = 0, etc. and ⟨ε_k+α+1- ε_k+α⟩ = 0, etc.
  
   For unbiased random walks A(s), B(t), C(u) (⟨ A(s) ⟩=⟨ B(t) ⟩=⟨ C(u) ⟩ = 0, and their displacements are independent for non-overlapping time sections), if we assume s<t <u,
  
      ⟨ A(s) B(t) C(u) ⟩ = ⟨ A(s) ⟩⟨ (B(t) - B(s))(C(t)-C(s)) ⟩ 
          + ⟨ A(s)C(s)⟩ ( ⟨ B(t) ⟩ - ⟨ B(s) ⟩) 
                + ⟨ A(s)B(s)⟩ ( ⟨ C(t) ⟩ - ⟨ C(s) ⟩)
          + ⟨ A(s) B(s) C(s) ⟩
          = ⟨ A(s) B(s) C(s) ⟩
  
  In general, with τ = min(s, t, u), ⟨ A(s) B(t) C(u) ⟩ = ⟨ A(τ) B(τ) C(τ) ⟩.
  Applying this to Z(t) and W(t) that satisfy ⟨ Z(t)^2 W(t) ⟩ = B t, we get ⟨ Z(t)Z(t')W(t”) ⟩ = B·min (t, t', t”).
  Using this, we evaluate
  
      ⟨[Z(t + (k + α)Δ t) - Z(t + (k + α-1)Δ t)] ×[Z(t' + (k + β)Δ t) - Z(t' + (k + β-1)Δ t)] 
                       ×[W(t” + (k + γ)Δ t) - W(t” + (k + γ-1)Δ t)]  ⟩
          =   ⟨  Z(t + (k + α)Δ t) · Z(t' + (k + β)Δ t) · W(t” + (k + γ)Δ t)   ⟩
               -⟨  Z(t + (k + α)Δ t) · Z(t' + (k + β)Δ t) · W(t” + (k + γ-1)Δ t)   ⟩
               -⟨  Z(t + (k + α)Δ t) · Z(t' + (k + β-1)Δ t) · W(t” + (k + γ)Δ t)   ⟩
               +⟨  Z(t + (k + α)Δ t) · Z(t' + (k + β-1)Δ t) · W(t” + (k + γ-1)Δ t)   ⟩
                -  ⟨  Z(t + (k + α-1)Δ t) · Z(t' + (k + β)Δ t) · W(t” + (k + γ)Δ t)   ⟩
               +⟨  Z(t + (k + α-1)Δ t) · Z(t' + (k + β)Δ t) · W(t” + (k + γ-1)Δ t)   ⟩
               +⟨  Z(t + (k + α-1)Δ t) · Z(t' + (k + β-1)Δ t) · W(t” + (k + γ)Δ t)   ⟩
               - ⟨  Z(t + (k + α-1)Δ t) · Z(t' + (k + β-1)Δ t) · W(t” + (k + γ-1)Δ t)   ⟩
         =     B · g(α, β, γ, t, t', t”),
            where
            g(α, β, γ, t, t', t”)  
         =   min(t+αΔ t,t'+βΔ t, t”+ γΔ t) - min(t+αΔ t,t'+βΔ t, t”+(γ-1)Δ t)
         -   min(t+αΔ t,t'+(β-1)Δ t, t”+γΔ t) + min(t+αΔ t,t'+(β-1)Δ t, t”+(γ -1)Δ t)
         -   min(t+(α-1)Δ t,t'+βΔ t, t”+γΔ t) + min(t+(α-1)Δ t,t'+βΔ t, t”+(γ-1)Δ t)
         +   min(t+(α-1)Δ t,t'+(β-1)Δ t, t”+γΔ t) - min(t+(α-1)Δ t,t'+(β-1)Δ t, t”+(γ-1)Δ t).
  
  g(α, β, γ, t, t', t”) can be calculated to be
  
      g(0, 0, 0, t, t', t”)    =   Δ t - (t + t' +t”) + min(t', t” )  + min(t, t” )  + min(t, t' ) 
      
      g(1, 0, 0, t, t', t”)    =   min(t', t”) - min(t, t', t”)
      
      g(0, 1, 0, t, t', t”)    =   min(t, t”) - min(t, t', t”) 
      
      g(0, 0, 1, t, t', t”)    =   min(t, t') - min(t, t', t”) 
      
      g(1, 1, 0, t, t', t”)    =    t” - min(t, t” )  - min(t', t”)  + min(t, t', t” )  
      
      g(1, 0, 1, t, t', t”)    =    t' - min(t, t' )  - min(t', t”)  + min(t, t', t” )  
      
      g(0, 1, 1, t, t', t”)    =    t - min(t, t' )  - min(t, t”)  + min(t, t', t” ) .
  
  Thus,
  
      ∑_α, β, γ=0 or 1 
      not α=β=γ= 1 g(α, β, γ, t, t', t”) = Δ t
  
  
  The error term 
  
      Ω(α, β, γ) =  ⟨ ( ε_k+α+1 - ε_k+α) ( ε_k+β+1 - ε_k+β) ( δ_k+γ+1 - δ_k+γ) ⟩
  
  is calculated to be
  
      Ω(0, 0, 0)    =    0 
      Ω(1, 0, 0)    =   Ω(0, 1, 0) = Ω(0, 0, 1) = -⟨ε_k^2 δ_k ⟩
      Ω(0, 1, 1)    =   Ω(1, 0, 1) = Ω(1, 1, 0) = ⟨ε_k^2 δ_k ⟩.
  
  Thus,
  
      ∑_α, β, γ=0 or 1
      not α=β=γ= 1Ω(α, β, γ) = 0
  
  Finally, we obtain
  
      ∑_α, β, γ=0 or 1
      not α=β=γ= 1⟨Δ Z_k+α·Δ Z_k+β·Δ W_k+γ⟩
                                =  ∫_0^Δ t∫_0^Δ t∫_0^Δ t s(t)s(t')s(t”) B Δ t dt dt' dt ” = BΔ t.
  
  Similarly, 
  
      ∑_α, β, γ=0 or 1
      not α=β=γ= 1⟨Δ Z_k+α·Δ W_k+β·Δ W_k+γ⟩ = CΔ t
  
  
  
  
   §.§ Covariance-based estimators for E
  
  Z and W are linked to E by
  
      ⟨ Z(t)^2 W(t)^2 ⟩ - 2⟨ Z(t)W(t)⟩ ^2 - ⟨ Z(t)^2 ⟩⟨ W(t)^2⟩ = Et
  
  With f(x, y, z, w) = ⟨ xyzw ⟩ - ⟨ xy⟩⟨ zw⟩ 
      - ⟨ xz⟩⟨ yw⟩  - ⟨ xw⟩⟨ yz⟩, this relation is written
  
      f(Z(t), Z(t), W(t), W(t)) = Et
  
  Thus, the relation that gives the covariance-based estimator of E is predicted to be
  
      ∑_α, β, γ, λ =0 or 1 
      not α=β=γ=λ= 1 f( Δ Z_k+α,  Δ Z_k+β  , Δ W_k+γ, Δ W_k+λ ) = EΔ t
  
  
  This can be proven by some algebra similar to the above proofs for the CVEs of A, B and C  (see appendix A for details), confirming that the recipe for constructing CVEs (<ref>) holds for the necessary cases to infer the anisotropic hopping rates based on the relation (<ref>) (Fig. <ref>).
  
  
  
   §.§ Inference of the hopping rates from the 2D trajectories simulated with motion blur and localization error
  
  We now examine the versatility of the covariance-based estimators of D_x, D_y, A, B, C,and  E (formulas <ref>, <ref>,  <ref>, <ref>, <ref>, and <ref>), the coefficients of the time development of the co-moments of the displacements, and whether we can infer the anisotropic hopping rates from the trajectory data of 2D random walks, using the relations that link the macroscopic coefficients to the microscopic hopping rates (formula <ref>). Let us first compare two examples of random walks, RW1 and RW2, respectively, that are generated by Monte Carlo simulations to have the same x- and y-velocities and diffusion coefficients (v_x, v_y, D_x, D_y), but have different hopping rates (Fig. <ref> a and b). The values of the velocities, v_x = 16 s^-1 and v_y = 1 s^-1 (hops per second) correspond, for example, to the helical motion of a kinesin-like motor protein with a helical pitch of ∼1.7 μm around a microtubule that has a 2D lattice of the discrete motor-binding sites on its surface consisting of parallel 13 protofilaments, i.e., linear arrays of tubulin subunits aligned at 8 nm periodicity <cit.>. To mimic a typical condition of image acquisition using an EM-CCD camera used for the single molecule/particle observation, a series of ‘true’ positions are generated by the Gillespie algorithm <cit.> (blue lines in Fig. <ref> c and d) and the ‘observed’ positions at regular time points (every 0.1 s for 200 time points, corresponding to observation for 20 sec) were calculated by averaging the positions during the open shutter (90% of the cycle) and by adding Gaussian noises of the standard deviation of half the size of the 
  lattice spacing (red lines in Fig. <ref> c and d).
  
  
  
  
  
  
  
  Fifty trajectories of RW1 and RW2 were generated, respectively, and the macroscopic parameters v_x, v_y, D_x, D_y, A, B, C,and  E for each trajectory were calculated based on the formulas <ref>, <ref>, <ref>, <ref>, <ref>, <ref>, <ref>, and <ref>. As shown in Fig. <ref> e, the covariance-based estimators calculated from the simulated trajectories (gray dots) distributed around the theoretical values predicted by the formula (<ref>) (solid lines). Although, due to the intrinsic stochasticity of the random walk, the calculated values from individual trajectories spread around the theoretical `true' values, the means of 50 trajectories were found close to the corresponding theoretical values within twice the standard error. Importantly, as expected, while there was no significant difference between RW1 and RW2 in v_x, v_y, D_x, and  D_y, a clear difference between them was found in the distribution of A. This suggests that our recipe gives reasonable estimates of the coefficients of the temporal development of the co-moments of the x- and y-displacements in a realistic scenario. This allows us to detect a difference between the two random walks that have distinct preferences in the hopping direction but look exactly the same if we only consider the velocities and the diffusion coefficients along the two axes, separately.
  
  
  
  
  We move on to examine the power of our approach to reveal the hopping behaviors of anisotropic random walks. With known macroscopic parameters (v_x, v_y, D_x, D_y, A, B, C, E), we can calculate (k_1, k_2, k_3, k_4, k_5, k_6, k_7, k_8), the hopping rates to the 8 surrounding sites, by solving by the formula (<ref>). However, due to the intrinsically stochastic behavior of the random walks, a simple solution might result in negative values of the hopping rates. To avoid this, we performed a Bayesian inference with a model that the observed macroscopic parameters are probabilistic variables that distribute around the true values that are defined by the microscopic hopping rates, which don't take negative values (Fig. <ref> a). This was implemented with Stan <cit.> on R <cit.> with the overall hopping rate, k as a positive real number and the hopping preferences p =(p_1, p_3, p_4, p_5, p_6, p_7, p_8) as simplex variables (0≤ p_i ≤ 1, ∑ p_i =1), where the individual hopping rates k_i = k · p_i with k=∑_i k_i. As priors, uniform probabilities were used.
  
  Fig. <ref> b shows the distributions of the posterior probabilities of the hopping rates, obtained with the fifty sets of macroscopic parameters (v_x, v_y, D_x, D_y, A, B, C, E) for each Brownian motion from Fig. <ref> e as the data for the Bayesian inference. As expected, RW1 (magenta) and RW2 (green) showed clearly distinct patterns in the hopping preferences, which are consistent with the theoretical values used for generation of the random walks by simulation. For example, the (posterior) probability of the hopping rate to the forward-left, k_1, of RW1 showed a distribution near the theoretical value 0 s^-1 while that of RW2 had a peak near 2 s^-1. This indicates that our approach can properly infer the microscopic parameters to the precision levels that are sufficient to distinguish the two example cases of random walk, which would look the same if we only analyze the movements in x- and y-directions separately.
  
  
  
  
  
  If we look closer, we realize that the peak positions of the posterior probabilities of the hopping rates don't exactly match with the theoretical values. This is likely due to the probabilistic uncertainties in the observed random walks, which are influenced by the length and the number of the trajectories and by the precision of the measurements (the motion blurs and localization errors), in combination with the natural condition that the hopping rates can't be less than zero. To assess the effects of the measurement conditions, the cycle of the generation of the observed trajectories by simulation, estimation of the coefficients for each trajectory, and inference of the rates and preferences of hopping was repeated a hundred times per condition. The average of the means and 95% credible intervals of the posterior distribution of the 100 independent trials were averaged and shown in Fig. <ref> c to e and Supplemental Figures.
  
  As expected, increasing the size of data by increasing the number of trajectories sharpened the distribution of the posterior probabilities of the rates and preferences of hopping (Fig. <ref> c and Supplementary Figure 1). Importantly, their means, which showed deviations from the theoretical values when the number of input trajectories was small as mentioned above, asymptotically approached their respective theoretical values as the available trajectories increase. Similar sharpening of the distribution and approaching to the theoretical value were observed also when the length of each trajectory was increased instead (Supplementary Figure 2).
  
  The localization errors don't explicitly appear in the formulas of the CVEs (<ref>, <ref>, <ref>, <ref>, <ref>, <ref>, <ref>, and <ref>) because we took averages of infinitely many combinations of the displacements. With the finite displacement data, however, the noise terms would not be completely canceled out and affect the estimations. Indeed, increasing the amplitude of the localization errors significantly broadens the posterior probability distributions of the hopping preferences (and the hopping rates) (Fig. <ref> d and Supplementary Figure 3). For example, with the data of the 50 trajectories with 200 observations each, although the p_4 values of RW1 and RW2 were distinguishable in the presence of the localization errors up to half the size of the lattice unit, they became indistinguishable when the amplitude of the errors was increased to be the same as the lattice unit (Fig. <ref> d). This is in contrast with little impact of the extent of correlation between the x- and y-measurement errors for each observation  (Fig. <ref> e, Supplementary Figure 4). Interestingly, the perturbation by the localization errors could be overcome by increasing the size of data, either by increasing the number of observations per trajectory or by observing more trajectories (Fig. <ref> f). The difference between RW1 and RW2 that was once obscured by the increased localization errors (1 grid size) became distinguishable again with a 10-fold increase in the amount of data either by increasing the number of observations made per trajectory or by increasing the number of trajectories observed.
  
  Although the hopping rates of RW1 and RW2 were set to have identical theoretical values of v_x, v_y, D_x, and  D_y, their theoretical values of A were different (0 vs 3 s^-1). Thus, we can't exclude the possibility that the above success in discriminating between the hopping patterns of RW1 and RW2 might solely rely on the distinct values of A, irrespective of the  coefficients for the other co-moments. To test whether our approach can distinguish two random walks that have identical v_x, v_y, D_x, D_y, and A, we finally consider another random walk, RW3 (Fig. <ref> a), whose theoretical A, in addition to the theoretical values of v_x, v_y, D_x, and  D_y, is identical to that of RW1. As expected, the coefficients v_x, v_y, D_x, D_y, and  A calculated from the 50 RW3 trajectories each with 200 observations showed distributions indistinguishable from the corresponding ones from the RW1 trajectories  (Fig. <ref> b). In contrast, the distributions of the calculated C and E values of RW3 were distinct from those of RW1. As expected, the posterior probabilities of the hopping rates of RW3 closely reproduced the values set for the simulation of the trajectories, exhibiting a difference from the corresponding one of RW1 (k_1, k_4, k_6, and k_7) (Fig. <ref> c). This suggests that our approach provides reasonable inference of the hopping behaviors under a realistic setting for the single particle observation even in a case where B, C or E, the coefficients for the higher-order co-moments of the drift-adjusted x- and y-displacements, are the sole observable clues.
  
  
  
  
  
  § DISCUSSION
  
  Here we studied the anisotropic random walk on a 2D lattice and derived the relationship between the macroscopic coefficients of the time-development of the drift-adjusted displacements and the microscopic hopping rates. We then extended the covariance-based estimators of the 1D diffusion coefficient to the higher-order co-moments of the drift-adjusted displacements and used them to infer the hopping rates from the trajectory data affected by motion blur and localization errors. The versatility of this novel approach was evaluated with the trajectory data generated by simulation. We could demonstrate that our approach can distinguish 2D random walks that have exactly the same x- and y-velocities and diffusion coefficients but have distinct hopping patterns. 
  
  An important advantage of the covariance-based estimators <cit.> is that the equations to calculate them don't explicitly contain terms for localization errors nor motion blur after taking the average of infinite terms of observed x- and y-displacements combined. This assumption, of course, is not true with the real-world data of finite size. With a finite size of data, the failure in cancellation of the interfering factors as well as the intrinsic uncertainly of the random process results in an error in the estimation of the coefficients that characterize the random walk. Although it had been demonstrated that the CVE-based approach is superior to the traditional ones in the case of 1D unbiased diffusion <cit.>, it was unclear how robust a 2D version of the CVE-based approach would be since it is involved with the higher-order terms of the displacements. However, through the analysis of the simulated model cases of 2D random walk, it was confirmed that the CVE-based approach combined with Bayesian inference can reasonably estimate the anisotropic hopping patterns, with a realistic size of data when the precision of the measurement is smaller than the size of the grid. Even with a lower precision of the measurements, increasing the size of data could restore the accuracy of the estimation. For real-world data, careful optimization will be necessary, considering rather complicated and intertwined influences of the frequency of image acquisition, exposure time, the intensity of illumination <cit.>. 
  
  It has recently been reported that some kinesin-like motor proteins show a helical motion around a microtubule <cit.>. The tubular structure of a microtubule consists of 13 or 14 protofilaments, i.e., the linear arrays of α- and β-tubulin heterodimers with the 8 nm periodicity, presenting the motor-binding sites as the 2D lattice on the surface <cit.>. Interestingly, the reported pitches of the helical are intermediate between the shortest helical pitch of the lattice due to the staggered alignment of the protofilaments and the longest helical pitch observed in the 14-protofilament tubule. The observed helical motion corresponds to the protofilament switch that occurs once per ∼10 forward steps on average, implying the stochastic stepping patterns to the neighboring sites. Our approach might be applicable to estimate such patterns based on the observed helical trajectories.
  
  On a 3D lattice, there are 26 (=9+8+9) choices for stochastic hopping to a nearby site. The relation between the macroscopic coefficients and the hopping rates for a 3D lattice analogous to the 2D version (<ref>) will contain 26 hopping rates and 26 coefficients for the time-development of the combinations of (higher-order co-)moments of the drift-adjusted displacements up to the 6th order. Whether the recipe to derive the covariance-based estimators of the coefficients for the higher-order co-moments developed here for a 2D lattice (<ref>) is applicable to a 3D lattice will be a future question.
  
  
  The author would like to thank Huong T. Vu (University of Warwick, UK), Matthew Turner (University of Warwick, UK), Junichiro Yajima (University of Tokyo, Japan) and Izumu Mishima (Durham University, UK) for insightful discussions and critical reading of the manuscript.
  
  
  
  
  
  
  
  
  § PROOF OF THE EQUATION FOR THE COVARIANCE-BASED ESTIMATOR OF E
  
  Here we will prove Equation <ref>.
  Let us define notations, ΔZ(t, k, α) = Z(t+(k+α)Δ t) - Z(t + (k+α -1)Δ t) and Δε_k, α =ε_k+α+1 - ε_k+α. With these, 
  
      f( Δ Z_k+α,  Δ Z_k+β  , Δ W_k+γ, Δ W_k+λ ) 
         =   ⟨[∫_0^Δ t s(t) ΔZ(t, k, α)dt + Δε_k, α] 
          ×[∫_0^Δ t s(t') ΔZ(t', k, β)dt' + Δε_k, β] 
            ×[ ∫_0^Δ t s(t”) ΔW(t”, k, γ)dt” + Δδ_k, γ]  
          ×[ ∫_0^Δ t s(t”') ΔW(t”', k, λ)dt”' + Δδ_k, λ] 
      ⟩
                  - ⟨[∫_0^Δ t s(t) ΔZ(t, k, α)dt + Δε_k, α]  ×[∫_0^Δ t s(t') ΔZ(t', k, β)dt' + Δε_k, β] ⟩
                       ×⟨[ ∫_0^Δ t s(t”) ΔW(t”, k, γ)dt” + Δδ_k, γ] ×[ ∫_0^Δ t s(t”') ΔW(t”', k, λ)dt”' + Δδ_k, λ]  ⟩
                 - ⟨[∫_0^Δ t s(t) ΔZ(t, k, α)dt + Δε_k, α] ×[ ∫_0^Δ t s(t”) ΔW(t”, k, γ)dt” + Δδ_k, γ]  ⟩
                       ×⟨[∫_0^Δ t s(t') ΔZ(t', k, β)dt' + Δε_k, β] ×[ ∫_0^Δ t s(t”') ΔW(t”', k, λ)dt”' + Δδ_k, λ] ⟩
                  - ⟨[∫_0^Δ t s(t) ΔZ(t, k, α)dt + Δε_k, α] ×[ ∫_0^Δ t s(t”') ΔW(t”', k, λ)dt”' + Δδ_k, λ] ⟩
                       ×⟨[∫_0^Δ t s(t') ΔZ(t', k, β)dt' + Δε_k, β] ×[ ∫_0^Δ t s(t”) ΔW(t”, k, γ)dt” + Δδ_k, γ]  ⟩
            = ∫_0^Δ t∫_0^Δ t∫_0^Δ t∫_0^Δ t s(t)s(t')s(t”) s(t”') ×
                        f(ΔZ(t, k, α), ΔZ(t', k, β), ΔW(t”, k, γ) , ΔW(t”', k, λ)  ) dt dt' dt” dt ”' 
                                    +  f(Δε_k, α, Δε_k, β,  Δδ_k, γ, Δδ_k, λ)
  
  In the above transformation, we used that the true displacements and the measurement errors are independent of each other and that the ensemble means of the true displacements and those of the measurement errors are both zero. 
  
  Now we evaluate f(ΔZ(t, k, α), ΔZ(t', k, β), ΔW(t”, k, γ) , ΔW(t”', k, λ)  ) in the first term. For unbiased random walks A(s), B(t), C(u), D(v) (⟨ A(s) ⟩=⟨ B(t) ⟩=⟨ C(u) ⟩ = ⟨ D(v) ⟩ = 0, and their displacements are independent for non-overlapping time sections), we can show, by a similar calculation to (<ref>), that
  
      f(A(s), B(t), C(u), D(v)) = f(A(τ), B(τ), C(τ), D(τ)),
  
  where τ = min(s, t, u, v). Thus, with τ = min(t +(k+α -i)Δ t, t' +(k+β -i')Δ t, t” +(k+γ -i”)Δ t,  t”' +(k+λ -i”')Δ t),
  
      f(ΔZ(t, k, α), ΔZ(t', k, β), ΔW(t”, k, γ) , ΔW(t”', k, λ)  )
         =   ∑_i=0, 1∑_i'=0, 1∑_i”=0, 1∑_i”'=0, 1 (-1)^i + i' + i” + i”'×
             f(Z(t +(k+α -i)Δ t), Z(t' +(k+β -i')Δ t), W(t” +(k+γ -i”)Δ t), W(t”' +(k+λ -i”')Δ t)) 
         =   ∑_i=0, 1∑_i'=0, 1∑_i”=0, 1∑_i”'=0, 1 (-1)^i + i' + i” + i”'  f(Z(τ), Z(τ), W(τ), W(τ)) 
         =   ∑_i=0, 1∑_i'=0, 1∑_i”=0, 1∑_i”'=0, 1 (-1)^i + i' + i” + i”'  E·τ
         =    E ·∑_i=0, 1∑_i'=0, 1∑_i”=0, 1∑_i”'=0, 1 (-1)^i + i' + i” + i”'×
                 min[ +(k+α -i)Δ t, t' +(k+β -i')Δ t, t” +(k+γ -i”)Δ t,  t”' +(k+λ -i”')Δ t]
  
  Here, let us consider
  
      g(α, β, γ, λ, t, t', t”, t”') = ∑_i=0, 1∑_i'=0, 1∑_i”=0, 1∑_i”'=0, 1 (-1)^i + i' + i” + i”'×
                 min[ +(k+α -i)Δ t, t' +(k+β -i')Δ t, t” +(k+γ -i”)Δ t,  t”' +(k+λ -i”')Δ t].
  
  
  The values of 15 combinations of α, β, γ, and λ = 0  or  1 (except for α=β=γ=λ = 1) are evaluated as 
  
      g(0, 0, 0, 0, t, t', t”, t”')    =   Δ t - (t + t' + t” + t”') 
      
      + min(t, t')     +   min(t, t”)  + min(t, t”') + min(t', t”) + min(t', t”')  + min(t”, t”') 
      
      - [ min(t', t”, t”')    +   min(t, t”, t”') + min(t, t', t”') + min(t, t', t”) ] +2 ·min(t, t', t”, t”') 
      
      g(1, 0, 0, 0, t, t', t”, t”')    =   min(t', t”, t”') - min(t, t', t”, t”')
      
      g(1, 1, 0, 0, t, t', t”, t”')    =   min(t”, t”') - min(t', t”, t”') - min(t, t”, t”') + min(t, t', t”, t”')
      
      g(1, 1, 1, 0, t, t', t”, t”')    =    t”' - [ min(t, t”') + min(t', t”') + min(t”, t”', t, t', t”, t”') ] 
         +   min(t, t', t”') + min(t, t”, t”') + min(t', t”, t”')  - min(t, t', t”, t”'),
  
  etc.
  
  All the terms except for Δ t cancel out (Table. <ref>).
  
  
  
  
  Thus,
  
      ∑_α, β, γ, λ =0 or 1 
      not α=β=γ=λ= 1 g(α, β, γ, λ, t, t', t”, t”')=Δ t
  
  i.e.,
  
      ∑_α, β, γ, λ =0 or 1 
      not α=β=γ=λ= 1 f(ΔZ(t, k, α), ΔZ(t', k, β), ΔW(t”, k, γ) , ΔW(t”', k, λ)  ) = EΔ t
  
  
  The second term of f( Δ Z_k+α,  Δ Z_k+β  , Δ W_k+γ, Δ W_k+λ ) (the formula <ref>) can be decomposed into a sum of the terms with various combinations of i, i', i”, i”'.
  
      f(Δε_k, α, Δε_k, β,  Δδ_k, γ, Δδ_k, λ) 
         =     ∑_i=0, 1∑_i'=0, 1∑_i”=0, 1∑_i”'=0, 1 (-1)^i + i' + i” + i”'  f(ε_k+α+1-i, ε_k+β+1-i', δ_k+γ+1-i”, δ_k+λ+1-i”')
  
  We need to evaluate 15 combinations of α, β, γ, and λ = 0  or  1 (except for α=β=γ=λ = 1) each for 16 combinations of i, i', i”,and  i”'. Most of the 15× 16 terms are zero. As in  Table. <ref>), the other non-zero terms are cancelled out each other. 
  
  
  
  Thus, 
  
      ∑_α, β, γ, λ =0 or 1 
      not α=β=γ=λ= 1 f(Δε_k, α, Δε_k, β,  Δδ_k, γ, Δδ_k, λ)=0
  
  Finally, by plugging (<ref>) and (<ref>) into (<ref>), 
  
      ∑_α, β, γ, λ =0 or 1 
      not α=β=γ=λ= 1 f( Δ Z_k+α,  Δ Z_k+β  , Δ W_k+γ, Δ W_k+λ ) 
         =   ∫_0^Δ t∫_0^Δ t∫_0^Δ t∫_0^Δ t s(t)s(t')s(t”) s(t”') ×
                  ∑_α, β, γ, λ =0 or 1 
      not α=β=γ=λ= 1 f(ΔZ(t, k, α), ΔZ(t', k, β), ΔW(t”, k, γ) , ΔW(t”', k, λ)  ) dt dt' dt” dt ”' 
                         + ∑_α, β, γ, λ =0 or 1 
      not α=β=γ=λ= 1  f(Δε_k, α, Δε_k, β,  Δδ_k, γ, Δδ_k, λ) 
         =   ∫_0^Δ t∫_0^Δ t∫_0^Δ t∫_0^Δ t s(t)s(t')s(t”) s(t”') E·Δ t dt dt' dt” dt”' = EΔ t
  
  we obtain the formula (<ref>).
  
  
  
  
  
  
  
  
  
  