A Lattice-Boltzmann scheme for the simulation of diffusion in intracellular crowded systems

The intracellular environment is a complex and crowded medium where the diffusion of proteins, metabolites and other molecules can be decreased. One of the most popular methodologies for the simulation of diffusion in crowding systems is the Monte Carlo algorithm (MC) which tracks the movement of each particle. This can, however, be computationally expensive for a system comprising a large number of molecules. On the other hand, the Lattice Boltzmann Method (LBM) tracks the movement of collections of molecules, which represents significant savings in computational time. Nevertheless in the classical manifestation of such scheme the crowding conditions are neglected. In this paper we use Scaled Particle Theory (SPT) to approximate the probability to find free space for the displacement of hard-disk molecules and in this way to incorporate the crowding effect to the LBM. This new methodology which couples SPT and LBM is validated using a kinetic Monte Carlo (kMC) algorithm, which is used here as our "computational experiment". The results indicate that LBM over-predicts the diffusion in 2D crowded systems, while the proposed coupled SPT-LBM predicts the same behaviour as the kinetic Monte Carlo (kMC) algorithm but with a significantly reduced computational effort. Despite the fact that small deviations between the two methods were observed, in part due to the mesoscopic and microscopic nature of each method, respectively, the agreement was satisfactory both from a qualitative and a quantitative point of view. A crowding-adaptation to LBM has been developed using SPT, allowing fast simulations of diffusion-systems of different size hard-disk molecules in two-dimensional space. This methodology takes into account crowding conditions; not only the space fraction occupied by the crowder molecules but also the influence of the size of the crowder which can affect the displacement of molecules across the lattice system.


Background
Several microorganisms are used in the conversion of raw materials to added-value products, for example Actinobacillus succinogenes has been used for the synthesis of succinic acid from crude bio refinery glycerol [1,2]. The analysis and simulation of the factors affecting the metabolism of these organisms allow the further identification of the strategies needed for its manipulation in order to increase the formation of the metabolite of interest over other by-products.
As it is known the environmental conditions and the properties of the medium play an important role in the metabolism. The intracellular processes are carried out in a complex, heterogeneous, and crowded medium composed by solid components (macromolecules, enzymes, etc.) in a fluid phase called cytoplasm (in 3D) or in cell membranes (2D) [3,4], where for prokaryotes the diffusion is the primary mean of intracellular motion.
According to a drawing proposed by Goodsell [5], if the cytoplasm of Escherichia coli is divided into 600 cubes of (100 nm) 3 , an average of 130 glycolytic enzymes and 100 from the Krebbs cycle are present in each cube in addition to the metabolites and other compounds, which all together comprise a very large number of molecules for the simulation of the intracellular environment. Henceforth, we will use the terms molecule and particle interchangeably to refer to the intracellular macromolecules, e.g., proteins.
The solid components of the cell occupy about 40 % of the total volume [6] and 25 % of the area of a typical membrane [7]. Due to the impossibility that two molecules occupy the same space at the same time (steric effects), these crowding conditions affect the intracellular process in two opposite ways: 1) decreasing the diffusion of the molecules [8], and 2) increasing the thermodynamic activity of the metabolites [6], and therefore enhancing the reaction rate, and modifying the thermodynamic feasibility of some reactions. A review of the crowding effects can be found in [9].
In particular the study and simulation of the diffusion process can be carried out using several methodologies at different levels of detail. One of the most popular is the Monte Carlo (MC) algorithm [10][11][12][13][14][15][16][17] (microscopic level) where each molecule is tracked during its journey through the cell so the restriction of the impenetrability of the molecules is satisfied in a straightforward way.
MC is a powerful technique and easy to implement, however it is limited to short simulation times, restricted lattice/domain sizes, and/or reduced number of molecules because of the large computational costs. Besides, due to the stochastic nature of MC, it requires several simulations to smoothen the noise of the results by computing average quantities. Moreover, in most cases the molecules are considered to be of the same size, so the size effect could be hidden [10,11].
On the other hand, Lattice Boltzmann (LBM) [18] is a mesoscopic method which allows the efficient simulation of the dynamics of collections across a defined lattice according to expressions that conserve mass and momentum [19]. Here, the solute transport is simulated either 1) treating the solute as another fluid and solving a multicomponent problem (active solute component) or 2) assuming that the molecules do not have velocity of their own so they are advected by the fluid (passive solute component). In both cases the volume of the solute is neglected. See Sukop and Thorne [20] for a review.
Alternatively, in particle suspensions simulations the motion of each molecule is described by a hard sphere model (with the drawback of being computationally expensive for large numbers of particles) while the fluid is described by LBM [21,22]. This is similar to other hybrid methods used for example to follow the enzymes' motion with MC, or the tumour growth with Cellular Automata [23], while the passive transport of the metabolites and the fluid is simulated by LBM.
Since LBM computes the evolution of the average molecules' density, it represents a good alternative to simulate the diffusion of large number of intracellular macromolecules or even metabolites. However, since classical LBM does not take into account the volume of the molecules, and therefore the effect of obstacles on the molecules' diffusion, it may overestimate the degree of mixing of the system analysed.
The displacement of a molecule depends on the probability P to find enough empty space to move at every step of its journey. Scaled Particle Theory (SPT) is a method that allows the estimation of this probability P which is a function of the radii and concentration/number of molecules present in the system. SPT also has been used to investigate the effect of macromolecular crowding on solvation [24], thermodynamic activity of proteins [25] and of metabolites [26].
The aim of this paper is to incorporate the crowding effect on the LBM simulation of the particles' diffusion. For this, a methodology is proposed for coupling SPT and LBM, allowing in this way faster simulations for systems with a large number of molecules of different size under crowded conditions, such as the intracellular environment. Here, we consider the diffusing molecules as passive solute components assuming that the fluid phase is at rest within the cell. In particular, this paper focuses on 2D simulations of macromolecules' diffusion which are relevant for the study and analysis of lateral diffusion of proteins in membranes [27][28][29]. For validation purposes the results are compared with those obtained from kinetic Monte Carlo (kMC) algorithm [14,15].

Methods
In the classical LBM [18], the system is represented by a regular lattice, where the molecules located at the same site or node at time t may interact with each other (collision step), and then according to a set of rules, some particles move to one of their neighbouring lattice sites (known as the streaming step), where they will interact with molecules from other nodes at time t + Δt and so on.
The methodology we propose here corrects the average number of molecules that enter a neighbouring lattice site, taking into account crowding effects and it can be summarised as follows: (1) Solve the classical LBM to find the number of molecules that will try to move into the d direction at time t + Δt (F d LB ). (2) Use F d LB to estimate the corrected number of molecules that actually enter the target site, F d , by solving the explicit formulation (see below) constrained by the size and number of the molecules, as well as the size of voxels or sites in which the lattice is divided.
(3) Use the F d values obtained in (2) for the streaming step (in the same way as in the classical LBM), and go back to point (1).
Here we use the term crowding-Lattice Boltzmann Method (cLBM) to distinguish this proposed methodology from the classical LBM, which in principle considers volumeless molecules.

Lattice Boltzmann Method
For comparison purposes with the kMC algorithm, in this paper the D2Q5 scheme [23] (Fig. 1b) is implemented, consisting of 5 possible directions in which the molecules can move, in a 2D lattice. The lattice is divided in squares of Δx [nm] side, called "voxels", whose position is identified by the index (i,j) (Fig. 1a).
The evolution of the distribution function of the species sp, F d,sp , in the diffusion system is given by the discrete Boltzmann equation [19]. We use the BGK approximation to estimate the nonreactive collision term Ω d,sp diff , which is given by [19].
where the weight factor w d is 0 for d = 0, and 1/4 for d = 1,2,3,4, while m sp is the mass of one single molecule of type sp. The macroscopic density of species sp, ρ sp , is expressed as: The expression for the relaxation parameter ω sp (indicated in Eq. (2)) can be deduced using the Enskog-Chapman procedure, and is given by [30] Due to the fact that the BGK model was formulated for non-crowded systems [19] D sp 0 takes the value of the diffusion coefficient for diluted solutions. Here, D sp 0 is considered independent of the position (i,j).

Crowding-adaptation of LBM
According to Fick's first law the diffusive flux of volumeless molecules (J) from one region to another is proportional to the gradient of the concentration. However, when the volume of the molecules is important and/or the solution is not considered diluted, J should be proportional to gradient of the activities [31].
Since the molecules in a defined system occupy a volume in space, not all the system's volume is available to the centre of mass of a test molecule [6].
The activity a is a term that describes the number of molecules per available volume, while the concentration C is the number of molecules per total volume [6]. Both variables are related by Eq. (6), where the activity coefficient γ indicates deviations from an ideal solution.
Since a sp is a dimensionless variable, for the purpose of dimensional consistency in this paper the standard concentration of the species sp (C sp st ) is considered equal to 1 molecule per voxel, i.e., diluted solutions where there is no crowding influence. Hence, the a sp value is not affected by C sp st , therefore Eq. (6) is simplified to a sp = γ sp C sp . γ is defined as the ratio of the total volume (or area in 2D systems) divided by the available volume (or area in 2D), in other words, it is inversely proportional to the probability to find free space in the system.
From the above and for the lattice-system analysed here, the flux of the molecules (J sp ) of the species sp in the x direction from voxel (i,j) to (i,j + 1) is defined in Eq. (7). The reverse flux is J sp(i,j + 1 → i,j) = − J sp(i,j → i,j + 1) .
where the concentration involved in the term a sp(i,j → i,j + 1) is proportional to the number of molecules (per voxel) trying to move from (i,j) to (i,j + 1) at time t, i.e., F 1,sp (2), and using as concentration the values F 1,sp LB (i, j, t) and F 3,sp LB (i, j + 1, t) yields: The use of a constant D sp independent of the concentration could be questionable for crowded and inhomogeneous systems since the presence of background/crowder molecules hinders the movement of a test molecule.
Due to the steric effects, the presence of background molecules leads to a reduction of the available space where the test molecule can move due to Brownian motion. Therefore the probability to find free space next to the test molecule also decreases.
Considering that the rate of Brownian displacement (or diffusion) is a function of the work required to free the target space from background molecules (ΔW), Muramatsu and Minton [32] proposed the relation where D sp 0 is the diffusion coefficient of species sp in diluted solutions. ΔW depends on the size and shape of the space required. The probability of observing the spontaneous formation of such free space as result of a fluctuation is [33,34] where k B is the Boltzmann constant, and T is the temperature. Assuming a well-mixed system (or subsystem), the probability P sp is defined as the available volume divided by the total volume, i.e., the inverse of the activity coefficient Because the diffusion coefficient D sp is a function of the available volume in the target voxel where the test molecule is trying to move in (and therefore dependent of the position (i,j)), D sp cannot be factored as it is in Eq. (8), so Substituting Eq. (9) and (10) in Eq. (12), the flux J sp(i,j → i,j + 1) can be rewritten as: where P sp is the probability of species sp to find available space in the target voxel. Eq. (13) looks like the Teorell formula for non-perfect systems [31]. In order to conserve the driving force, the corrected values F 1,sp (i,j,t) and F 3,sp (i,j + 1,t) should also be related to the flux J sp(i,j → i,j + 1) , so that Comparing Eq. (13) and Eq. (14), the new F 1,sp (i,j,t) and F 3,sp (i,j + 1,t) values are found as The molecules that could not move into their corresponding target voxel (i next ,j next ), i.e., the voxel next to (i,j) in the direction d, will remain in the current (i,j). Therefore in order to conserve mass F 0,sp becomes The activity coefficient γ sp (i,j,t) is a function of the size and shape of the molecules present in the voxel/site (i,j,t).
In this paper, the Scale Particle Theory [33,34] is used to approximate γ sp in a mixture of (non-overlapping) hard disk-molecules (Eq. (20) and (21)) which could be of different radii r [31]. For this we assume that the volume and temperature of the subsystems/voxels are held constant. In order to simplify the notation, the position-time dependence (i,j,t) of the variables γ sp and S x has been dropped from Eq. (20) and (21).
where S x (0 ≤ x ≤ 2) is given by Note that if all the particles simulated are considered to be point-like molecules (r i → 0, where the space fraction occupied by the molecules tends to zero S 2 → 0), and/or they are in a non-crowded system ( (1) in Eq. (20) so that the probability P sp (Eq. (11)) would be equal to 1, therefore F d,sp becomes equal to F d,sp LB (Eq. (17)). The use of P sp (or its equivalent γ sp − 1 ) in Eq. (17) restricts the maximum number of molecules in a voxel, and also avoids the exchange of positions between molecules. When two molecules moving in opposite directions are adjacent and have volume (or area in 2D systems), one acts as an obstacle for the other, so exchanging positions between them should be prohibited. In other words the displacement of a molecule is limited by the probability to find empty space P sp .

Streaming step
The new F d ,sp values are used in the streaming step in the same way as in LBM (Eq. (22) − (26) below). Then the scheme goes back to step 1 and the procedure is repeated until the simulation end time is reached.
Unlike (microscopic-) MC methods that track the movements of every single particle, (mesoscopic-) cLBM simulates the movement of collections of molecules across a (lattice) system. Hence, when a large number of molecules is considered and/or for long time simulations, cLBM simulations are computationally more efficient than MC methods. However, since each (cLBM-) voxel can fit more than one molecules and assuming that the voxel is well mixed, which is not necessarily true especially for large voxels, some discrepancies could be found between the diffusion results computed by cLBM and MC (in this paper we use kMC).
In order to validate the crowding-adaptation of LBM presented above, the simulation of a diffusion system was carried out by an on-lattice kMC algorithm [14,15], whose results are used as our computational experiments. A brief description of the kMC algorithm is given below.

Lattice Kinetic Monte Carlo (kMC) algorithm
In the on-lattice kMC [14] algorithm each site of the lattice can be occupied by at most one molecule. Each molecule can move to the one of the 4 neighbouring sites (top, bottom, right, or left), as long as they are free, i.e., there is no other molecule in the target site.
The basic idea of the procedure is the following: At every time step 1) Identify the classes of species, i.e., the combinations of adjacent species that can react or diffuse (if there is an empty site next to a molecule). 2) The rates (including the diffusion rate) of all the possible events or processes are listed and their cumulative value is computed. 3) An event is probabilistically chosen, e.g., the next molecule to move and the site where it will take place, using a random number and the cumulative value of the rates listed in point 2. 4) The diffusion (or reaction) of the chosen molecule(s) takes place. 5) A variable time step value is also probabilistically estimated using another random number and the cumulative value of the rates. 6) The number of species and classes in a region and/ or in the whole lattice is updated as well as the time. 7) Go back to point 1 until the end simulation time is reached.
See [14,15] for detailed information of the algorithm and the corresponding equations for each step of the process.
MC algorithms have been widely used for the simulation of processes in crowded media [10-13, 16, 17]. The main difference between kMC and MC is that the former use a constant Δx and a variable Δt, while in MC both parameters are constant.
However for both cases, kMC and MC, a molecule only can move if the target (adjacent) site, which was randomly chosen is free/empty. Hence, the excluded volume restriction is straightforward satisfied in both kMC and MC. A comparison of the results of diffusion simulation obtained from these two algorithms is shown in the Additional file 1.

Results and discussion
For the diffusion examples presented in this section the following assumptions have been made: 1) The system analysed consists of a square lattice of (1000 nm) 2 , divided in equal-sized voxels or sites of Δx per side whose value is indicated in each example. 2) Each voxel is well-mixed.
3) The fluid phase is considered as a continuum and it is at rest. 4) The system has periodic boundaries. 5) The mass of all species are assumed equal to 1 g per molecule, therefore the value of ρ sp is equal to the number of molecules, i.e., The methodologies LBM and cLBM were programed in MATLAB R2011a (The MathWorks, Natick, MA), while the lattice kMC was implemented in Fortran 90.

Validation of cLBM: a lattice-model
For validation purposes, we consider a lattice model (Fig. 2) for the diffusion of molecules here represented by non-rotating squares of (10 nm) 2 which have a square uniform packing order. Taking into account that the system' size is (1000 nm) 2 the maximum number of molecules that can fit inside is 10,000. This type of system is consistent with the ones commonly used in on-lattice kMC simulations where a molecule can move to one of its neighbouring sites at every time step.
According to the Δx value used in the cLBM simulations, for comparison purposes, the average results obtained (after 1000 repetitions) at every time step Δt by kMC were coarse grained in lattice regions equivalent to voxels of side length equal to Δx.
Due to the square uniform packing order of the same size molecules in this lattice model the activity coefficient γ sp (which is required in cLBM) given by Eq. (20) and (21) can be simplified to where A i is the area of a molecule of species i. The diffusion of two types of molecules has been simulated using the parameters shown in Example 1 of Table 1. For this, the lattice is divided into 10 vertical regions, where the voxels (of Δx = 100 nm) of the first column of the lattice have been filled with 100 molecules of type A and B as indicated in the Fig. 3a.
Once the molecules begin to diffuse, it is possible to compare their movement summing the number of molecules A (used as tracer molecule) in each vertical region of the lattice at different times. Figure 3b shows that LBM and cLBM predict the same diffusion profile, which in turn is close to the results estimated by kMC. This is because there are no obstacles in the horizontal direction and therefore cLBM results are the same as the ones from If we compare the molecules' movement in the vertical direction (summing the number of molecules A but now in each horizontal region), then it is possible to see differences between the diffusion profile of cLBM and LBM (Fig. 3c). Results show that the system evolves faster, i.e., LMB computes that more molecules move or diffuse to their neighbouring voxels, than cLBM. Hence, the profile predicted by cLBM in the vertical direction (where the obstacles are initially located) is closer to the one obtained by kMC than that predicted by LBM.
The relative error of the results obtained by cLBM and/ or LBM compared with kMC is calculated as the Frobenius norm of the difference of the matrixes containing the number of molecules sp in each voxel predicted by kMC (ρ kMC ) and cLBM (ρ cLBM ), or by kMC and LBM (ρ LBM ), divided by the total number of molecules sp simulated, i.e., Since m sp of all sp is assumed to be equal to 1 g per molecule, we keep the term ρ sp to represent the number of molecules in the system. A comparison of the error(ρ kMC − ρ cLBM ) and error(ρ kMC − ρ LBM ) for the diffusion of species A in the Example 1 indicates that cLBM predictions are more accurate than those of LBM (Fig. 3d). For example, the relative error estimated by cLBM at time 30 ms is error A = 0.93 % (which indicates that from 100 molecules A simulated, cLBM predicts that 0.93 molecules are allocated in different voxel's positions than those predicted by kMC, in other words there are deviations in the distribution of 0.93 molecules A from the 100 simulated in this example), while the error from LBM amounts to error A = 16.39 %. This is because LBM considers point-like molecules so that they will always find enough space to fit in the target voxel, causing LBM to over-predict diffusion in the vertical direction.
Nevertheless, if we assume that both molecules A and B are of the same type then the relative error between kMC and LBM decreases, becoming equal to that between kMC and cLBM (Fig. 3e). This means that LBM is a good alternative for quick simulations of one type of molecules, but when two or more species are present then it can over predict the system's degree of mixing.
As was pointed out by Li et al. [35], the relaxation parameter value affects the accuracy of LBM. According to Eq. (5) and the chosen parameters Δx and Δt in Example 1, both ω A and ω B were estimated equal to one. In order to show the dependency between the error A and ω A , we tested different Δt values, maintaining constant Δx = 50 nm, which corresponds to ω A ranging from 0.5 to 1.  The order of magnitude of diffusion coefficients used agrees with that reported by Elowitz et al. [41] The results indicate that ω A = 1 gives the lowest error A (Fig. 4a), therefore in the following we use Δx and Δt values that allow setting ω sp close to one. From Eq. (5) when ω sp = 1, then the diffusion coefficient , is equal to the stability limit, Δt≤ Δx 2 4D sp , given by the finite difference (FD) approximation of the diffusion equation ∂ρ sp ∂t ¼ D sp ∇ 2 ρ sp . Under this condition (ω sp = 1), LBM estimates F LB d;sp i next ; j next ; t þ Δt (1)-(3)), that substituting in Eq. (4) gives the same equation obtained by FD at the stability limit [19], i.e.,  Table 1 Moreover, since cLBM is formulated as an explicit method that requires information of the neighbouring voxels at time t to estimate ρ sp of a voxel (i,j) for the next time t + Δt, decreasing Δt also diminishes the error predicted between kMC and cLBM for Example 1 as shown in Fig. 4b, where 3 different sets of coefficient Δt and the corresponding Δx that keeps constant ω sp = 1 were tested.
In all cases, the LB methods required s time steps of length Δt, before reaching a quasi-stable error. Therefore the smaller the chosen Δt (and the corresponding Δx) the faster the system reached that state. Since good results were obtained with Δt = 1.25 ms and Δx = 50 nm (where a voxel allows to fit a maximum of 25 molecules) compared with the more accurate but slower Δt = 0.2 ms and Δx = 20 nm, in the following we use Δx = 50 nm.
In the previous example the crowding conditions are in the vertical direction so that reduced displacement of species A is in this direction (Fig. 3c). Continuing the analysis of the crowding effect on the molecules' displacement, a second example is proposed for the diffusion of three types of molecules initially allocated in the region indicated in Fig. 5a (the parameters are shown in Example 2 of Table 1).
For the purpose of this example, a region consists of two columns of 20 voxels with Δx = 50 nm, where each voxel of region 3 (Fig. 5a) is filled with 25 molecules of species A. Region 4 and 5 are filled with equal numbers of species B and C, respectively. The results of the diffusion of species A show that cLBM predicts the same behaviour as kMC as can be seen in Fig. 5b.
On the other hand, LBM predicts a symmetric movement of A despite the fact that species B acts as an obstacle in its way (Fig. 5b). The same symmetric profile is obtained for species B and C ( Fig. 5c and d).
A comparison of the relative error between kMC − cLBM and kMC − LBM is shown in Fig. 6. As was also observed in Example 1, the error computed with the cLBM results is smaller than that computed through the LBM results. In fact, at time 30 ms the errors given by cLBM (error A = 1.103 %, error B = 0.833 %, error C = 1.172 %) represent deviations in no more than 2 % of the total number of molecules simulated for each species.
Additionally, a comparison of the mean squared displacement (MSD) of a tracer molecule, i.e., the displacement from one voxel to another, was carried out under different crowding conditions (the parameters are given in Example 3 of Table 1). The total number of tracer molecules represents 1 % of the lattice area of the 2D system simulated by cLBM and LBM.
The results show ( Fig. 7a and b) that only cLBM is sensitive to increments in the concentration of the background/crowder molecules, which is reflected in a reduction of the displacement of the tracer particle.
The comparison of the MSD computed by cLBM (lines and dashed lines in Fig. 7a) and kMC (circles in Fig. 7a) reveals a very good agreement between both methods. The error estimated between MSD cLBM and ρ sp i next ; j next ; t þ Δt À Á ¼ ρ sp i next −1; j next ; t À Á þ ρ sp i next þ 1; j next ; t À Á þ ρ sp i next ; j next −1; t À Á þ ρ sp i next ; j next þ 1; t À Á 4 ð30Þ Fig. 4 Relative error of the molecules distribution estimated by cLBM − kMC for Example 1. a Different relaxation parameter ω A were tested but maintaining constant Δx = 50 nm. b Different Δx and Δt values were tested maintaining ω A = 1 100, was found to be lower than 0.7 % for all the crowding conditions tested (Fig. 7c). These small differences in the tracer's displacement may be due to the fact that LB methods are unable to identify and quantify the motion of the molecules that remain in the same voxel at time t + Δt, i.e., LB methods only quantify the "effective" displacement made when the molecules pass from one to another voxel, but not how many movements have to be executed before entering to the next voxel as is the case with lattice kMC.
Despite the fact that kMC and cLBM were implemented in different computing languages (Fortran and Matlab, respectively, using an Intel Xeon 5160, CPU 3.00 Hz processor) a comparison of the execution time   Certainly, the speed of cLBM-simulations depends on the time simulated and the chosen Δt, but also other factors can affect its execution time, e.g., the number of different species analysed and the number of voxels in which the system is divided (i.e., the chosen Δx). This is because F d,sp (i,j,t) is estimated for each species sp at every voxel position.

Different size molecules in cLBM
Up to this point we have analysed a lattice-model were only molecules of the same size can be considered, however as was pointed out by Vilaseca et al. [10,11], the size (and also the shape) of the molecules could be an important parameter in their movement or diffusion.
In order to study the influence of the size of the molecules on the diffusion process we simulate the motion of a tracer molecule in a 2D crowding system composed by 5 types of background particles (all of them of circular shape with different radii and concentrations) which are randomly located and all together occupy 30 % of the total lattice space (the corresponding parameters are given in Example 4 of Table 1).
According to the SPT assumption of a well-mixed voxel, here the molecules can be anywhere unlike the previous lattice model where only a square uniform packing order is allowed. Figure 8a shows the displacement across the lattice predicted by cLBM for a tracer molecule of different radius magnitude whose total concentration (i.e., the total number of molecules in the lattice) is 796 molecules. As it can be seen (Fig. 8) more movements from one voxel to another were detected when a small radius is assumed. This is because the probability of the tracer species to move to the target voxel increases inasmuch as the second and third term of Eq. (20) disappear when r tracer = 0 nm. In other words, the available space for point-like molecules equals the free space in the target voxel, i.e., the total area not occupied by other molecules represented by the term 1 − S 2 in Eq. (20).  Table 1 Moreover, the comparison of the displacement of the tracer molecule with r tracer = 0 nm obtained by cLBM (broken pink line in Fig. 8a) and that estimated by the classical LBM for the same Example 4 (broken pink line in Fig. 8b) indicates that even if the metabolites are considered as point-like particles they are still affected by the molecular crowding conditions, unless all the molecules simulated in the system are also point-like molecules.
Finally, the effect of the size of the crowder molecules on diffusion was investigated. For this, 3 different diffusion simulations of a system composed of tracer molecules of r tracer = 1.5 nm which occupies 1 % of the lattice area, and crowder molecules having different radii: 2, 1.5, and 1 nm, representing 30 % of the total lattice space were performed.
Since the area covered by the crowder is kept constant despite the fact that the radius is modified, more particles of r crowder = 1 nm are simulated compared with those used if r crowder is 1.5 or 2 nm. A comparison of the mean squared displacement of a tracer molecule (Fig. 9) under such conditions indicates that the particles' motion decreases when the number of crowder molecules increases.
This suggests that not only the space fraction occupied by crowders is important but also the way it is covered,  Table 1 i.e., the size and number of molecules. Similar behaviour was found by Vilaseca et al. [10] for the diffusion simulation of molecules having different sizes using a Monte Carlo algorithm.
The (2D-) cLBM methodology presented in this paper is useful to simulate the diffusion in systems such as cell membranes. Furthermore, cLBM can be easily extended for diffusion simulations in three dimensions. Following the same reasoning for the estimation of the corrected value F d,sp (Eq. (17)) there are only small differences between the equations used for 3D and those for 2D. For example the variables F d,sp LB , F d,sp , ρ sp , Ω d,sp diff , and P sp will depend on the position (i,j,k), where the index k denotes the z-coordinate of the space. Besides, in the calculation of ω sp (Eq. (5)) the coefficient 4 used in the second term of the denominator would change to the value 6, see [30] for the deduction of the relaxation parameter equation in 3D. Finally, the 2D SPT equation used to estimate γ sp (Eq. 20) would also change to consider spherical molecules in 3D systems, see [33,34] for more details.
Therefore, it is expected that 3D-cLBM will maintain the computational efficiency of 2D-cLBM. However, since the number of voxels will increase due to the dimensionality of the system, the time required to simulate the motion of molecules will increase accordingly.
On the other hand, reactions have been simulated by the classical LBM [35][36][37][38] through the law of mass action, which assumes that the reaction rate is proportional to the concentration of the reactants. However, more general speaking the reaction rate is proportional to the activity of the reactants (Eq. 6).
Considering the above, cLBM can be also extended to simulate reaction-diffusion systems in crowded media. For this, the reaction process takes place at every time step and in every voxel. The evolution (consumption/ production) of the species in a well-mixed voxel is a function of the activity of the reactants in the voxel, i.e., the species' concentration (or number of molecules per voxel equivalent to ρ sp ) and γ sp estimated by SPT.
Since the reaction processes have to be evaluated at every time step and in every voxel, an increase in the execution time is expected for the simulation reactiondiffusion systems proportional to the parameters Δx and Δt. Therefore strategies and simplifications will be required to optimise the simulation run time.
When the difference between the diffusion and reaction rates is important, strategies like the time splitting method [30] can be adopted into LBM and cLBM methods. For example, in the scheme DRD splitting the time step Δt is divided in two Δt = Δt 1 + Δt 2 . First, the diffusion is carried out during Δt 1 , then the reactions take place using the number of molecules present in the voxel after time Δt 1 . Finally, the updated number of molecules (counted after the reactions took place) are allowed to diffuse during time Δt 2 .
This paper focuses on the diffusion of macromolecules. However, when the simulation of a mixture of metabolites and macromolecules (where the diffusion coefficients could be of different order of magnitude) is required, then some strategies are needed to balance the computational cost with accuracy of cLBM. For example, the use of a very small Δx will increase the number of voxels and therefore the computation effort required for the simulation. The use of the time splitting method (or similar strategies) will help in the simulation of diffusion of a mixture of metabolites and macromolecules. In this way cLBM may be extended for the reaction-diffusion simulation in more realistic (and crowded) scenarios.

Conclusions
The LBM predicts with great accuracy the diffusion of particles under ideal conditions, i.e., considering point-like molecules, and/or non-crowding systems, and/or when only one type of molecules is simulated. However, if conditions change, for example, a system involving more than two species at crowding conditions, LBM predictions increasingly deviate from our on-lattice kMC-based computational experiment.
Although small discrepancies were found between the cLBM and kMC results (differences that were expected due to the level of detail inherent in each method), the proposed crowding adaptation of LBM is able to predict the same behaviour in the species diffusion profile. This suggests that the coupled SPT-LBM can be considered as a computational alternative for fast simulations of diffusion systems with a large number of molecules of different size and/or for long times, under crowding conditions at a fraction of the computational cost compared to a molecular (microscopic) method such as kMC.
Nevertheless, as in other mesoscopic methods, the saving in the execution time is accompanied by a reduction in the information that cLBM can provide compared with that obtained from a microscopic method. For example, the quantification of the total displacement of the molecules on time for the parameters estimation of anomalous diffusion [10,11].
The accuracy of cLBM is influenced by the chosen voxel size Δx and time increment Δt, both related by the relaxation parameter ω sp . It was found that the use of ω sp values close to one gives better results than any other in the range 0 ≤ ω sp ≤ 1. Moreover, cLBM being an explicit method, small values of Δt (maintaining ω sp = 1, which involves reducing Δx) also reduce the error between the proposed methodology (cLBM) and kMC.
Other factors that can affect the accuracy of cLBM is the area (or volume in 3D) fraction occupied by molecules, and the presence of immobile species. In particular, SPT (used in cLBM to estimate the probability to find available space P sp ) works well for low to moderate area (or volume in 3D) fraction occupied by molecules [39,40]. At high area (or volume) fraction occupied by molecules or when immobile particles are considered, the spatial distribution of the molecules inside a voxel could form free space "pockets" which will be not available for the incoming molecules from neighbouring voxels but that SPT takes into account in the calculation of P sp . This would lead to an overestimation of the molecules' diffusion.
Regarding the influence of the size of the particles on the diffusion process, a reduction in the mean squared displacement of a tracer molecule when its size is increased was observed, as well as when the size of the crowders is decreased (but maintaining constant the lattice fraction occupied by them). Hence, the incorporation of small molecules, e.g., metabolites, in the simulation system can affect the diffusion profile predicted for macromolecules.
Even though cLBM requires the species' concentration of the neighbour voxels at time t to compute the results at t + Δt, therefore the LBM's local feature is lost, the correction for the crowding effects is external to the estimation of LBM distributions, i.e., F d,sp LB , so that alternative LBM schemes can potentially be implemented within cLBM, e.g., for the simulation of reaction-diffusion systems.