Bridging Continuous and Lattice-Based Models of Two-Dimensional Diffusion: A Systematic Approach for Estimating Transition Probabilities, Grid Size and Diffusivity

: Lattice-based models have been broadly applied in mathematical and computational modeling of biological and biomedical systems for which spatial effects are important. These discrete models commonly include diffusion of mobile constituents as a key underlying mechanism. While the direct simulation of diffusion in continuous (off-lattice) domains is possible, it is computationally intensive, particularly when multiple coupled mechanisms are involved. This study presents a systematic approach for connecting continuous models of two-dimensional diffusion with internal obstacles to discrete, lattice-based (surrogate) models of diffusion. Results from continuous model simulations on a representative domain, and over many realizations, are used to develop accurate lattice-based surrogate models by exploiting internal symmetries. Probabilities determined for the lattice-based surrogate models are also connected to theoretical diffusivities for 2D random walks on a square lattice, necessitating the calibration of a spatial grid size. This approach can facilitate the inclusion of more accurate diffusive transport models of complex media within the general framework of lattice-based models that incorporate multiple coupled mechanisms.


Introduction
Lattice-based models have seen broad application in the study of a variety of mechanisms and dynamics through simulation in environments where incorporation of spatial effects is important.For example, they have been commonly used for modeling biotransport, cell dynamics and tissue dynamics with applications in many biological and biomedical areas [1][2][3][4][5][6].One example is the cellular Potts model, in which a biological cell is defined over a region comprising multiple lattice sites [2,5,6].This lattice-based modeling framework enables the study of coupled mechanisms in biological cell dynamics and has been applied to various biomedical problems, including epidermal biology, cancer, angiogenesis and vasculogenesis [1,4].Such discrete, lattice-based models commonly involve diffusion (e.g., of mobile solutes) as a key transport mechanism.Alternatively, agent-based models [7][8][9] often restrict movement to a lattice, with random walks commonly used as a simple or convenient approach to advance individual agent locations at each time increment.In such cases, diffusion of nutrients and cell-synthesized solutes is an important coupled mechanism.
As computing power has advanced, direct simulation in continuous-domain (offlattice) models is a useful tool for investigating systems where particles or agents interact with a background medium or with each other.However, the direct approach can become tedious or computationally intensive in situations involving several coupled mechanisms or complex geometry or in the presence of many particles or agents.The ability to systematically connect off-lattice (continuous) models to lattice-based (discrete) models, accounting for a few key underlying model features or mechanisms, has the potential to improve the overall accuracy of discrete models.Moreover, such discrete models can also be connected to estimates of theoretical diffusivities for random walks on a two-dimensional (2D) square lattice.An appropriate spatial length scale (grid size) must be determined for the discrete model, especially when diffusion is one of several coupled mechanisms in a 2D lattice-based model.
In this study, we evaluate the extent to which 2D lattice-based surrogate models can preserve measures of particle diffusivity observed in direct simulation of continuous models based on 2D random walks with internal obstacles.Our overall approach is illustrated in Figure 1.Results from continuous model simulations (Figure 1A) on a 2D representative volume element (RVE), and over many realizations, are used to develop lattice-based surrogate models (Figure 1B) that reproduce key features of the continuous model.This is accomplished by exploiting internal symmetries in the continuous model representation.Probabilities in the lattice-based surrogate models are also connected to theoretical diffusivities for 2D random walks on a square lattice (Figure 1C).This final step requires calibration of the discrete length scale (grid size) in the square lattice surrogate model.The random walker is prescribed a fixed particle radius, and obstacle radii are varied in the context of two different obstacle configurations, each having distinct internal symmetries.Particle dynamics observed for the continuous model simulations are used to develop probabilistic lattice-based surrogate models for transitions across RVE subdomains.The subdomains are chosen to preserve internal symmetries for each of the two obstacle arrangements.A primary advantage of this approach is the capability of the surrogate models to bypass the need to directly represent the obstacles in the simulation domain.
The well-known linear relationship between the mean squared particle displacement and time for the continuous model [10][11][12] is used to estimate and compare diffusivities for both the continuous and lattice-based surrogate models.This comparison motivates the introduction of a simulation parameter (commitment index) that is calibrated to minimize adverse effects of noise on achieving stationarity of transition probabilities in the latticebased surrogate model.Theoretical derivations for isotropic and anisotropic diffusion on a 2D lattice are used to connect the transition probabilities in the lattice-based surrogate models to theoretical estimates of the diffusivity for a square lattice surrogate model.Relations enforcing consistency between diffusivity estimates, based on continuous model simulations and analytical results on the lattice, are used to estimate the effective spacing (grid size) in a (square) lattice-based model for both obstacle configurations.
This paper is organized as follows.In Section 2, effective diffusivity is first outlined in terms of the mean squared displacement of a single particle undergoing a 2D random walk in a continuous domain or on a 2D lattice.Section 3 presents our equally spaced obstacles model (Section 3.2) and our multisize obstacles model (Section 3.3).Within each model, quantities are defined to determine probabilities connecting the continuous model with its lattice-based surrogate model.In Section 4, methods for calibrating the lattice-based surrogate models based on the continuous model simulation results are first presented (Section 4.1).A framework is then presented to connect the lattice-based surrogate models to square lattice surrogate models (Section 4.2), necessitating the calibration of the square lattice grid size parameter.Our results are presented for both models in Section 5, including the calculation of all model probabilities and calibration of other model parameters, followed by a brief discussion and conclusions in Section 6. Supplementary derivations are presented in Appendices A and B.

Effective Diffusivity
We first summarize approaches for estimating the effective diffusivity of a particle undergoing 2D random walks based on observations of its mean squared displacement.

Effective Diffusivity in a Continuous Domain
Normal diffusion describes diffusion via random walks of a particle in a continuous domain.The mean-squared displacement ⟨r 2 ⟩, relative to a fixed starting location at t = 0, is calculated by averaging displacement trajectories over a large number of realizations.In prior studies of 2D single particle random walks, it was shown that ⟨r 2 ⟩ grows linearly with time (t) after a sufficiently large number of time steps have passed [10][11][12].The general formula for normal diffusion is written in terms of an effective diffusivity (κ) as [11][12][13] where d is the dimension of the medium.In our study, all simulation results are sampled in a stationary time regime where this observed relationship is linear; the effective diffusivity is estimated by rearranging (1) to obtain where N is the number of time steps (of size ∆t) in the simulation.

Effective Diffusivity on a 2D Lattice
To develop our lattice-based surrogate models, we estimate the diffusivity analogous to (2) but based on analyzing ⟨r 2 ⟩ on the lattice.In this section, we demonstrate, via an analytical derivation, that a particle undergoing random walks on a 2D lattice exhibits normal diffusion, to leading order.Overall, our derivations of mean squared displacements and theoretical diffusivities for random walks on 2D lattices are based on analytical approaches outlined in Holmes [14].
We can track ⟨r 2 ⟩ systematically through a set of transition probabilities (Figure 2).Here, p x and p y denote the probability of moving, at the next time step, to adjacent horizontal (by ∆x) or vertical (by ∆y) lattice sites, respectively.The particle remains in its current location with probability 1 − p x − p y .In this section, we denote the mean squared displacement as ⟨r 2 K (N)⟩, with K referring to the number of realizations of a single particle random walk on the 2D lattice and N denoting a particular time step.

Lattice Representation of Mean Squared Displacement
The mean squared displacement after N time steps and for K realizations is where x i (N) and y i (N) are the horizontal and vertical particle locations, respectively.The mean squared displacement is then ⟨r 2 ⟩ = lim K→∞ ⟨r 2 K (N)⟩.For brevity, let x N i and y N i denote the horizontal and vertical positions of the i th particle at time step N, respectively, (the full notation is utilized in the Appendices).We write these positions as where h N xi and h N yi are the possible displacements with associated step sizes of ∆x and ∆y, respectively.In any one time step ∆t, the five possible transitions are where cases (i) (horizontal transition) and (ii) (vertical transition) occur with equal probability, and case (iii) represents no transition.Substituting (4) into (3), ⟨r 2 ⟩ can be written in recursive form as In the next subsection, we briefly discuss limits used to simplify (6).

Useful Limits for Simplifying Equation (6)
We use two main limit results in simplifying Equation (6).Direct evaluation of these limits is summarized in Appendices A and B, resulting in the following values associated with case (i) (horizontal transition) above: Analogous results also hold for the vertical transition case (ii).

Reduced Form of Mean Squared Displacement
Using (7), we rewrite (6) as to simplify (8) to the relation, Lastly, the recurrence property enables (9) to be simplified to On a square lattice (∆x = ∆y), (10) can be reduced to the linear relation where we denote q as the probability of staying in the current subdomain, i.e., since q + p x + p y = 1.Equation (11) will serve as the basis for connecting the square lattice surrogate models with the lattice-based surrogate models obtained from the continuous model simulations (see Section 4.3).

Models
We consider direct simulation of 2D random walks on a continuous domain.Results from these simulations are used to determine transition probabilities in the lattice-based surrogate model.The number of such independent probabilities, as well as their values, change for different obstacle arrangements.We consider a single particle undergoing random walks in a 2D representative volume element (RVE).Within the RVE, we place a small number of rigid obstacles with which the particle collides; the collisions are idealized to be perfectly elastic with the angle of reflection equal to the angle of incidence.The RVE is assumed to have periodic boundary conditions, further capturing internal symmetries in the background medium.Two different obstacle arrangements are considered (Figure 3): (i) equally spaced obstacles (Figure 3a) and (ii) multisize obstacles (Figure 3b).While both cases model inhomogeneous, anisotropic materials, case (i) exhibits a greater degree of (rotational) symmetry as compared to case (ii).The second case also serves as a simple analogue of a composite with two different obstacle sizes.

Continuous Model Simulation Details
We consider symmetries inherent to each case that are evident when the RVE is partitioned into nine equally sized subdomains for case (i) (Figure 3a) and four equally sized subdomains in case (ii) (Figure 3b).Each (square) subdomain is prescribed a length L = 10; thus, to capture internal symmetries, the RVE is 30 × 30 for case (i) and 20 × 20 for case (ii).The state of the system at time t is the subdomain in which the particle is located (1-9 in case (i) and 1-4 in case (ii), Figure 3).State changes, which are vertical or horizontal transitions between neighboring subdomains, are tracked to determine transition probabilities for this surrogate model (see Sections 3.2.2 and 3.3.2).
Additional assumptions are as follows: (a) our single particle has prescribed radius R part and moves with a speed v = 8; (b) at t = 0, the particle is placed at the center of subdomain 5 for the equally spaced obstacles model and at the center of subdomain 3 for the multisize obstacles model; (c) at each time step of prescribed temporal step size ∆t = 0.05, a new direction is chosen, via uniform sampling, from θ ∈ [0, 2π], and the particle is advanced a distance of ∆x = v∆t in this direction; (d) collisions with the obstacles are perfectly elastic, i.e., the particle bounces off the obstacle at the (reflected) incidence angle relative to the local tangent line at the point of collision.
We simulate many realizations over N t = 1 × 10 7 time steps to generate plots of the mean-squared displacement ⟨r 2 ⟩ versus time.For all simulations considered in this study, 2000 realizations on a high-performance cluster computer were sufficient for the slopes of these linear curves to reach stationary values.

Continuous Model
For the equally spaced obstacles model, the circular obstacles (radius R obs ) are centered on the subdomain intersection points.The particle radius is fixed at R part = 1.0, and the obstacle radius is varied from 0.5 to 3.0 by increments of 0.5.This range was chosen to enable particle movement between obstacles; an obstacle radius greater than 3.0 inhibits transitions between subdomains.

Lattice-Based Surrogate Model
Figure 4 illustrates the construction of the 2D lattice-based surrogate model for the equally spaced obstacles case (Figure 3a).A periodic extension of the RVE is shown, with a primary transition from the center subdomain illustrated.Due to extensive symmetries in this case, all internal transitions across boundary walls between subdomains are equivalent.The probability of diagonal transitions is zero due to the presence of obstacles at the subdomain corners.Hence, the lattice model has only two distinct transition probabilities, that the particle (1) stays in its current subdomain (p E 1 ) or ( 2) moves to a neighboring subdomain in a cardinal direction (p E 2 ).The corresponding state transitions are illustrated for subdomain 5 in Figure 4b.Since transition probabilities within a subdomain must sum to 1, only one surrogate model probability needs to be calculated from the results of continuous model simulations.Results of the continuous model simulations (Figure 4a) are tracked via the time course η(t k ) (t k = k∆t, k = 0, . . ., N), where N is the total number of time steps.η(t k ) records the state of a particle as the index of the ith subdomain (i = 1, . . ., 9) at time t k if its center is located in that subdomain at time t k .A 9 × 9 matrix S (k) is used to track all transitions (including remaining in the same subdomain) up to the current time step (S (0) is the zero matrix).S (k) is updated according to the formula S (k) is then used to build a vector ⃗ p E of transition probabilities.Within each time step, a vector ⃗ T E i (i = 1, . . ., 9) is created to count transitions that occurred for each subdomain.In particular, the rows of S (k) are used to create vector ⃗ T E i by grouping transition sums by transition type; since the equally spaced obstacles model has two transition types, ⃗ T E i has dimensions of 1 × 2. At the current time t k , we construct ⃗ T E i from S (k) as follows, Due to the symmetries in the equally spaced obstacles case, we create the transition probability vector ⃗ p E via the average:

Continuous Model
For the multisize obstacles model, the RVE has one set of obstacles with a fixed obstacle radius R obs 1 and a second set of obstacles of radius R obs 2 that is varied across cases.Symmetry for this model is preserved by partitioning the RVE using a 2 × 2 grid of subdomains (Figure 3b).We view the domain as having "columns" of obstacles; obstacles in the same column have identical radii.The particle radius and first obstacle radius are fixed at values of R part = 1.0 and R obs 1 = 1.0, respectively.The second obstacle radius (R obs 2 ) is then varied between 0.5 to 3.0 in increments of 0.5.

Lattice-Based Surrogate Model
Figure 5 illustrates the construction of the 2D lattice-based surrogate model for the multisize obstacles case.A periodic extension of the RVE is shown, with primary transitions from subdomain 1 (top-left corner) illustrated.Due to inherent symmetries, all subdomains have the same number of distinct transition types.There are four types of transition with corresponding probabilities (Figure 5b).In type 1, the particle stays in its current subdomain (p M 1 ).In types 2-4, it transitions to a neighboring subdomain (2,3 or 4) having two obstacles.We delineate the three latter types into probabilities as follows: the two obstacles are of different sizes (R obs 1 , R obs 2 ) along the boundary (p M 2 ); the two obstacles are the same size R obs 2 along the boundary (p M 3 ); or the two obstacles are the same size R obs 1 along the boundary (p M 4 ).Again, the probabilities must sum to one; the three surrogate model probabilities need to be calculated from the results of the continuous model simulations.Again, the time course η(t k ) records the state of the particle as the index of the ith subdomain (i = 1, ..., 4) at time t k if its center is located in that subdomain at time t k .For the multisize obstacles case, we differentiate between a horizontal transition within the RVE versus a horizontal transition into an extended domain by introducing an additional time course 1 if a horizontal transition occurred into an extended domain 0 otherwise (includes all vertical transitions). .
In this case, the transition matrix S (k) (for time steps k = 1, . . ., N) is a 4 × 5 (augmented) matrix; the fifth column tracks horizontal transitions occurring when the particle moves into an extended domain.For a large majority of time steps, ζ(t k ) = 0 and S (k) are updated in a similar manner as in the equally spaced obstacles case via Equation (12) (Section 3.2.2).Alternatively, when ζ(t k ) = 1, only the new column in the augmented matrix is updated via the formula and S (k) is then used to build the vector ⃗ p M of five transition probabilities.Within each time step and for each of the four subdomains, a vector ⃗ T M i (i = 1, 2, 3, 4) is created to count transitions, according to the following rules: (i) the first component counts the number of times a particle stayed in its current subdomain; (ii) the second component counts the number of times a particle moved vertically; (iii) the third component counts the number of times a particle moved horizontally but remained in the original RVE; (iv) the fourth component counts the number of times a particle moved horizontally but moved to the extended domain.
Due to the symmetries in this case, we create the transition probability vector ⃗ p M via the average:

Model Calibration and Square Lattice Surrogate Models 4.1. Commitment Index
For random walks in our continuous model, a particle may cross, return and then recross a subdomain boundary several times before "committing" to a new subdomain.This leads to inaccuracy in the calculation of transition probabilities as stationary quantities in the lattice-based surrogate models.Consequently, we introduce a commitment index, M C , to indicate the number of successive time steps after which a particle is deemed to commit to a new subdomain.The role of this parameter is illustrated in Figure 6 for the case M C = 4.While the value of M C requires calibration via direct comparison of simulated results for mean-squared displacement in the continuous and lattice-based surrogate models, we now explain its narrow range of variation.(i) Consider a particle, committed to a subdomain, that has just crossed the boundary into a neighboring subdomain (first-step location).A (black) ring indicates all possible particle locations at the next step, including some locations that result in recrossing the same boundary (Figure 6a).(ii) We next illustrate a few possible steps relative to the particle's first-step location, thus indicating possible second-step locations.A (red) ring of possible 3rd-step locations is shown around each such second-step location (Figure 6b).(iii) Choosing one of the possible second-step locations that is in the new subdomain (gray ring), the resulting set of possible third-step locations is shown in red (Figure 6c).(iv) Lastly, the larger the number of successive steps for which a particle remains in the new subdomain, the smaller the probability that it will recross back into the subdomain to which it was previously committed; in Figure 6d, blue rings depict fourth-step locations.
For each of our two obstacle arrangements, M C is calibrated by constructing transition probabilities from continuous model simulations over a range of M C values and by comparing mean-squared displacement for the continuous and lattice-based surrogate models to identify the case with the best agreement (see Section 5.1).

Diffusivities for the Square Lattice Surrogate Models
Recall from Section 2 the following estimates of the effective diffusivity (e.g., (2) with d = 2) for the continuous model and the surrogate model, respectively, , and In Section 4.2, we use 2D Taylor series expansions to derive the (continuum) diffusion equation and the associated theoretical diffusivities on a lattice for both obstacle arrangements.

Equally Spaced Obstacles Model
Consider a 2D rectangular lattice with horizontal spacing h(= ∆x) and vertical spacing k(= ∆y).Assuming a fixed volume and a uniform particle mass, let u(x, y, t) be the density of particles and let q be the probability that a particle stays in its current subdomain; (1 − q) is the probability of transitioning to a neighboring lattice point.The resulting particle density is given by u(x, y, t + ∆t) = qu(x, y, t) Consider the Taylor series expansions about the point P(x, y, t), u(x, y, t Using these series expansions, (19) becomes Assuming a 2D square lattice (∆x = ∆y) and dividing by ∆t, we can rewrite (23) as resulting in the approximate equation (assuming Recognizing (25) as the 2D diffusion equation, the theoretical diffusivity for the equal obstacles model (D) is given in terms of square lattice model properties by

Multisize Obstacles Model
In the multisize obstacles model, the particle has different diffusivities in the vertical and horizontal directions.Yet, the large majority of lattice-based models in computational biology employ a square lattice.Thus, while our derivation below allows for distinct values of ∆x and ∆y, we ultimately set these quantities to be equal in order to reconcile the two viewpoints.
We follow the same procedure as in Section 4.2.1, but for this case, note the following: (i) q 1 is the probability that the particle stays at its current lattice point; (ii) q x is the probability that the particle moves to a neighboring lattice point in the horizontal direction; (iii) q y is the probability that the particle moves to a neighboring lattice point in the vertical direction.Note that q 1 + q x + q y = 1.We can then write the following relation for the particle density: Via the Taylor series expansions in (20)-( 22), the above equation can be expanded and simplified to + q y 2 u yy (∆y) 2 + O((∆y) 4 ) . (28) Since q 1 + q x + q y = 1, upon dividing by ∆t we can rewrite (28) as Imposing our restriction to a uniform 2D lattice (∆x = ∆y), the simplified diffusion equation for the multisize obstacles case is thus (assuming q x (∆x) 4 ≪ ∆t and q y (∆x) 4 ≪ ∆t) Hence, when a 2D square lattice is used as a surrogate model for the multisize obstacles case, there are two theoretical diffusivities that we refer to as the horizontal diffusivity (D H ) and the vertical diffusivity (D V ), which are given by We note that in the case q x = q y , the diffusion equation for the multisize obstacles model (30) reduces to the diffusion equation for the equally spaced obstacles model (25).

Relating the Diffusivities
We now relate our effective diffusivities (κ cont , κ surr ) in (18) to the theoretical diffusivities for each of the two obstacle models, i.e., (26) and (31).

Equally Spaced Obstacles Model
The effective diffusivity for the surrogate model κ surr is first used to approximate the effective diffusivity for the continuous model κ cont .Specifically, slopes of the mean squared displacement response are matched by calibrating the commitment index M C in the surrogate model (see Section 4.1).
Substituting ⟨r 2 ⟩ from (11) for the surrogate model into the second relation of (18 ), we obtain the relations, Since ∆t is the same across all representations, it remains to calibrate ∆x, i.e., the effective spatial grid size for the square lattice model with theoretical diffusivity D.
Although L = 10 is a characteristic length scale for both the continuous and (surrogate) lattice models, it is not an appropriate choice for ∆x in (32).A particle at an arbitrary location in a subdomain has many different net distances it can travel before crossing into and committing to a neighboring subdomain.We determine the spatial step size via the representation ∆x = L α , where α > 1 will be calibrated based on our numerical simulations.Substituting ∆x = L α into (32), we obtain Hence, the ratio of the two diffusivities can be used to approximate α via the formula Using this coordinated approach, we ensure that all diffusivity estimates obtained via simulations using the continuous model and the surrogate model are consistent with theoretical estimates for 2D random walks on a lattice based on the choice ∆x = L α .

Multisize Obstacles Model
Considering the horizontal and vertical effective diffusivities independently, and using (2) with d = 1, we can define κ H and κ V as where ⟨r 2 H ⟩ and ⟨r 2 V ⟩ are the horizontal and vertical mean squared displacements.In the multisize obstacles model, we identified horizontal and vertical components of diffusivity in terms of the corresponding probabilities (see (31)).
In contrast to the equally spaced obstacles model, here we do not directly identify the probabilities q x (of a horizontal transition) and q y (of a vertical transition) in the square lattice surrogate model.In particular, q x is not explicitly identified from our simulations, since our results delineate (1) a horizontal transition in the original RVE and (2) a horizontal transition to the extended domain.We, instead, relate the following ratios: Note that in the case where all obstacles are the same size, this yields D H ≈ D V .For the lattice-based surrogate model, the four probabilities p 1 , p 2 , p 3 , p 4 correspond to the transitions depicted in Figure 5b.Recall that q 1 + q x + q y = 1 in the square lattice surrogate model and p 1 + 2p 2 + p 3 + p 4 = 1 in the lattice-based surrogate model.Equating these two relations and assuming that q 1 = p 1 and that 2p 2 = q y , we obtain the reduced equation Assuming κ H ≈ D H and κ V ≈ D V , we can use simulation results to calibrate ∆x using each of the two relations in (31) (see Section 5.2.2).Moreover, we note that the ratio of diffusivities in (36) is independent of the square lattice grid size.

Results
For both the equally spaced obstacles model (Section 3.2) and the multisize obstacles model (Section 3.3), we use results of the continuous model simulations to determine transition probabilities for the two surrogate models, as described in Section 3.2.2 and Section 3.3.2,respectively.As outlined in Section 3.1, we simulate a single particle (radius R part ) undergoing a 2D random walk (speed v = 8) with time step ∆t = 0.05.At each step, the direction is sampled from a uniform distribution and all collisions with internal obstacles are perfectly elastic.All continuous model simulations were carried out on a highperformance cluster computer for N = 2000 realizations.Once transition probabilities were determined, the corresponding lattice-based surrogate model simulations were carried out on a laptop computer (MacBook Pro) in under 2.5 h.By contrast, the continuous model simulations were completed at a rate of 400 realizations per 27 h, distributed over eight processors (50 realizations per processor), for a specified set of model parameters.A total of N t = 1 × 10 7 time steps were sufficient for all quantities of interest to exhibit stationary responses.

Equally Spaced Obstacles Case
Obstacle radii: Continuous model simulations were performed on the RVE for the equally spaced obstacles model (see Section 3.2.1 and Figure 3a).The particle radius was first fixed at R part = 1.0, and the obstacle radius (R obs ) was varied from 0.5 to 3.0 by increments of 0.5.Stationarity and linearity of ⟨r 2 ⟩ versus time (Section 2.1) for an increasing number of realizations is illustrated, for a representative case (R part = 1.0,R obs = 1.0), in Figure 7. Determining transition probabilities: Based on the approach illustrated in Figure 4, continuous model results were used to determine transition probabilities for the equally spaced obstacles model via (14).These probabilities are illustrated for transitions from subdomain 6 over the time course of one realization in Figure 8b-d; rescaling the range more effectively illustrates the stay and leave probabilities (Figure 8c,d, respectively).This figure also demonstrates the stationary nature of these probabilities at later times.In determining the final stay (p E (1)) and leave (p E (2)) transition probabilities for the latticebased surrogate model, values of each type were pooled together (see ( 14)) and then averaged across all realizations.Note that our results are consistent with the improbability of diagonal transitions (Figure 8d, blue curve).Commitment index calibration: Lattice-based surrogate models were built using the transition probabilities obtained from the continuous model (see (14), Section 3.2.2) by first calibrating the commitment index M C (Section 4.1).M C was varied from 3 to 8 and calibrated by comparing the continuous and surrogate plots of ⟨r 2 ⟩ versus time.This procedure was repeated for varying obstacle radii, with the particle radius fixed at R part = 1.0, illustrated in Figure 9 for R part = 1.0 and R obs = 2.0.Here, M C = 5 yields the best result; 2000 realizations were sufficient to observe stationary values sufficient to calibrate M C (Figure 9-right).Probabilities and commitment index values are shown in Table 1.We observe that the sum of the first probability and four times the second probability is very close to 1 (based on symmetries in the RVE); the largest discrepancy was observed to be 3.6 × 10 −6 .We also observe that the probability of staying in the current subdomain increases with the obstacle radius.Comparisons of slopes for ⟨r 2 ⟩ versus time for the continuous and surrogate models are compared in Table 2 and Figure 10.We observe differences in slopes ranging from 0.004% to roughly 4%.

Multisize Obstacles Case
Obstacle radii: Continuous model simulations were also performed for the multisize obstacles model (Figure 5) on an RVE partitioned into four subdomains as described in Section 3.3.1.The particle radius was fixed at R part = 1.0, with one obstacle radius fixed at R obs 1 = 1.0 and the other obstacle radius (R obs 2 ) varied from 0.5 to 3.0 in increments of 0.5.Stationarity and linearity of ⟨r 2 ⟩ versus time is illustrated, for a representative case (R part = 1.0,R obs 1 = 1.0 and R obs 2 = 1.0), in Figure 11.Determining transition probabilities: Based on the approach illustrated in Figure 5, continuous model results were used to determine transition probabilities for the multisize obstacles model using (17).
Commitment index calibration: Lattice-based surrogate models were built using the transition probabilities obtained from the continuous model (Section 3.  We observe the manner in which the slope decreases with increasing obstacle radius R obs 2 , which differs from that observed in the equally spaced obstacles case (Figure 10).A summary of the transition probabilities and values of M C needed to build the lattice-based surrogate model for all cases considered is provided in Table 3.Note that the sum of the probabilities in each row is 1 when the second probability is weighted by a factor of 2 due to symmetries in the RVE.
Comparisons of slopes for the relation between ⟨r 2 ⟩, including the best values for M C , are shown for the continuous and surrogate models in Table 4 and Figure 12-right.These results exhibit a difference in slopes ranging from 0.6% to roughly 7%.

Estimating Diffusivities
We use the continuous and lattice-based surrogate model results to estimate effective and theoretical diffusivities (see Figure 1).

Equally Spaced Obstacles Case
We estimate the effective diffusivities from simulated results for the mean squared displacement in both the continuous (κ cont ) and surrogate (κ surr ) models via (18) and for each obstacle radius (Figure 13).
Time averages for both effective diffusivities over the last 100,000 values are shown in Table 5, along with the commitment index M C and the stay probability q.The continuous and surrogate diffusivities exhibit a correlation coefficient of R = 0.9784 and percent differences between 0.02% and roughly 6%.Estimated diffusivities for the surrogate models were used to compute values of α 2 and ∆x(= L α ), for each case.Comparisons of the effective diffusivities to the resulting theoretical diffusivities, determined via (33), are shown in Table 6.For comparison, the theoretical estimate of the diffusivity for the case L = 10, corresponding to α 2 = 1, is also shown in the table.The estimated values of both α 2 (mean 5.0466, st.dev.0.0888) and ∆x (mean 4.4519, st.dev.0.0397) exhibited little variation or correlation (R = 0.5532) across the range of obstacle sizes considered.This finding suggests that the coefficient α 2 can be approximated as constant across the range of obstacle radii considered.
The 2D lattice-based surrogate model for the equally spaced obstacles case is thus fully represented by the set of parameters q, p(= 1 − q), ∆x and α shown in Table 6 for each obstacle radius.We also estimated effective diffusivities from simulated results for the continuous and surrogate models via (18) and for each obstacle radius combination (Figure 14).Time averages for each of these effective diffusivities over the last 100,000 values are shown in Table 7 along with the commitment index M C and the stay probability q.The continuous and surrogate diffusivities exhibit a correlation coefficient of R = 0.8503 and differences ranging between roughly 1.1% and 6.8%.
Here, ∆x M denotes the effective grid size for the multisize obstacles case, and its value remains to be calibrated.Note, however, that the ratio of these two diffusivities is independent of this parameter.Consequently, we first examined this ratio using (36) by evaluating the relevant mean squared displacements (Figure 15) and their ratio (Figure 16).The multisize obstacles model has the property that the RVE is more open in either the horizontal (R obs 2 < 1) or vertical (R obs 2 > 1) direction, depending on the value of R obs 2 .This effect is subtle for R obs 2 < 1 due to the large amount of available space for particle diffusion (Figures 15a and 16a).Significant effects on the ratio of diffusivities are observed for larger values in the regime R obs 2 > 1, where available space in the horizontal direction becomes more limited (Figures 15e,f and 16e,f).An increased time to achieve a stationary response is also observed for these cases.The effective grid size ∆x M can be computed in one of two ways, by approximating D H and D V by κ H and κ V , respectively, and solving the corresponding relation in (38) for ∆x M (Table 8).The values of ∆x M computed in this manner are 4.2268 ± 0.1444 (mean ± SD) using κ H and 4.0684 ± 0.2099 using κ V .A single-factor ANOVA analysis results in p-values of 0.159 for R obs 2 between 0.5 and 3.0, 0.226 for R obs 2 between 0.5 and 2.5 and 0.394 for R obs 2 between 0.5 and 2.0.These findings suggest that ∆x M can be approximated as constant, with some caution needed for the case with the largest obstacle sizes.The 2D lattice-based surrogate model for the multisize obstacles case is thus fully represented by the set of parameters q x , q y , q 1 (= 1 − q x − q y ) and ∆x M , shown for each obstacle radius in Table 8.In reality, using an average of the two calculated values of ∆x M is recommended.

Discussion and Conclusions
In this study, we developed a systematic approach (outlined in Figure 1) for bridging two-dimensional continuous models of diffusive transport in the presence of rigid internal obstacles with perfectly elastic collisions with two types of surrogate models.The first type is a lattice-based surrogate model that exploits symmetries in the RVE of the continuous domain used for direct numerical simulation.The second type is a square lattice surrogate model that can enable the accurate estimation of diffusivities to be subsequently used in more general lattice-based modeling frameworks where diffusion is one of several important coupled mechanisms.Here, diffusion is often the fastest spatiotemporal mechanism.Our systematic approach can be used, first, to establish the discrete model probabilities and grid size for building 2D lattice models in applications that couple diffusion with additional mechanisms in a background medium having specific geometric properties.Conversely, in other applications, the ability to calibrate the grid size such that the associated probabilities are consistent with an observed or prescribed diffusivity may also be beneficial.
While the level of noise and error increases in the more complex (multisize obstacles) model considered in this study, it remains in a range that is comparable to the level of noise or variation typical for 2D lattice-based simulations with underlying stochastic or probabilistic mechanisms.Our approach is versatile in the sense that it can be extended to other geometric configurations of obstacles in an RVE exhibiting several degrees of symmetry.Similarly, the techniques used in our overall approach can be extended to three-dimensional lattice-based models for materials with internal geometric symmetries.Nevertheless, it should be noted that the computational time needed to achieve stationarity of the underlying continuous and lattice-based surrogate model probabilities can increase substantially with increasing model complexity.Consequently, the extent to which computation time is practical (versus prohibitive) will depend on the degree of internal symmetry in the background medium as well as availability of ultra-high-performance computing resources.
For the equally spaced obstacles problem, moving left or right is equally likely.Hence, for very large values of K, it follows that k x2 K = p x 2 and that k x3 K = p x 2 .We then write lim Similarly, in the vertical direction, by analogous steps we obtain the result lim which corresponds to the first relation in (7).We can again follow a similar line of reasoning to obtain a similar result in the vertical direction.

Appendix B. Derivation of Type 2 Limits
Consider a particle's possible horizontal locations after one time step, x i (N) = x i (N − 1) + h xi (N), (A16) We use induction to derive the second equation in (7).
When N = 0, we obtain lim

Figure 1 .
Figure 1.Illustration of the approach used for bridging continuous and lattice-based models in this study.

Figure 2 .
Figure 2. Illustration of a particle undergoing a random walk on a 2D lattice.The particle moves one step (size ∆x) in the horizontal direction with probability p x (left) and one step (size ∆y) in the vertical direction with probability p y (right).The probability of staying in place is 1 − p x − p y .

Figure 3 .
Sample domains (RVEs) for continuous model simulations (R part = 1.0), with numbered subdomains: (a) equally spaced obstacles model (with R obs = 2.0); (b) multisize obstacles model (with R obs 1 = 1.0 and R obs 2 = 2.0).Obstacles are shown as black circles and the particle undergoing random walks is shown as a red circle.

Figure 4 .
Building a lattice-based surrogate model for the continuous model in the equally spaced obstacles case: (a) nine subdomains are used to exploit symmetry of the obstacle arrangement.and transitions between subdomains are tracked.(b) Lattice-based surrogate representation is shown, with red circles denoting corresponding subdomains in the continuous model.

Figure 5 .
Building a lattice-based surrogate model for the continuous model in the multisize obstacles case: (a) four subdomains are used to exploit symmetries of the obstacle arrangement, and transitions between subdomains are tracked.Smaller green circles indicate obstacles of fixed radius R obs 1 = 1.0, and larger blue circles indicate obstacles with varying obstacle radius R obs 2 ; (b) lattice-based surrogate representation with red circles denoting corresponding subdomains in the continuous model.

Figure 6 .
Figure 6.Illustration of the role of the commitment index (calibration) parameter M C : (a) The particle (solid black circle) has just crossed the boundary between subdomains (dashed gray line) to reach its 1st-step location.Possible 2nd-step locations are indicated by the solid black ring.(b) Three possible 2nd-step locations are indicated by solid red circles.Some possible next steps (3rd-step locations) are shown with red rings.(c) A red ring illustrates the particle's remaining possible 3rd-step locations when it has remained in the current subdomain at the prior step.(d) Possible 4th-step locations are indicated with blue rings.

Figure 7 .
Figure 7. Mean squared displacement versus time for the equally spaced obstacles model.Effects of an increasing number of realizations are shown for the case R part = 1.0 and R obs = 1.0:(a) curves; (b) linear fits to the curves.

Figure 8 .
Figure 8. Transition probabilities for building the lattice-based surrogate model from the equally spaced obstacles (continuous) model.Transitions from subdomain 6 are illustrated: (a) simulation domain showing the 9 subdomains; (b) plot of all transition probabilities (here, only the stay probability is discernible); (c) magnified plot to illustrate the stay probability; (d) magnified plot to illustrate the leave probabilities.

Figure 9 .
Figure 9. Mean squared displacement and its slope for the equally spaced obstacles model.Results are shown for particle radius R part = 1.0 and obstacle radius R obs = 2.0: (left) ⟨r 2 ⟩ versus time for continuous model and resulting surrogate model as the commitment index M c is varied; (right) slope as a function of obstacle radius for an increasing number of realizations.

Figure 10 .
Figure 10.Slope of ⟨r 2 ⟩ versus time for equally spaced obstacles (continuous) model and the latticebased surrogate model with best M c values.Results are shown for particle radius R part = 1.0 and varying obstacle radius R obs .

Figure 11 .
Figure 11.Mean squared displacement versus time for the multisize obstacles model.Effects of an increasing number of realizations are shown for the case R part = 1.0,R obs 1 = 1.0 and R obs 2 = 1.5:(a) curves; (b) linear fits to the curves.
3.2) by first calibrating the commitment index M C (Section 4.1).For each case considered, M C was varied from 3 to 8 and calibrated by comparing the continuous and surrogate plots of ⟨r 2 ⟩ versus time.A sample plot illustrates increasing stationarity of the slope of ⟨r 2 ⟩ versus obstacle radius as the number of realizations increases (Figure 12-left).

Figure 12 .
Figure 12.Slope of the mean squared displacement for the multisize obstacles model: (left) variation with obstacle radius as the number of realizations is increased; (right) comparison of the continuous and lattice-based surrogate models for the best M C values.

Table 1 .
Transition probabilities in the lattice-based surrogate model for the equally spaced obstacles (continuous) model for the best values of M c with particle radius R part = 1.0.

Table 2 .
Slope values for ⟨r 2 ⟩ vs. t in the equally spaced obstacles (continuous) model (R part = 1.0) and the lattice-based surrogate model.The best values of M C are also shown.

Table 4 .
Slope values for ⟨r 2 ⟩ vs. t in the multisize obstacles (continuous) model (R part = 1.0,R obs 1 = 1.0) and the lattice-based surrogate model.The best values of M C are also shown.

Table 5 .
Effective diffusivities in the equally spaced obstacles case for the continuous model (κ cont ) and the lattice-based surrogate model (κ surr ).The stay probability (q) and best commitment index (M C ) are also shown for each obstacle radius (R obs ).

Table 6 .
Theoretical diffusivity (D) and grid size (∆x) in the equally spaced obstacles case.The stay probability (q), best commitment index (M C ) and theoretical diffusivity with L = 10 (D| ∆x=10 ) are also shown for each obstacle radius (R obs ).

Table 7 .
Effective diffusivities in the multisize obstacles case for the continuous model (κ cont ) and the lattice-based surrogate model (κ surr ).The resulting stay probability (q 1 ) in the square lattice surrogate model and the best commitment index (M C ) are also shown for each obstacle radius (R obs ).

Table 8 .
Estimated values of grid size ∆x M based on the lattice-based surrogate model simulation results for the multisize case using the effective horizontal diffusivity (κ H ) and the effective vertical diffusivity (κ V ).The resulting transition (leave) probabilities in the square lattice surrogate model are also shown.