Computational analysis of the influence of the model boundary conditions on the bacterial self-organization*

Abstract

This paper deals with the effects of the type of boundary conditions on the spatiotemporal pattern formation in the computational modelling the bacterial pattern formation in a onedimensional-in-space domain. The computational model was derived from a mathematical model used in another research. By running tests with different boundary conditions, the output results were analyzed with a special emphasis on the edges of formed patterns. The numerical simulation, based on the governing equations of the reaction-diffusion-chemotaxis type, was carried out using the finite difference technique. The developed numerical simulator was validated by using published experimental data and known numerical solutions.

Keywords

• eol>Chemotaxis
• Pattern formation
• Mathematical modeling

1. Introduction

Organisms are affected by chemical stimuli in their environment. Whether it is nutritious and attractive or harmful and repulsive, movement based on such stimuli is called chemotaxis [23]. One of the most analyzed organisms when dealing with chemotaxis is Escherichia coli. E. coli exhibit chemotaxis to self-excreted attractant which results in the swarming of bacteria [2, 24]. Cultures of E. coli are capable of self-organization forming growth patterns of millimeter-scale and exhibiting strong inhomogeneities in their density, mainly near the contact lines and surfaces [2].

Starting from pioneer work of Keller and Segel [5], the mathematical modelling plays a substantial role in understanding the chemotaxis [7, 27]. Although various mathematical models based on partial differential equations have been developed, the system introduced by Keller and Segel remains among the most widely used [6, 7, 18, 27].

According to the Keller and Segel approach [5], the dynamics of the bacteria and chemoattractant are mathematically modelled by a system of nonlinear equations of the reaction-diffusion-chemotaxis type [7]. The processes in bacterial cultures are rather complicated and need to be modelled in threedimensional space [3, 14, 27, 28]. Nevertheless, due to the accumulation of cells near the contact line at the top surface of the fluid cultures the essentially three-dimensional processes may be approximated in one-dimension [1, 19].

The bioluminescence pattern formation in suspensions of lux gene engineered E. coli in right circular containers (test-tubes) has been experimentally and numerically studied [1, 8, 9, 28]. When modeling the bacterial migration near the three-phase contact line, the one-dimensional in space mathematical model is usually defined in an interval – a continuous circle and periodic boundary conditions are applied for boundaries of the interval [1, 26, 28]. One dimensional modelling can be also used for central axis of the circular container, using different boundary conditions [14].

The aim of this work was to investigate the effects of the type of the boundary conditions on the spatiotemporal pattern formation in the computational modelling the bacterial self-organization in a one-dimensional domain. Two types of the boundary conditions were applied in the analysis, the zero flux (non-leakage) and the Dirichlet boundary condition when the bacterial density is held at a fixed value. The numerical simulation, based on the governing equations of the reaction-diffusionchemotaxis type [1, 5, 6, 7], was carried out using the finite difference technique [29].

2. Mathematical model

The dynamics of bacterial population and chemoattractant is modeled mathematically by a coupled system of the reaction-diffusion-chemotaxis equations together with some initial and boundary conditions [1, 5, 6, 7]. 2.1.

Governing equations

Keller and Segel proposed two conservation equations to model the biological and physical processes in a bacterial population [5, 6, 7, 13]: ∂ n ∂ t ∂ c ∂ t =∇ ( Dn ∇ n−h ( n , c ) n ∇ c )+ f ( n , c ) , ( 1 )

=∇ ( Dc ∇ c )+ g p ( n , c ) n−gd ( n , c ) c , x∈ Ω , t >0 , where x is space, t is time, n(x, t) is cell density, c(x, t) is the chemoattractant concentration, Dn and Dc are the diffusion coefficients (assumed to be constant), f(n, c) is cell growth and death, h(n, c) is the chemotactic sensitivity, g p is the production of the chemoattractant and gd is the degradation of the chemoattractant [5].

For this chemotaxis model f, g p, gd, h variables are instantiated with concrete expressions. Cell growth is often assumed to be logistic f ( n , c )=k3 n ( 1−n / n0 ) [7], where k3 is the growth rate of cell population and n0 is the cell density under steady-state conditions (both constant) [2]. Chemoattractant production is expressed as g p ( n , c )=k 4 n ( k5+n2 ) specifically for E. coli rapid increase in attractant production (introduced by Tyson et al. [13]), where values of k 4 , k5 are not exactly known. The degradation of chemoattractant is marked as gd ( n , c )=k6, where k6 is also not exactly known [13]. Chemotactic sensitivity is noted as h ( n , c )=k1/( k 2+ c )2 (expression derived by Lapidus and Schiller), where k1 and k2 are constants [10, 13, 20].

The mathematical model was applied to investigate the bacterial self-organization in a small, rounded container as detected by bioluminescence imaging [8, 9]. For modeling the experimentally observed space-time plots of quasi-one-dimensional bioluminescence measured along the three-phase contact line, the mathematical model was defined in a one-dimensional domain – the circumference of the vessel. Applying the instantiated expressions of the functions f , g p, gd and h, the following system of partial differential equations describes the bacterial self-organization:

= Dn ∆n−∇ ( ∂ n ∂ t

k1 n ( k2+ c ) 2 ∇ c)+ k3 n(1− n ), n 0 ( 2 ) ∂ c ∂ t = Dc ∆c + k 4 n2 k5+n2 −k6 c , x∈ (0 , l ) , t >0 , where l is the length of an interval, ∆ is the Laplace operator in the one-dimensional Cartesian coordinate. In the case of modelling the bacterial self-organization along the three-phase contact line [1, 26, 28], l corresponds to the circumference of the vessel (the length of the contact line).

Initial and boundary conditions

The initial cell and chemoattractant distribution are assumed to be non-uniform: n ( x ,0)=n0 x ( x ) , c ( x ,0)=c0 x ( x ) , x∈ [ 0 , l ] , ( 3 ) where n0 x ( x ) and c0 x ( x ) are the initial (t = 0) cell density and chemoattractant concentration, respectively.

Three types of boundary conditions were analyzed to model bacterial self-organization in a cylindrical test-tube container (vessel).

Periodic boundary conditions are used to model the bacterial self-organization at the edge (a continuous circle) of the top surface of a cylindrical container (t > 0), n (0 , t )=n (l , t ) , c (0 , t )=c (l , t ) , ∂∂ nx ¿x=0= ∂∂ nx ¿x=l , ∂∂ cx ¿x=0= ∂∂ cx ¿x=l , ( 4 ) where l is the circumference of the cylinder. The bacterial pattern formation applying the periodic boundary conditions ( 4 ) have been already investigated [1, 26, 28].

Zero flux boundary conditions are used on both sides of the interval [0, l] in modelling the onedimensional vertical bacterial self-organization in the central axis of the vessel, ∂ n ∂ n ∂ c ∂ c ∂ x ¿x=0=0 , ∂ x ¿x=l=0 , ∂ x ¿x=0=0 , ∂ x ¿x=l=0. where l is the height of the cylinder. The mathematical model defined by ( 2 ), ( 3 ) and ( 5 ) can be also used to model the bacterial self-organization in a transverse section of a glass channel [14]. Zero flux (non-leakage) boundary conditions defined on both sides of the interval restrict any migration of bacteria outside the interval [21].

The Dirichlet boundary condition is used when the boundary is held at a fixed density of cells[22]. Assuming that the upper part of the liquid bacterial culture is stirred, a thin bottom layer usually remains unstirred [4]. The density of bacteria in the stirring layer can be assumed constant if the stirring layer is relatively tick in comparison with the stagnant layer. Since the bacterial density remains constant above the stagnant layer, the Dirichlet boundary condition can be applied on the upper boundary of the stagnant layer [22]. Thus, mixed boundary conditions, the zero-flux and the Dirichlet, are used when modelling the bacterial self-organization in the central axis of the vessel [22]: | = ∂ n ∂ x x=0

| ∂ c ∂ x x=0

=0 , n (l , t )=a , c (l , t )=b , where l is the height of the stagnant layer, a and b are the constant concentrations of cells and chemoattractant in the bulk, respectively. This condition simulates bacterial movement in a constantly refreshed and mixed solution, that has a constant bacterial density and chemoattractant concentration. ( 5 ) ( 6 )

Dimensionless model

A dimensionless mathematical model is derived from the mathematical model ( 2 ) – ( 4 ) to define the main governing parameters: ( 7 ) u= n , a¿= a , v = n0 n0

D=

D

n , χ = Dc ∂ u ∂2 u ∂ χu ∂ v ∂ t = D ∂ x2 − ∂ x ( (1+ αν )2 ∂ x

)+ γru (1−u ) , ( 8 ) ∂ v ∂2 v

= ∂ t ∂ x2

Here u is the dimensionless cell density, v is the dimensionless chemoattractant concentration, α is the receptor sensitivity, β is the suturing of the signal production, and γ is the spatial and temporal scale.

The following initial conditions for the dimensionless model were used:

u ( x ,0)=1+ ε ( x ) , v ( x ,0)=0 , x∈ [ 0 ,1] , where ε ( x ) is a 20% random uniform spatial perturbation.

The periodic boundary conditions ( 4 ) for the dimensionless model becomes:

u (0 , t )=u (1 , t ) , v (0 , t )=v (1 , t ) , ∂∂ ux ¿x=0= ∂∂ ux ¿x=1 , ∂∂ vx ¿x=0= ∂∂ vx ¿x=1 .

The zero flux boundary conditions ( 5 ) for the dimensionless model becomes: The mixed boundary conditions ( 6 ) take the following dimensionless form: ∂ u ∂ u ∂ v ∂ v ∂ x ¿x=0=0 , ∂ x ¿x=1=0 , ∂ x ¿x=0=0 , ∂ x ¿x=1=0.

| = ∂ u ∂ v ∂ x x=0

| ∂ x x=0 =0 , n (1 , t )=a , c (1 , t )=b .

( 9 ) ( 10 ) ( 11 ) ( 12 )

3. Numerical simulation

To solve the initial boundary value problem ( 2 )-( 12 ) numerically a uniform discrete grid in space and time was introduced. An explicit finite difference scheme was developed as a result of the difference approximation [30]. The governing equations ( 2 ) were approximated as follows: vi , j+1=vi , j+ ∆ t ( vi+1, j−2 vi , j+ vi−1, j + γ ( ∆ x2

2 1+ β ui2, j −vi , j)),

ui , j i=1 , … , N −1 , j=1 , … , M −1 , where ui,j  u(xi, tj), xi = xi-1 + x, tj = tj-1 + t, for i = 1,…, N and j = 1,…, M, x0 = 0, xN = 1, t0 = 0, tM = T, where T is the dimensionless duration of simulated process [29]. The initial ( 9 ) and boundary conditions ( 10 )-( 12 ) were approximated accordingly. To obtain an accurate and stable numerical solution, the dimensionless size of time step t varied between 10-7 and 10-6 depending on the boundary conditions. The size of space step x varied from 0.002 to 0.004.

The numerical simulator was developed in Python language. The numerical solution of the mathematical model was validated by using published experimental data and known numerical solutions [1, 8, 9, 26, 28]. The following typical values of the model parameters were used in all the numerical experiments [1]:

Values ( 14 ) of the model parameters were determined experimentally by changing input parameters and aiming to simulate the spatiotemporal patterns comparable to those observed experimentally in the liquid cultures of luminous E. coli [8, 9]. The patterns were fitted semi-formally by counting the spatial and temporal accumulations.

For an analysis of the process dynamics the average cell density ū(t) and chemoattractant concentration v1(t) were calculated as functions of time [1],

D=0.1 , χ =9.2 , r =1 , α =0.7 , β =1.4 , γ =625. ( 14 ) ( 13 ) ( 15 ) ( 16 )

To evaluate the bacterial average distribution in a space, the cells density u(x, t) was integrated over time and then averaged, where ~u stands for average cell density, taken across all the time moments. The function ~u(x) expresses the generally favored bacteria distribution and any specific accumulation of bacteria, especially in the edges.

4. Results and discussion

To investigate the effects of the type of the boundary conditions on the spatiotemporal pattern formation the bacterial self-organization was simulated using two types of the boundary conditions, the zero flux ( 11 ) and the mixed ( 12 ) boundary conditions. Different boundary conditions correspond to different experimental conditions. The simulated patterns were compared with each other, with experimental data [8, 9] and with known simulated patterns obtained at the periodic boundary conditions ( 10 ) [1, 28]. 4.1.

Zero flux boundaries

The zero flux (non-leakage) boundary conditions ( 11 ) were used on both sides of the space interval. The simulation results are depicted in Fig. 1.

a) c)

Fig. 1. Simulated space-time plots of the dimensionless cell density u (a) ant chemoattractant concentration v (b), the cell distribution ~u averaged over the time (c) and the corresponding values ū and v1averaged on the space interval (d). The simulation was performed using the zero flux boundary conditions ( 11 ).

The simulated spatiotemporal patterns (Fig. 1a and Fig. 1b) are very similar to the corresponding experimental patterns [8, 9] and those simulated using the periodic boundary conditions [1, 28]. The evolution of the average cell density and the chemoattractant concentration (Fig. 1c) are practically identical to those simulated at periodic boundary conditions [8, 9].

Although, the spatiotemporal patterns simulated using the zero flux boundary conditions ( 11 ) are very similar to those obtained applying the periodic boundary conditions ( 10 ), one can see noticeable difference near the interval borders x = 0 and x = 1. Fig. 1c easily shows accumulations of bacteria at those edges of the interval, while in the case of the periodic boundary conditions no specific accumulations were observed due to the continuity of the circle [8, 9].

As the application of the zero flux boundary conditions corresponds to modelling the vertical bacterial self-organization in the central axis of the vessel, Fig. 1c shows that the bacterial population tends to accumulates near the boundaries that restrict the cell migration outside the vessel. Similar peaks of the bacterial densities have also been observed in snapshots of cell density on the inner lateral surface of the vessel in two-dimensional-in-space modelling the bacterial self-organization [30]. 4.2.

High external bacterial density

The Dirichlet boundary condition was applied on the right boundary of the interval (x = 1) assuming the other boundary as zero flux boundary (x = 0). These mixed boundary conditions ( 12 ) are applied when modelling the bacterial self-organization in the central axis of the vessel keeping the constant bacterial density at the top boundary (x = 1). Fig. 2 presents the simulation results at the outer bacterial density a = 1 and the outer concentration b = 1 of the chemoattractant.

Fig. 2. Simulated space-time plots of the dimensionless cell density u (a) and chemoattractant concentration v (b), the cell distribution ~u averaged over the time (c) and the corresponding values ū and v1averaged on the space interval (d). The simulation was performed using the mixed boundary conditions ( 12 ) at a = 1 and b = 1.

Fig. 2 shows that the higher right-side values of the cell density and the chemoattractant concentration do not last far dropping instantly to below background values in places next to the boundary. Consequently, the lower values also mean that bacteria close to the edge do is not incentivized to stay long and that pushes the whole structure away from that boundary. The structure in the middle shows some signs of branching out to the right side, but those branches quickly fade away, as getting close to the edge brings worse conditions for the bacteria.

The large spatial changes in the bacterial density can be explained by large difference between the outer bacterial density (u(1, t) = a = 1) and the average inner density ū(t) when t > 1 (Fig. 2d). Fig. 2d also shows large difference in the chemoattractant concentration. These differences lead to formation of patterns (Figs. 2a and 2b) very different from the patterns simulated at zero flux boundary conditions (Figs. 1a and 1b).

Moderate external bacterial density

Since the patterns, simulated using mixed boundary conditions ( 12 ) at relatively high external cell density (a = 1) and zero concentration (b = 0) of the chemoattractant, extremely differs from the patterns simulated using zero flux boundary conditions ( 11 ) and using periodic boundary conditions ( 10 ) [1, 28], the bacterial self-organization was also simulated at values of a and b comparable with the average (background) values of the bacterial density u and the chemoattractant v, respectively. The simulation results at a = 0.7 and b = 0.25 are depicted in Fig. 3.

Fig. 3. Simulated space-time plot of the dimensionless cell density u (a) and chemoattractant concentration v (b), the cell distribution ~u averaged over the time (c) and the corresponding values ū and v1averaged on the space interval (d). The simulation was performed using the mixed boundary conditions ( 12 ) at a = 0.7 and b = 0.25.

One can see that the simulated spatiotemporal patterns depicted in Figs. 3a and 3b are very different from the patterns shown in Figs. 2a and 2b though the same mixed boundary conditions ( 12 ) were used. Only values of a and b were different.

Values of a and b similar to the background provide stability to the bacterial structure, similar to the results obtained using the zero flux boundary conditions (Fig. 1). Only a slight difference can be noticed near the right (outer) boundary, x ≥ 0.95. Further form this boundary (x < 0.95), the spatiotemporal patterns and the averaged values are approximately the same as in the case of zero flux boundary conditions (Fig. 2).

Zero-flux boundary conditions are usually used for theoretical analysis of bacterial pattern formation [7, 11, 18, 19]. The results of the numerical simulation obtained with different boundary conditions show that type of boundary conditions may be crucial for pattern formation and should be taken into consideration in the theoretical analysis of the chemotaxis model.

5. Conclusions

We have shown that numerical simulation based on the Keller–Segel mathematical model ( 2 ) can be used as a tool to study the formation of patterns representing the self-organization of the bacteria. An application of different boundary conditions ( 10 )-( 12 ) correspond to different experimental conditions.

The application of the zero flux boundary conditions ( 11 ) can be used to model the vertical bacterial self-organization in the vessel, where the bacterial population tends to accumulates near the boundaries that restrict the cell movement outside the vessel (Fig. 1).

The simulated spatiotemporal patterns inside a one-dimensional-in-space domain (an interval) are similar for all three types of the boundary conditions, periodic ( 10 ), zero flux ( 11 ) and mixed ( 12 ), only if the outer fixed bacterial density a and the chemoattractant concentration b are similar to the corresponding average values.

The spatiotemporal patterns change boldly (Fig. 2) when the mixed boundary conditions ( 12 ) are used at values of a and b notably different from the corresponding average values. [18] K.J. Painter, T. Hillen, Spatio-temporal chaos in a chemotactic model, Physica D, 2011, 240( 4–5 ), pp. 363–375. [19] M.R. Myerscough, P.K. Maini, K.J. Painter, Pattern formation in a generalized chemotactic model,

Bull. Math. Biol, 1998., 60( 1 ), pp. 1–26. [20] D.E. Woodward, R. Tyson, M.R. Myerscough, J.D. Murray, E.O. Budrene, H.C. Berg, Spatiotemporal patterns generated by Salmonella typhimurium, Biophys. J., 1995, 68( 5 ), pp. 2181– 2189. [21] M. Di Francesco, A. Lorz and P. A. Markowich, Chemotaxis–fluid coupled model for swimming bacteria with nonlinear diffusion: global existence and asymptotic behavior, Discrete Contin.

Dyn. Syst., 2010, 28( 4 ), pp. 1437–1453. [22] Yulan Wang, Michael Winkler, Zhaoyin Xiang. Global masspreserving solutions to a chemotaxisfluid model involving Dirichlet boundary conditions for the signal. Anal. Appl. 2022, Vol. 20, No. 01, pp. 141-170. [23] M. Eisenbach. Chemotaxis. London, England: Imperial College Press, 2004. [24] M. P. Brenner, L. S. Levitov, E. O. Budrene. Physical mechanisms for chemotactic pattern formation by bacteria. Biophys. J. 1998, 74( 4 ), p. 1677–1693. [25] T. Hillen, J. Zielinski, K. J. Painter. Mergingemerging systems can describe spatiotemporal patterning in a chemotaxis model. Discrete Contin. Dyn. Syst., 2013, 18( 10 ), p. 2513–2536. [26] L. Litvinas, R. Baronas, R. Šimkus. Švytinčių bakterijų kolonijos struktūros formavimosi kompiuterinis modeliavimas. XV kompiuterininkų konferencijos mokslo darbai, Žara, 2011, p. 99-107. [27] Kevin J. Painter. Mathematical models for chemotaxis and their applications in self-organisation phenomena. J. Theor. Biol., 2019, 481:162–182. [28] R. Baronas, Z. Ledas, and R. Šimkus. Computational modeling of the bacterial self-organization in a rounded container: The effect of dimensionality. Nonlinear Anal. Model. Control, 2015, 20( 4 ):603–620. [29] A.A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, Basel, 2001. [30] R. Šimkus, R. Baronas, Ž. Ledas. A multi-cellular network of metabolically active E. coli as a weak gel of living Janus particles. Soft Matter, 2013, 9 ( 17 ):4489-4500.

References

[1] R. Baronas , R. Šimkus , Modeling the bacterial self-organization in a circularcontainer along the contact line as detectedby bioluminescence imaging , Nonlinear Anal. Model. Control , 2011 , Vol. 16 , No. 3 , 270 - 282 .

[2] E.O. Budrene , H.C. Berg , Dynamics of formation of symmetrical patterns by chemotactic bacteria , Nature , 1995 , 376 ( 6535 ), pp. 49 - 53 .

[3] M.P. Brenner , L.S. Levitov , E.O. Budrene , Physical mechanisms for chemotactic pattern formation by bacteria, Biophys . J., 1998 , 74 ( 4 ), pp. 1677 - 1693 .

[4] J.D. Murray , Mathematical Biology: II. Spatial Models and Biomedical Applications , Springer, Berlin, 3rd edition, 2003 .

[5] E.F. Keller , L.A. Segel , Model for chemotaxis, J. Theor. Biol ., 1971 , 30 ( 2 ), pp. 225 - 234 .

[6] D. Horstmann , From 1970 until present: the Keller-Segel model in chemotaxis and its consequences, I . Jahresberichte DMV , 2003 , 105 ( 3 ), pp. 103 - 165 .

[7] T. Hillen , K.J. Painter , A user's guide to PDE models for chemotaxis , J. Math. Biol. , 2009 , 58 ( 1-2 ), pp. 183 - 217 .

[8] R. Šimkus , Bioluminescent monitoring of turbulent bioconvection , Luminescence , 2006 , 21 ( 2 ), pp. 77 - 80 .

[9] R. Šimkus , V. Kirejev , R. Meškiene , R. Meškys , Torus generated by Escherichia coli , Exp. Fluids , 2009 , 46 ( 2 ), pp. 365 - 369 .

[10] I.R. Lapidus , R. Schiller , Model for the chemotactic response of a bacterial population, Biophys . J., 1976 , 16 ( 7 ), pp. 779 - 789 .

[11] P.K. Maini , M.R. Myerscough , K.H. Winters , J.D. Murray , Bifurcating spatially heterogeneous solutions in a chemotaxis model for biological pattern generation , Bull. Math. Biol. , 1991 , 53 ( 5 ), pp. 701 - 719 .

[12] K. Kawasaki , A. Mochizuki , M. Matsushita , T. Umeda , N. Shigesada , Modeling spatiotemporal patterns generaged by Bacillus subtilis , J. Theor. Biol ., 1997 , 188 ( 2 ), pp. 177 - 185 .

[13] R. Tyson , S.R. Lubkin , J.D. Murray , Model and analysis of chemotactic bacterial patterns in a liquid medium , J. Math. Biol. , 1999 , 38 ( 4 ), pp. 359 - 375 .

[14] N. Perry , Experimental validation of a critical domain size in reaction-diffusion systems with Escherichia coli populations , J. R. Soc. Interface , 2005 , 2 ( 4 ), pp. 379 - 387 .

[15] M.P. Zorzano , D. Hochberg , M.T. Cuevas , J.M. Gomez-Gomez , Reaction-diffusion model for pattern formation in E. coli swarming colonies with slime , Phys. Rev. E , 2005 , 71 ( 3 ), 31908 .

[16] A. Rabner , E. Martinez, R. Pedhazur , T. Elad , S. Belkin , Y. Shacham, Mathematical modeling of a bioluminescent E. coli based biosensor , Nonlinear Anal. Model. Control , 2009 , 14 ( 4 ), pp. 505 - 529 .

[17] M.D. Egbert , X.E. Barandiaran , E.A. Di Paolo , A minimal model of metabolism-based chemotaxis , PLoS Comput. Biol ., 2010 , 6 ( 12 ), e1001004 .