Applied Mathematics and Computation 217 (2010) 23392347

Contents lists available at ScienceDirect

Applied Mathematics and Computation

The effect of boundary conditions on the mesoscopic lattice Boltzmann method: Case study of a reactiondiffusion based model for Min-protein oscillation
Waipot Ngamsaad a, Paisan Kanthang a, Charin Modchang a, Somchai Sriyab b, Wannapong Triampo a,c,*
a b c

R&D Group of Biological and Environmental Physics, Department of Physics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand Department of Mathematics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand Center of Excellence for Vector and Vector-Borne Diseases, Faculty of Science, Mahidol University, Nakhon Pathom 73170, Thailand

a r t i c l e

i n f o

a b s t r a c t
Min-protein oscillation in Escherichia coli has an essential role in controlling the accurate placement of the cell division septum at the middle-cell zone of the bacteria. This biochemical process has been successfully described by a set of reactiondiffusion equations at the macroscopic level. The lattice Boltzmann method (LBM) has been used to simulate Minprotein oscillation and proved to recover the correct macroscopic equations. In this present work, we studied the effects of LBM boundary conditions (BC) on Min-protein oscillation. The impact of diffusion and reaction dynamics on BCs was also investigated. It was found that the mirror-image BC is a suitable boundary treatment for this Min-protein model. The physical signicance of the results is extensively discussed. 2010 Elsevier Inc. All rights reserved.

Keywords: Lattice Boltzmann method Boundary conditions Min-protein oscillation Cell division

1. Introduction Cell division is a crucial event in the life of every organism. Most bacteria divide symmetrically so that both of the newly formed daughter cells contain a copy of the chromosome. In Escherichia coli (E. coli) and other rod-shaped bacteria, the accurate partitioning process and the proper placement of the septum are precisely restricted to midcell by the effect of nucleoid occlusion [14] and the MinCDE protein system [1,2,5]. The nucleoid occlusion mechanism blocks Z-ring formation from all sites except the nucleoid-free-zones, while the Min system prevents the septum proteins at the cell poles [6]. The Min-protein system consists of MinC, MinD, and MinE proteins [5]. MinC is recruited to the membrane by MinD, and inhibits the formation of the septum proteins in terms of the MinCD complex [710]. MinE is also recruited to the membrane by MinD. In the absence of MinE, MinCD is distributed homogeneously over the entire membrane [11], which blocks septum protein formation at all sites. MinD performs oscillation pole-to-pole, driven by MinE [12]. MinC also colocalizes and cooscillates with MinD [7,8], so that the time-averaged concentration of MinC is lowest at midcell [12]. As MinC inhibits septum protein formation, the septum can assemble only near the midcell area. According to this protein mechanism, Howard et al. [13] proposed the coarse-grained nonlinear reactiondiffusion model for studying the self-organized protein pattern and identifying midcell topological markers. Although this model is simple, it can demonstrate the oscillation pattern and positional information of the cell division process via the Min-protein system;

* Corresponding author at: R&D Group of Biological and Environmental Physics, Department of Physics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand. E-mail addresses: scwtr@mahidol.ac.th, wtriampo@gmail.com (W. Triampo). 0096-3003/$ - see front matter 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.07.033

2340

W. Ngamsaad et al. / Applied Mathematics and Computation 217 (2010) 23392347

however, the dynamics is a continuum version at the macroscopic level. Recently, the mesoscopic dynamics approach assumed the implementation of the lattice Boltzmann method (LBM) [14] for the Min-protein system [15,16]. The method was shown to recover the Howard model at the macroscopic level [15,16]. In addition, the numerical results from this technique were consistent with experimental results [12,17,18]. Hence the LBM approach is an effective alternative tool for studying Min-protein oscillation, which naturally deals with particle dynamics. LBM uses an alternative numerical scheme for simulating uid dynamics. This method evolved from the lattice gas model, which involves many stochastic or random processes resulting in noise or uctuations. It has provided insights into the underling particle dynamics of the physical system, whereas most other approaches focus only on the solution of the macroscopic equation. The advantage of LBM is that it is easy to implement for several shape boundaries due to the nature of particle kinetics. It is perhaps reasonable to say that LBM operates at the mesoscopic level, which is in between the microscopic (lattice gas) and macroscopic (continuum) level. Though LBM consumes tremendously large computing resources, the amount is still relatively small compared to those consumed with other approaches, like molecular dynamic and Monte Carlo simulationswhich could the main reason why LBM has become the more favored method. LBM has been particularly successful in simulating uid ow for a broad variety of complex physical systems, such as multiphase and multi-component uids [19], advectiondispersion [20], and blood ow [2123]. However, the LBM boundary treatments, in some situations, have critical effects on the solution at the macroscopic level. As pointed out by Zhang et al. [24], the standard boundary conditions are not accurate for the class of dispersion transport modeled by LBM. They suggested that the mirror-image method is a suitable boundary treatment to address this problem. The boundary treatment for LBM simulation is an important issue and progress is still being made [2426]. In this present work, we re-examined the LBM boundary effects on Min-protein oscillation. We applied the impermeable boundary condition suggested by Zhang et al. and the mid-grid bounce-back boundary condition [27] to our reactiondiffusion lattice Boltzmann model. We then compared the results from two boundary conditions in both one- and two-dimensional cases. This paper is the rst one to demonstrate the effect of boundary conditions in the LBM reaction diffusion model, especially on the intracellular process. This paper is outlined as follows: In Section 2, we briey review the reactiondiffusion model for MinCDE oscillation by Howard et al. In Section 3, the details for the LMB algorithm and boundary conditions are explained. Finally, the numerical results are discussed in Section 4. 2. Model The dynamics of the interaction between MinD (in the complex with MinC) and MinE proteins have been proposed by Howard et al. [13] through a set of reactiondiffusion equations (RDEs). In this model, they describe the time rates of density change due to the diffusions of MinD and MinE, and the mass transfer between the cell membrane and the cytoplasm. The local mass densities for each Min-protein are denoted by: qD for cytoplasmic MinD, qd for membrane-bound MinD, qE for cytoplasmic MinE, and qe for membrane-bound MinE. In general form, the RDEs can be written as

@ qc @ 2 qc Rc ; Dc @t @x2 a

2:1

where qc qc ~ t is the local mass density of protein species c = {D, d, E, e} at position ~ (with xa as the components) and x; x time t; Rc = Rc({qc}) is the reaction term which depends on the density set of the species {qc}; Dc is the isotropic diffusion coefcient. The reaction terms for each Min-protein are given by:

RD Rd RE

r1 qD
1 r01 qe

r2 qe qd ;

2:2 2:3 2:4 2:5

r1 qD r4 qe

1 r01 qe 1 r04 qD

r2 q e q d ; r3 qD qE ; r3 qD qE :

Re

r4 qe

1 r04 qD

The constant r1 represents the association of MinD to the membrane [11]. r01 corresponds to the membrane-bound MinE, suppressing the recruitment of MinD from the cytoplasm. r2 reects the rate that MinE on the membrane drives MinD on the membrane into the cytoplasm. Based on the evidence of the cytoplasmic interaction between MinD and MinE [28], we let r3 be the rate that cytoplasmic MinD recruits cytoplasmic MinE to the membrane, while r4 corresponds to the rate that MinE disassociates from the membrane to the cytoplasm. Finally, r04 corresponds to the cytoplasmic MinD, suppressing the release of the membrane-bound MinE. The time scale of the diffusion on the membrane is much slower than that in the cytoplasm. It is therefore reasonable to set Dd and De to zero. In this dynamics, we allow for the Min-protein to bind/unbind from the membrane, but not for it to be degraded in the process. Thus, the total amount of each type of Min-protein is conserved. The E. coli cell is a close-system, therefore the appropriated boundary condition is reecting or hard-wall. The zero-ux J  Dc@ qc/@xa = 0 is imposed at the boundary regime.

W. Ngamsaad et al. / Applied Mathematics and Computation 217 (2010) 23392347

2341

We apply the minimal linear stability analysis specically to Min-protein RDEs [29] to determine in which conditions the oscillatory pattern is expected to generate, using the case of the one-dimensional for the sake of simplicity. For practical analysis of other RDEs systems, please refer to [30,31]. The obtained results could lead to parameter estimations, but to do this we rst write our set of Eq. (2.1) in the matrix form:

@~ q @ 2~ q f q D 2 ~~; @t @x

2:6

where ~ is the Min-protein density vector and its elements are qc = qc(x, t); D is a diagonal matrix of diffusion coefcients Dc; q and ~~ is a vector of reaction terms Rc. By perturbing with a small variation vector d~ from the homogeneous state ~ , we f q q q can dene

~  ~ d~; q q q
where q can be obtained by evaluating rst order, we get: ~
@ ~ q @t

2:7

@~ q @xa

~~ 0. Substituting (2.7) into (2.6) and taking a Taylor expansion to the f q

 @~ f where J @~ is the Jacobian matrix evaluated at the x point ~ . Applying the general solution for d~ d~0 ext eiqx to Eq. q q q q ~~ q q (2.8), we obtain the Eigen equation

@d~ q @ 2 d~ q ~ J dq; D @t @x2

2:8

Ad~ xd~; q q
2

2:9

where A  Dq + J*, x is the Eigen-value and q is the wave number. The real part of x will then determine if the equation is linearly stable under a small spatial perturbation, whereas its imaginary part will determine the period of the oscillation,

2p : Imx

2:10

If at least one Eigen-value x has a real part that is positive, then the homogeneous x point is linearly unstable under a small spatial perturbation. Next, we apply linear stability analysis to the reactiondiffusion Eqs. (2.1)(2.5) and obtain the matrix:

1 e 6 6 r1 6 1r0 q 6 1 e 6 A6 r4 r04 q e 6 r3 q E 6 1r04 q 2 4 r r0 q D 4 4 e r3 q E 1r0 q 2 4 D

1r1 q DD q2 r0

r2 q e
r2 q e 0 0

0 0 r3 q DE q2 D

3 r1 r01 q D 0 q 2 r2 qd 1r1 e 7 7 r1 r0 q 1r01qD2 r2 q 7 d7 1 e

r4

r3 q D

1r4 q r0

7: 7 7 7 5

2:11

4 D

We solve the Eigen Eq. (2.9) numerically. First, we use an iterative nite-different method to nd the homogeneous x points of our system of Eq. (2.6). We carried out the iterations several times with different starting points, using parameter values given by Howard et al. [13], and found only one xed point, namely, q 116:873; q 3:36183; q 1383:13; q D E e d 81:6382 (in the unit of lm1). Then, with this xed point, we are able to nd the Eigen-value x of Eq. (2.9). It was found that the positive real part of x is maximized at q = 1.5 lm1. This indicates the presence of a maximum linear unstable oscillating mode, with a wavelength of 4.2 lm [13,16,29]. 3. Lattice Boltzmann method The lattice Boltzmann method (LBM) for reactiondiffusion equations [32] is applied in order to simulate the Min-protein system [15,16]. The lattice Boltzmann equation (LBE) [14,27] with the BhatnagarGrossKrook collision (LBGK) [33] can be written as

x e x; fic ~ ~i Dt; t Dt fic ~ t

Dt

sc

x; fic ~ t fieq;c qc ; ~ DtSc ; u i

3:1

where fic x ; t is the single-particle distribution function of species c associated with lattice velocity e i at position x and time * t; fieq;c qc ; u is the local equilibrium corresponding Maxwellian distribution function with the local density qc and local uid * speed u ; Sc is the source distribution function and sc is the relaxation parameter. The macroscopic quantities can be obtained i from the moments of the distribution function as follows:

X X
i i

fic ~ t x; Sc R c : i

X
i

fieq;c ~ t qc ~ t ; x; x;

3:2 3:3

2342

W. Ngamsaad et al. / Applied Mathematics and Computation 217 (2010) 23392347

We consider the zero mean velocity (~ 0 ow. In addition, to obtain the correct reactiondiffusion equation, we need the u quantities:

~i fieq;c ~ t qc~ 0; e x; u i X X ~i~i fieq;c ~ t ei;a ei;b fieq;c ~ t qc c2 dab ; ee x; x; s i i X c ~i Si 0; e i X ~i~i Sc 0; ee i

3:4 3:5 3:6 3:7

where cs is the speed of sound. The choice of the equilibrium distribution function and the source term is dependent on the constraint in (3.2)(3.7). We chose the equilibrium distribution

fieq;c wi qc ;
where wi is the lattice weight factor. For the source term, we chose the simple form

3:8 3:9

Sc wi Rc ; i

where Rc is macroscopic reaction as dened in (2.1), as well as in (2.2)(2.5). The choices given in (3.8) and (3.9) are the simplest ones that can provide the correct macroscopic equations. We implement the LBM algorithm with two steps: x; x;  Collision step: ~c ~ t fic ~ t Dct fic ~ t fieq;c qc DtSc , f i x; i s c ~ ~c ~ t .  Streaming step: fi x ~i Dt; t Dt f i x; e The boundary condition is computed before streaming the resulting collision distribution function. Here we selected the impermeable boundary condition [24] and the mid-grid bounce-back boundary condition [27], which is equivalent to the hard-wall boundary at the macroscopic level. For simplicity, we set Dx = Dt = 1. Now let us consider the boundary domain * in the 2D system for illustration, located at b x 0; y. The distribution functions in this domain are not permitted to compute the collision step, but the streaming step is allowed. We dene the adjacent uid domain closest to the boundary do* main located at g x 1; y. If the boundary condition is mid-grid bounce-back (we will call this simply bounce-back), we calculate the unknown distribution functions as

fic b fic b ei ;
* * e i Dt

3:10

where b must be inside the uid domain. The accuracy of this boundary condition is known as second order [27]. The impermeable boundary condition via the mirror-image method (we will call this simply mirror-image) [24] can be calculated by

fic b fkc g ;

3:11

Fig. 1. The ow chart demonstrates LBM implementation.

W. Ngamsaad et al. / Applied Mathematics and Computation 217 (2010) 23392347

2343

where index k is the opposite direction of i corresponding to lattice velocity e i . The owchart of LBM implementation is summarized in Fig. 1. It has been shown, through the Chapman-Enskog analysis, that the continuum limit of LBE (3.1) recovers the correct reactiondiffusion equations [32]. For the D1Q3 and D2Q9 lattice model, the reactiondiffusion equation (up to the second order correction) is given by [15,16]

@ qc @ 2 qc Rc O @ 3 ; Dc 2 @t @xa
where Dc c2 sc Dt=2 is the diffusion coefcient. s

3:12

Fig. 2. The density prole q(x, t) arrived at from LBM simulation with two different boundary conditions, bounce-back (s) and mirror-image (), compared to the exact solution (dashed line) from Eq. (4.1). The LBM results were simulated on 100 grids with D = 0.0028 and data was selected at iteration 10.

Fig. 3. The partial time evolution of normalized-densities qD+d (left) and qE+e (right) at the selected positions, including x = 2(top), x = 99 (bottom) and x = 50 (middle), for the two polar zones and middle-cell zone, respectively.

2344

W. Ngamsaad et al. / Applied Mathematics and Computation 217 (2010) 23392347

4. Numerical results and discussion We simulated the dynamics of Min-proteins oscillation with LBM, as described in the previous section, in both one- and two-dimensional systems. For the one-dimensional system, we selected a one-dimensional lattice model with three allowed n* o velocities (D1Q3) where the lattice velocities are dened by e i f0; v ; v g, and i = {0, 1, 2} is the lattice link index. The weight factors for this case are given by w0 = 2/3 and w1 = w2 = 1/6. For two-dimensional system, we selected a two-dimensional square lattice model with nine allowed velocities (D2Q9) where the lattice velocities are dened by n* o e i f0; 0; v ; 0; 0; v ; v ; v g, and i = {0, 1, 2, . . . , 8} is the lattice link index. The weight factors for this case are given by w0 = 4/9, w1 = w2 = w3 = w4 = 1/9 and w5 = w6 = w7 = w8 = 1/36. Here, v  Dx/Dt is the lattice speed magnitude; Dx and Dt are the lattice spacing and time step intervals, respectively. Usually, we selected Dx = Dt = 1, which leads to v = 1 p (for 2D, Dy = Dx). Therefore our system is dimensionless. The (dimensionless) speed of sound is cs 1=3 for both D1Q3 and D2Q9 [34]. The relaxation time is calculated from sc = 3Dc + 0.5 for both cases. We implemented the LBM on a PC using C programming. First, we validated our LBM on the pure diffusion equation by dropping the reaction term (by setting Rc = 0) in the onedimensional system. The system was simulated on 100 grids, which represented a bacterium 2 lm long. In this situation, we simulated only a single species of density: q(x, t)  qc(x, t). The dimensionless diffusion constant was set to be D  Dc = 0.0028. We initialized the density with a delta-like function: q(x0 = 50, t = 0)  q0 = 1000. We show the density prole q(x, t) from our LBM simulation with two different boundary conditions: bounce-back and mirror-image. We also compare these results with the exact solution:

q0 x x0 2 : qx; t p exp 4Dt 4pDt

"

4:1

The results are shown in Fig. 2. It was seen that both numerical results from the LBM and the exact solution (4.1) are in good agreement. This indicates that the two different boundary conditions, namely bounce-back and mirror-image, have no effect whatsoever on the pure

Fig. 4. The normalized-time-averaged MinD (left) and MinE (right) densities in 1D, relative to their respective time-average maxima as functions of the position along the bacterium cell, which compares the results from LBM with two different boundary conditions: bounce-back (s) and mirror-image (). Each sub-plot magnies as non-smooth due to the boundary effects at the polar zone.

W. Ngamsaad et al. / Applied Mathematics and Computation 217 (2010) 23392347

2345

diffusion equation. However, the boundary conditions will affect the system when we take the reaction terms into account, as we will see in the next simulation. It should be noted that since reaction terms of any given reactiondiffusion based model can greatly differ, depending on model details and dynamics, great care must be taken when making a general or universal statement about the results. To simulate Min-protein oscillation in a one-dimensional system, we initialized the distribution function from an equilibrium state where MinD and MinE were randomly distributed as 3000 for qD, 170 for qE, and qd = qe = 0. We used the diffusion constants and the reaction rate parameters as used in Ref. [13], which were tested for linear stability following Section 2. However, to put them into the LBM, we need to make all parameters dimensionless. The dimensionless parameters we chose are DD = 0.28, DE = 0.6, r1 = 8 103, r2 = 2.52 102, r3 = 1.6 101, r4 = 3.2 104, r01 280, and r04 270. Because no diffusion occurs on the membrane, we set Dd = De = 0. We simulated the system up to 2.5 107 iterations. We computed the total concentration in the cytoplasm and on the membrane of each Min-protein, qD+d = qD + qd for MinD and qE+e = qE + qe for MinE, by excluding the boundary node. The partial oscillations of Min-proteins in the one-dimensional system are shown in Fig. 3. Here, we used qD+d and qE+e at the selected positions, including x = 2, x = 99 (at the polar zones) and x = 50 (in the middle-cell zone), and then we normalized them. The (normalized) densities of Min-proteins evolve from the random-distributed initial state to the spatiotemporal oscillation state with some transient time (see Fig. 3). We noticed that the transient times of the result when using the mirror-image boundary condition are faster than when using the bounce-back boundary condition, approximately 1.5 times faster for all sets. To observe the spatial topology of Min-proteins, we calculated the normalized-time-averaged MinD and MinE densities as functions of position, as shown in Fig. 4. The MinE normalized-averaged-density peaks at the cell middle and has a minimum at the cell rims, with MinD being nearly out of phase with MinE. This implies that MinE tries to prevent septum proteins from being formed at the cell middle site and tries to suppress the inhibitor, MinD, from this region. To validate this conclusion, we compared the numerical results to the experimental data compiled by our group members [17,18]. When we compared the LBM results with the experimental results, they were in qualitative agreement. These results are also in qualitative agreement with the experimental data by another group [12], as well as the numerical results obtained in another work [13].

Fig. 5. The normalized-time-averaged MinD (left column) and MinE (right column) densities in 2D, relative to their respective time-average maxima as functions of the position (x, y) inside the bacterium cell. Comparison results from LBM with two different boundary conditions: bounce-back (upper row) and mirror-image (lower row).

2346

W. Ngamsaad et al. / Applied Mathematics and Computation 217 (2010) 23392347

We then extended our simulation to a two-dimensional system. The system was simulated on 50 100 grids, which represented typically a bacterium 1 2 lm in size. We initialized the distribution function from an equilibrium state where MinD and MinE were randomly distributed as 3000 for qD, 170 for qE, and qd = qe = 0 for each grid point. The dimensionless parameters were chosen in the same way as in the one-dimensional case. The normalized-time-averaged MinD and MinE densities as functions of position are shown in Fig. 5 for this case. It was found that the average densities of Min-proteins are nearly homogenous in x-direction; therefore the dominant dynamics oscillation of Min-proteins were expressed in ydirection. Due to this reason, the dynamics in 2D are exactly identical to that in 1D. Regarding the effect of boundary conditions, namely bounce-back and mirror-image, on dynamic oscillation, differences in simulation results via LBM were noticeably found. We observed in the Min-protein densities at the boundary node a portion of a non-smooth curve when using the bounce-back boundary conditions, whereas the densities were found to be very smooth at the boundary node when using the mirror-image boundary condition for both one- and two-dimensional systems (see Figs. 4 and 5). This happened for both MinD and MinE, but was more prominent in the former. The difference in the degree of non-smoothness is mainly due to the difference in protein copy concentrations. Interestingly, it evidently appears that reaction terms respond differently to different boundary conditions. This may suggest that one should take greater care in boundary conditions than one might expect when applying LBM for the RDE based model, especially for the model where reactions crucially occur at the boundary. It is worth mentioning that there are several published methods for boundary treatments via LBM for the hydrodynamics system. Among them, the bounce-back seems to be the most popular for non-slip type boundary implementation. This perhaps is due to the fact that it ensures that the components of the uid velocity vanish at the boundary. However, this method cannot be applied to solute transport, especially the reactiondiffusion system, because the nature of dispersion-transport boundaries differs from hydrodynamics boundaries [27]. The diffusive system requires that the ux vanish (J  Dc@ qc/ @xa = 0) at the boundary and the mirror-image implementation satises this constraint [35]. Lastly, it is obviously seen from our results that the mirror-image boundary condition is better than the bounce-back boundary condition for the simulation of this model via the LBM. 5. Concluding remarks In this paper, we investigated the effects that boundary conditions have on the LBM when simulating the dynamic oscillation of Min-proteins and describing MinD/MinE interaction [15,16]. We applied the mirror-image boundary condition suggested by Zhang et al. and the mid-grid bounce-back boundary condition [24] to our reactiondiffusion lattice Boltzmann model. We compared the results from two boundary conditions and found that the mirror-image method is a suitable boundary condition for this Min-protein problem, as far as the smoothness of the data curves is concerned. The degree of non-smoothness varied, depending on Min-protein species or concentrations. However, more work is required to understand which boundary condition is the most appropriate for a given model. It is reasonable to say that the LBM approach could provide an alternative fast computational tool for studying Min-protein oscillation. It provides numerical results which are in qualitative agreement with experimental results [12,17,18]. However, optimization treatment, especially regarding boundary conditions, could be a very important issue. Still, we believe that the LBM is a useful scheme for simulating at the cellular level those biological systems governed by reactiondiffusion equations. Acknowledgements We are grateful to the anonymous reviewer for truly helpful comments. The authors thank Mr. David Blyler for reading and editing the manuscript. This work was supported by the National Center for Genetic Engineering and Biotechnology (BIOTEC), the Thailand Research Fund (TRF), the Commission on Higher Education (CHE), the Development and Promotion of Science and Technology Talents program (DPST), and the Thailand Center of Excellence in Physics. References
[1] J. Lutkenhaus, S.G. Addinall, Bacterial cell division and the Z ring, Annu. Rev. Biochem. 66 (1997) 93116. [2] L. Rotheld, S. Justice, J. Garcia-Lara, Bacterial cell division, Annu. Rev. Genet. 33 (1999) 423448. [3] C.L. Woldringh, E. Mulder, P.G. Huls, N. Vischer, Toporegulation of bacterial division according to the nucleoid occlusion model, Res. Microbiol. 142 (1991) 309320. [4] X.C. Yu, W. Margolin, FtsZ ring clusters in min and partition mutants: role of both the Min system and the nucleoid in regulating FtsZ ring localization, Mol. Microbiol. 32 (1999) 315326. [5] P.A. de Boer, R.E. Crossley, L.I. Rotheld, A division inhibitor and a topological specicity factor coded for by the minicell locus determine proper placement of the division septum in E. coli, Cell 56 (1989) 641649. [6] K. Suefuji, R. Valluzzi, D. Ray Chaudhuri, Dynamic assembly of MinD into lament bundles modulated by ATP, phospholipids, and MinE, Proc. Natl. Acad. Sci. USA. 99 (2002) 1677616781. [7] P.A. de Boer, R.E. Crossley, A.R. Hand, L.I. Rotheld, The MinD protein is a membrane ATPase required for the correct placement of the Escherichia coli division site, EMBO J. 10 (1991) 43714380. [8] D.M. Raskin, P.A. de Boer, MinDE-Dependent pole-to-pole oscillation of division inhibitor MinC in Escherichia coli, J. Bacteriol. 181 (1999) 64196424. [9] Z. Hu, J. Lutkenhaus, Topological regulation of cell division in Escherichia coli involves rapid pole to pole oscillation of the division inhibitor MinC under the control of MinD and MinE, Mol. Microbiol. 34 (1999) 8290.

W. Ngamsaad et al. / Applied Mathematics and Computation 217 (2010) 23392347

2347

[10] A.L. Marston, J. Errington, Selection of the midcell division site in Bacillus subtilis through MinD-dependent polar localization and activation of MinC, Mol. Microbiol. 33 (1999) 8496. [11] S.L. Rowland, X. Fu, M.A. Sayed, Y. Zhang, W.R. Cook, L.I. Rotheld, Membrane redistribution of the Escherichia coli MinD protein induced by MinE, J. Bacteriol. 182 (2000) 613619. [12] C.A. Hale, H. Meinhardt, P.A.J. de Boer, Dynamic localization cycle of the cell division regulator MinE in Escherichia coli, EMBO J. 20 (2001) 15631572. [13] M. Howard, A.D. Rutenberg, S. de Vet, Dynamic compartmentalization of bacteria: accurate division in E. coli, Phys. Rev. Lett. 87 (2001) 278102. [14] S. Chen, G.D. Doolen, Lattice Boltzmann method for uid ows, Annu. Rev. Fluid Mech. 30 (1998) 329364. [15] W. Ngamsaad, W. Triampo, P. Kanthang, I.-M. Tang, N. Nuttawut, C. Modjung, Y. Lenbury, A lattice Boltzmann method for modeling the dynamic poleto-pole oscillations of Min proteins for determining the position of the midcell division plane, J. Korean Phys. Soc. 46 (2005) 10251030. [16] S. Sriyab, J. Yojina, W. Ngamsaad, P. Kangthang, C. Modchang, N. Nuttavut, Y. Lenbury, C. Krittanai, W. Triampo, Mesoscale modeling technique for studying the dynamics oscillation of Min protein: pattern formation analysis with lattice Boltzmann method, Comput. Biol. Med. 39 (2009) 412424. [17] U. Junthorn, S. Unai, P. Kangthang, W. Ngamsaad, C. Modchang, W. Triampo, C. Krittanai, D. Triampo, Y. Lenbury, Single-particle tracking method for quantitative tracking and biophysical studies of the MinE protein, J. Korean Phys. Soc. 52 (2008) 639648. [18] S. Unai, P. Kanthang, U. Junthon, W. Ngamsaad, W. Triampo, C. Modchang, C. Krittanai, Quantitative analysis of time-series uorescence microscopy using a spot tracking method: application to Min protein dynamics, Biologia 64 (2009) 2742. [19] N.S. Martys, H.D. Chen, Simulation of multicomponent uids in complex three-dimensional geometries by the lattice Boltzmann method, Phys. Rev. E 53 (1996) 743750. [20] R.G.M. van der Sman, M.H. Ernst, Advectiondiffusion lattice Boltzmann scheme for irregular lattices, J. Comput. Phys. 60 (2000) 766782. [21] M. Hirabayashi, M. Ohta, D.A. Rufenacht, B. Chopard, A lattice Boltzmann study of blood ow in stented aneurism, FGCS 20 (2004) 925934. [22] C. Migliorini, Y.H. Qian, H.H. Chen, E. Brown, R. Jain, L. Munn, Red blood cells augment leukocyte rolling in a virtual blood vessel, Biophys. J. 84 (2002) 18341841. [23] C.H. Sun, C. Migliorini, L. Munn, Red blood cells initiate leukocyte rolling in postcapillary expansions: a lattice Boltzmann analysis, Biophys. J. 85 (2003) 208222. [24] X. Zhang, J.W. Crawford, A.G. Bengough, I.M. Young, On boundary conditions in the lattice Boltzmann model for advection and anisotropic dispersion equation, Adv. Water Resour. 25 (2002) 601609. [25] S. Chen, D.O. Martinez, R. Mei, On boundary conditions in lattice Boltzmann methods, Phys. Fluids 8 (1996) 25272536. [26] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (1997) 15911598. [27] S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Oxford University Press, New York, 2001. p. 82. [28] J. Huang, C. Cao, J. Lutkenhaus, Interaction between FtsZ and inhibitors of cell division, J. Bacteriol. 178 (1996) 50805085. [29] C. Modchang, W. Triampo, P. Kangthang, U. Junthorn, S. Unai, W. Ngamsaad, N. Nuttavut, D. Triampo, Y. Lenbury, Stochastic modeling of external electric eld effect on Escherichia coli Min protein dynamics, J. Korean Phys. Soc. 53 (2008) 851862. [30] Rui Dilo, Turing instabilities and patterns near a Hopf bifurcation, Appl. Math. Comput. 164 (2005) 391414. [31] V. Gaychuk, B. Datsko, Inhomogeneous oscillatory solutions in fractional reactiondiffusion systems and their computer modeling, Appl. Math. Comput. 198 (2008) 251260. [32] S.P. Dawson, S. Chen, G.D. Doolen, Lattice Boltzmann computations for reactiondiffusion equations, J. Chem. Phys. 98 (1993) 5141523. [33] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511525. [34] Y.H. Qian, D. dHumieres, P.A. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17 (1992) 479484. [35] G. Drazer, J. Koplik, Tracer dispersion in two-dimensional rough fractures, Phys. Rev. E 056104-1 (2001) 11.