 Research article
 Open Access
 Published:
A LatticeBoltzmann scheme for the simulation of diffusion in intracellular crowded systems
BMC Bioinformatics volume 16, Article number: 353 (2015)
Abstract
Background
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.
Methods
In this paper we use Scaled Particle Theory (SPT) to approximate the probability to find free space for the displacement of harddisk 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".
Results
The results indicate that LBM overpredicts the diffusion in 2D crowded systems, while the proposed coupled SPTLBM 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.
Conclusions
A crowdingadaptation to LBM has been developed using SPT, allowing fast simulations of diffusionsystems of different size harddisk molecules in twodimensional 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 addedvalue 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 byproducts.
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–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–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 ^{LB}_{ d } ).

(2)
Use F ^{LB}_{ d } 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 crowdingLattice Boltzmann Method (cLBM) to distinguish this proposed methodology from the classical LBM, which in principle considers volumeless molecules.
CrowdingLattice Boltzmann Method (cLBM)
1. 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].
The superscript LB is used to distinguish F ^{LB}_{ d,sp } that is calculated from the classical LBM from the crowdingcorrected value F _{ d,sp }. Both F ^{LB}_{ d,sp } and F _{ d,sp } represent the concentration or number of molecules of the species sp in a voxel that try to move to a neighbouring site, so they are given in [molecules per voxel]. Other concentration units e.g., [mol per voxel] can be used, but dimensional changes in other variables are required for consistency.
We use the BGK approximation to estimate the nonreactive collision term Ω ^{diff}_{ d,sp } , which is given by [19].
Assuming that the fluid phase is at rest, the equilibrium distribution function F ^{eq}_{ d,sp } takes the form [23]
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 EnskogChapman procedure, and is given by [30]
Due to the fact that the BGK model was formulated for noncrowded systems [19] D ^{0}_{ sp } takes the value of the diffusion coefficient for diluted solutions. Here, D ^{0}_{ sp } is considered independent of the position (i,j).
2. Crowdingadaptation 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 ^{st}_{ sp } ) 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 ^{st}_{ sp } , 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 latticesystem 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 ^{LB}_{1,sp } (i, j, t), while a _{ sp(i,j + 1 → i,j)} is proportional to F ^{LB}_{3,sp } (i, j + 1, t).
Substituting Eq. (1) in Eq. (2), and using as concentration the values F ^{LB}_{1,sp } (i, j, t) and F ^{LB}_{3,sp } (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 ^{0}_{ sp } 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 wellmixed 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 nonperfect 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
Generalizing, the corrected values F _{ d,sp } (d = 1,2,3,4) are calculated as
Or their equivalent:
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 (nonoverlapping) hard diskmolecules (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 positiontime 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 pointlike molecules (r _{ i } → 0, where the space fraction occupied by the molecules tends to zero S _{2} → 0), and/or they are in a noncrowded system (\( {\displaystyle \sum_{i=1}^l{\rho}_i}\to 0 \), hence S _{ x } → 0), then ln(γ _{ sp }) → − ln(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 ^{LB}_{ d,sp } (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 }.
3. 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 crowdingadaptation of LBM presented above, the simulation of a diffusion system was carried out by an onlattice 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 onlattice 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 equalsized voxels or sites of ∆x per side whose value is indicated in each example.

2)
Each voxel is wellmixed.

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., \( {\displaystyle \sum_d{F}_{d,sp}} \). Notice that m _{ sp } can take other values if it is required.
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 latticemodel
For validation purposes, we consider a lattice model (Fig. 2) for the diffusion of molecules here represented by nonrotating 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 onlattice 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 LBM. The same happens when the volume (or area in 2D) of the molecules are zero (data not shown).
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 pointlike molecules so that they will always find enough space to fit in the target voxel, causing LBM to overpredict 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 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 \( {D}_{sp}=\frac{\Delta {x}^2}{4\Delta t} \), is equal to the stability limit, \( \Delta t\le \frac{\Delta {x}^2}{4{D}_{sp}} \), given by the finite difference (FD) approximation of the diffusion equation \( \frac{\partial {\rho}_{sp}}{\partial t}={D}_{sp}{\nabla}^2{\rho}_{sp} \). Under this condition (ω _{ sp } = 1), LBM estimates \( {F}_{d,sp}^{LB}\;\left({i}_{next},{j}_{next},t+\Delta t\right)=\frac{\rho_{sp}\left(i,j,t\right)}{m_{sp}}{w}_d \) (see Eq. (1)–(3)), that substituting in Eq. (4) gives the same equation obtained by FD at the stability limit [19], i.e.,
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 quasistable 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 MSD _{ kMC }, \( erro{r}_{MSD}(t)=\frac{norm\left(MS{D}_{kMC}(t)MS{D}_{cLBM}(t)\right)}{MS{D}_{kMC}(t)}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 (CPU time) for the diffusion Example 1 (t _{ kMC } = 315 s per run simulation or repetition vs t _{ CLBM } = 7.13 s), and Example 2 (t _{ kMC } = 441 s per repetition vs t _{ CLBM } = 18.31 s) reveals the potential use of cLBM for faster simulations of larger systems. The total time required by kMC depends on the number of repetitions performed. The minimum error estimated between kMC and cLBM (Eq. (29)) for Example 1 after different kMC simulation repetitions was 1.43 % (50 repetitions ~ t _{ kMC } = 4.37 hr), 1.32 % (100 repetitions ~ t _{ kMC } = 8.75 hr), 0.71 % (500 repetitions ~ t _{ kMC } = 43.75 hr), and 0.51 % (1000 repetitions ~ t _{ kMC } = 87.5 hr). While in Example 2 the minimum error computed was 0.89 % (50 repetitions ~ t _{ kMC } = 6.12 hr), 0.871 % (100 repetitions ~ t _{ kMC } = 12.25 hr), 0.61 % (500 repetitions ~ t _{ kMC } = 61.25 hr), and 0.54 % (1000 repetitions ~ t _{ kMC } = 122.5 hr).
Certainly, the speed of cLBMsimulations 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 latticemodel 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 wellmixed 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 pointlike 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).
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 pointlike particles they are still affected by the molecular crowding conditions, unless all the molecules simulated in the system are also pointlike 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, 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 ^{ LB}_{ d,sp } , F _{ d,sp }, ρ _{ sp }, Ω ^{diff}_{ d,sp } , and P _{ sp } will depend on the position (i,j,k), where the index k denotes the zcoordinate 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 3DcLBM will maintain the computational efficiency of 2DcLBM. 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–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 wellmixed 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 reaction–diffusion 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 pointlike molecules, and/or noncrowding 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 onlattice kMCbased 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 SPTLBM 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 ^{ LB}_{ d,sp } , so that alternative LBM schemes can potentially be implemented within cLBM, e.g., for the simulation of reaction–diffusion systems.
Nomenclature
a _{ sp } Activity of the species sp [dimensionless].
A _{ sp } Area of a molecule of the species sp [nm^{2}].
C _{ sp } Concentration of the species sp [molecules nm^{−2}].
C ^{st}_{ sp } Standard concentration of the species sp [molecules nm^{−2}].
d Direction chosen by the molecules to jump to neighboring voxel [dimensionless].
D _{ sp } Diffusion coefficient [nm^{2} ms^{−1}].
D ^{0}_{ sp } Diffusion coefficient in dilute solutions [nm^{2} ms^{−1}].
error _{ sp } Relative error of molecules’ distribution predicted by kMC − cLBM or kMC − cLBM of the species sp [%].
error _{ MSD }Relative error of tracer’s MSD predicted by kMC − cLBM [%].
F _{ d,sp } Distribution function of the species sp in the direction d predicted by cLBM [molecules per voxel].
F ^{eq}_{ d,sp } Equilibrium distribution function of the species sp in the direction d [molecules per voxel].
F ^{LB}_{ d,sp } Distribution function of the species sp in the direction d predicted by LBM [molecules per voxel].
i Index that identify the position of a voxel [dimensionless].
j Index that identify the position of a voxel [dimensionless].
J _{ sp } Diffusive flux of molecules sp in 2D [molecules nm^{−1} ms^{−1}].
m _{ sp } Mass of a molecule sp [g molecule^{−1}].
next Subscript that indicates the target voxel where the molecules will move in the t + Δt [dimensionless].
P _{ sp } Probability to find available space for species sp in the target voxel [dimensionless].
r _{ sp } Radii of the species sp [nm].
k _{ B } Boltzmann constant equivalent to 1.3806 × 10^{−23} [J K^{−1}].
sp Index that identify the molecule species sp [dimensionless].
t Time [s].
T Temperature of the medium [K].
w _{ d }Weight factor for the calculation of the equilibrium function [dimensionless].
Δt Time increment [s].
∆WWork required to free the target space from background molecules [J molecule^{−1}].
ΔxSize of the voxel in which the lattice is divided [nm].
γ _{ sp }Activity coefficient of a molecule sp [dimensionless].
Ω ^{diff}_{ d,sp } Nonreactive collision term [molecules per voxel].
ω _{ sp } Relaxation parameter [dimensionless].
ρ _{ sp } Macroscopic density of species sp in a voxel [g per voxel].
ρ _{ sp } Matrix with the macroscopic density of species sp in all voxels of the lattice [g voxel].
Abbreviations
 LBM:

Lattice Boltzmann Method
 cLBM:

crowding adaptation to Lattice Boltzmann Method
 kMC:

kinetic Monte Carlo algorithm
 MC:

Monte Carlo algorithm
References
Beauprez JJ, De Mey M, Soetaert WK. Microbial succinic acid production: natural versus metabolic engineering producers. Process Biochem. 2010;45:1103–14.
Vlysidis A, Binns M, Webb C, Theodoropoulos C. Glycerol utilisation for the production of chemicals: conversion to succinic acid, a combined experimental and computational study. Biochem Eng J. 2011;58–59:1–11.
Axelrod D. Lateral motion of membrane proteins and biological function. J Membr Biol. 1983;75:1–10.
McCloskey M, Poo MM. Protein diffusion in cell membranes: some biological implications. Int Rev Cytol. 1984;87:19–81.
Goodsell DS. Inside a living cell. Trends Biochem Sci. 1991;16:203–6.
Minton AP. The influence of macromolecular crowding and macromolecular confinement on biochemical reactions in physiological media. J Biol Chem. 2001;276:10577–80.
Grasberger B, Minton AP, DeLisi C, Metzger H. Interactions between proteins localized in membranes. Proc Natl Acad Sci. 1986;83:6258–62.
Verkman AS. Solute and macromolecule diffusion in cellular aqueous compartments. Trends Biochem Sci. 2002;27:27–33.
Chebotareva NA, Kurganov BI, Livanova NB. Biochemical effects of molecular crowding. Biochem Mosc. 2004;69:1239–51.
Vilaseca E, Isvoran A, Madurga S, Pastor I, Garces JL, Mas F. New insights into diffusion in 3D crowded media by Monte Carlo simulations: effect of size, mobility and spatial distribution of obstacles. Phys Chem Chem Phys. 2011;13:7396–407.
Vilaseca E, Pastor I, Isvoran A, Madurga S, Garces JL, Mas F. Diffusion in macromolecular crowded media: Monte Carlo simulation of obstructed diffusion vs. FRAP experiments. Theor Chem Acc. 2011;128:795–805.
Berry H. Monte Carlo simulations of enzyme reactions in two dimensions: fractal kinetics and spatial segregation. Biophys J. 2002;83:1891–901.
Aranda JS, Salgado E, MuñozDiosdado A. Multifractality in intracelular enzymatic reactions. J Theor Biol. 2006;240:209–17.
Reese JS, Raimondeau S, Vlachos DG. Monte Carlo algorithms for complex surface reaction mechanisms: efficiency and accuracy. J Comput Phys. 2001;173:302–21.
Fragkopoulos IS, Theodoropoulos C. Modelling of electrochemically promoted systems. Electrochim Acta. 2014;150:232–44.
Grima R, Schnell S. A systematic investigation of the rate laws valid in intracellular environments. Biophys Chem. 2006;124:1–10.
Schnell S, Turner TE. Reaction kinetics in intracellular environments with macromolecular crowding: simulations and rate laws. Prog Biophys Mol Biol. 2004;85:235–60.
McNamara GR, Zanetti G. Use of the Boltzmann equation to simulate latticegas automata. Phys Rev Lett. 1988;61:2332–5.
WolfGladrow DA. LatticeGas Cellular Automata and Lattice Boltzmann Models – An Introduction. Berlin, Heidelberg, New York: Springer; 2005.
Sukop MC, Thorne DT. Solute transport. In: Lattice Boltzmann Modeling. An Introduction for Geoscientists and Engineers. Heidelberg, Berlin, New York: Springer; 2006. p. 117–44.
Wang L, Zhou G, Wang Q, Ge W. Direct numerical simulation of particlefluid system by combining timedriven hardsphere model and lattice Boltzmann method. Particuology. 2010;8:379–82.
Aidun CK, Clausen JR. LatticeBoltzmann method for complex flows. Annu Rev Fluid Mech. 2010;42:439–72.
Alemani D, Pappalardo F, Pennisi M, Motta S, Brusic V. Combining cellular automata and lattice Boltzmann method to model multiscale avascular tumor growth couple with nutrients diffusion and immune competition. J Immunol Meth. 2012;376:55–8.
Tang KES, Bloomfield VA. Excluded volume in solvation: sensitivity of scaledparticle theory to solvent size and density. Biophys J. 2000;79:2222–34.
Minton AP. Excluded volume as a determinant of macromolecular structure and reactivity. Biopolymers. 1981;20:2093–120.
Angeles–Martinez L, Theodoropoulos C. The influence of crowding conditions on the thermodynamic feasibility of metabolic pathways. Biophys J. 2015; doi:10.1016/j.bpj.2015.09.030.
Saxton MJ. Lateral diffusion in an archipelago. The effect of mobile obstacles. Biophys J. 1987;52:989–97.
Saxton MJ. Lateral diffusion in a mixture of mobile and immobile particles. A Monte Carlo study. Biophys J. 1990;58:1303–6.
Minton AP. Lateral diffusion of membrane proteins in proteinrich membranes. A simple hard particle model for concentration dependence of the twodimensional diffusion coefficient. Biophys J. 1989;55:805–8.
Alemani D. A Lattice Boltzmann numerical approach for modelling reaction–diffusion processes in chemically and physically heterogeneous environments. Switzerland: PhD thesis. University of Geneva; 2007.
Gorban AN, Sargsyan HP, Wahab HA. Quasichemical models of multicomponent nonlinear diffusion. Math Model Nat Phenom. 2011;6:184–262.
Muramatsu N, Minton A. Tracer diffusion of globular proteins in concentrated protein solutions. Proc Natl Acad Sci U S A. 1988;85:2984–8.
Reiss H, Frisch HL, Lebowitz JL. Statistical mechanics of rigid spheres. J Chem Phys. 1959;31:369–80.
Lebowitz JL, Helfand E, Praestgaard E. Scale particle theory of fluid mixtures. J Chem Phys. 1965;43:774–9.
Li Q, Zheng CG, Wang NC, Shi BC. LBGK simulations of Turing patterns in CIMA model. J Sci Comput. 2001;16:121–34.
Dawson SP, Chen S, Doolen GD. Lattice Boltzmann computations for reaction–diffusion equations. J Chem Phys. 1993;98:1514–23.
Chen S, Dawson SP, Doolen GD, Janecky DR, Lawniczak A. Lattice methods and their applications to reacting systems. Comput Chem Eng. 1995;19:617–46.
Theodoropoulos C, Qian YH, Kevrekidis IG. “Coarse” stability and bifurcation analysis using timesteppers: A reaction–diffusion example. PNAS. 2000;97:9840–3.
Heying M, Corti DS. Scaled particle theory revisited: new conditions and improved predictions of the properties of the hard sphere fluid. J Phys Chem. 2004;108:19756–68.
Grima R. Intrinsic biochemical noise in crowded intracellular conditions. J Chem Phys. 2010;132:185102.
Elowitz MB, Surette MG, Wolf PE, Stock JB, Leibler S. Protein mobility in the cytoplasm of Escherichia coli. J Bacteriol. 1999;181:197–203.
Acknowledgements
LAM would like to acknowledge the financial support of CONACYTMexico. LAM and CT would like to thank Prof. Alexander Gorban for useful discussions as well as his suggestions in the development of the crowding LBM methodology. The authors would also like to acknowledge the contribution of Dr. Ioannis Fragkopoulos in the adaptation of the kMC code for the biological systems in this work.
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
LAM developed the crowdingadaptation to LBM, performed the diffusionsimulations, and drafted the manuscript. CT supervised the development of the methodology, reviewed and edited the manuscript. All authors read and approved the final manuscript.
Additional file
Additional file 1:
A LatticeBoltzmann scheme for the simulation of diffusion in intracellular crowded systems. (PDF 363 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
AngelesMartinez, L., Theodoropoulos, C. A LatticeBoltzmann scheme for the simulation of diffusion in intracellular crowded systems. BMC Bioinformatics 16, 353 (2015). https://doi.org/10.1186/s1285901507698
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1285901507698
Keywords
 Lattice Boltzmann Method
 Scaled Particle Theory
 Crowding conditions
 Diffusion systems