Nabil Mabrouk
email: nabil.mabrouk@cemagref.fr

Guillaume Deffuant
email: guillaume.deffuant@cemagref.fr

Tim Tolker-Nielsen

Claude Lobry

Bacteria can form interconnected microcolonies when a self-excreted product reduces their surface motility: evidence from individual-based model simulations

Keywords: Biofilm, Pseudomonas aeruginosa, spatial pattern, extracellular DNA

Recent experimental observations of Pseudomonas aeruginosa, a model bacterium in biofilm research, reveal that, under specific growth conditions, bacterial cells form patterns of interconnected microcolonies. In the present work, we use an individual-based model to assess the involvement of bacteria motility and self-produced extracellular substance in the formation of these patterns. In our simulations, the pattern of interconnected microcolonies appears only when bacteria motility is reduced by excreted extracellular macromolecules. Immotile bacteria form isolated microcolonies and constantly motile bacteria form flat biofilms. Based on experimental data and computer simulations we suggest a mechanism that could be responsible for these interconnected microcolonies.

Introduction

Many bacteria have an innate propensity to form biofilms: they build structured multicellular communities attached to solid surfaces. Microscopic examination of biofilms formed by Pseudomonas aeruginosa, a model bacterium in biofilm research [START_REF] Costerton | Microbial biofilms[END_REF][START_REF] Davey | Microbial biofilms: from ecology to molecular genetics[END_REF], reveals a wide diversity of spatial patterns which, depending on the growth conditions, range from a flat thin layer of cells to a patchy pattern with interconnected microcolonies having complex tower or mushroom like shapes (Klausen et al., 2003b[START_REF] Barken | Roles of type IV pili, flagellum mediated motility and extracellular DNA in the formation of mature multicellular structures in Pseudomonas aeruginosa biofilms[END_REF]. At least some explanation for the different biofilm patterns formed under different conditions relates to surface motility. Conditions that promote extensive surface motility can lead to the formation of flat, homogenous biofilms, whereas biofilms characterized by aggregates result from at least a subpopulation of the community ceasing to move at an early stage of biofilm formation (Parsek and Tolker-Nielsen, 2008). In P. aeruginosa biofilms, formation of initial microcolonies (stalks) that subsequently become colonized by cap-forming bacteria is necessary for the formation of mushroom-shaped multicellular structures (Klausen et al., 2003a;2003b). These multicellular structures often provide important benefits such as a higher tolerance to adverse conditions [START_REF] Parsek | Pattern formation in Pseudomonas aeruginosa biofilms[END_REF], and can be crucial in industrial processes.

Therefore, much experimental and theoretical effort is currently devoted to understand their mechanisms of formation.

The formation of spatial patterns in P. aeruginosa biofilms involves a complex interplay between cell proliferation, surface-associated motility and the production of extracellular macromolecules that form a structural matrix (for a review see [START_REF] Parsek | Pattern formation in Pseudomonas aeruginosa biofilms[END_REF]). Evidence has been provided that arrest of type IV pili-mediated motility (twitching motility) plays a role in the formation of the initial microcolonies in P. aeruginosa biofilms, whereas flagelladriven motility (swarming motility) plays a role in the subsequent formation of the cap-portion of the mushroom-shaped structures (Klausen et al., 2003a;2003b; [START_REF] Barken | Roles of type IV pili, flagellum mediated motility and extracellular DNA in the formation of mature multicellular structures in Pseudomonas aeruginosa biofilms[END_REF]. The early stages of biofilm development by P. aeruginosa is dependent on extracellular DNA [START_REF] Whitchurch | Extracellular DNA required for bacterial biofilm formation[END_REF], which is known to bind with high affinity to type IV pili [START_REF] Aas | Competence for natural transformation in Neisseria gonorrhoea: components of DNA binding and uptake linked to type IV pilus expression[END_REF][START_REF] Van Schaik | DNA binding: a novel function of Pseudomonas aeruginosa type IV pili[END_REF], a fimbrae extending from the cell body and mediating surface-associated twitching motility. Hence it may be assumed that initially motile bacteria stop and form microcolonies in the regions of abundant extracellular DNA. Recently, Allesen-Holm et al. ( 2006) visualized the spatial distribution of extracellular DNA and bacteria in P. aeruginosa biofilms, grown in flow chambers on minimal glucose medium. In 2-day-old biofilms the extracellular DNA was present inside the small microcolonies, but accumulated mainly in the outer layer of the microcolonies and between the microcolonies forming a grid-like structure. Zoomed views of the microcolonies revealed that they were often interconnected with thin strands of extracellular DNA covered with bacteria (see Fig. 1). In 4-day-old biofilms caplike multicellular structures had formed on top of the initial small microcolonies, and the highest concentration of extracellular DNA was present between the stalkportion and the cap-portion of the mushroom-shaped multicellular structures. Individual-based modeling is a bottom-up approach as it attempts to model a population or community by describing the actions and properties of the individuals comprising the population or the community [START_REF] Grimm | Ten years of individual-based modelling in ecology: what have we learned and what could we learn in the future?[END_REF][START_REF] Kreft | Individual-based modelling of biofilms[END_REF]. Spatially explicit IBMs, which include a representation of the spatial configuration of the system, allow the simulation of the spatial pattern that emerges from the interactions at the level of the individuals. IBMs are not the only bottom-up simulating approach of studying emergent spatial patterns.

Cellular automaton (CA) modeling is another bottom-up approach used for simulating microcolonies patterns [START_REF] Barker | A cellular automaton model of microbial growth[END_REF] and biofilms (Wimpenny and Colasanty, 1997, Picioreanu et al., 1998a, 1998b[START_REF] Hunt | A computer investigation of chemically mediated detachment in bacterial biofilms[END_REF]. CA models are discrete grid-based models, and differ from IBMs in that they operate on the spatial cells instead of considering the individual settled in them [START_REF] Ferrer | Individual-based Modelling : An Essential Tool for Microbiology[END_REF]. Comparison of IBMs and CA models showed that differences in the simulated biofilm structure may occur for multi-species biofilms due to different approaches in redistributing the surplus biomass (van Loosdrecht et al., 2002) but that both approaches yield similar biofilm structure for a simple mono-species and mono-substrate system [START_REF] Van Loosdrecht | Mathematical modelling of biofilm structures[END_REF]. IBMs are however usually preferred to CA when the individual's variability need to be considered explicitly [START_REF] Kreft | Individual-based modelling of biofilms[END_REF][START_REF] Ferrer | Individual-based Modelling : An Essential Tool for Microbiology[END_REF].

IBMs have been widely used for modeling spatial organization of bacteria within colonies [START_REF] Kreft | Bacsim, a simulator for individual-based modeling of bacterial colony growth[END_REF][START_REF] Ginovart | Individual-based simulations of bacterial growth on agar plates[END_REF] and biofilms [START_REF] Kreft | Individual-based modelling of biofilms[END_REF][START_REF] Xavier | A framework for multidimensional modeling of activity and structure of multispecies biofilms[END_REF] (for a review see [START_REF] Ferrer | Individual-based Modelling : An Essential Tool for Microbiology[END_REF] and [START_REF] Hellweger | A bunch of tiny individuals: individual-based modeling for microbes[END_REF]). Recently, an individual-based model of surface associated populations of P. aeruginosa has been presented [START_REF] Picioreanu | Microbial mobility involvement in biofilm structure formation -a 3d modelling study[END_REF]. The model involves a three-dimensional space and aims to provide a proof-ofprinciple of the implication of motility in the formation of biofilm structure. It reproduces qualitatively the tendency of motile bacteria to form flat biofilms and that of immotile bacteria to form microcolonies by clonal growth, and proposes detachment and reattachment processes of the motile bacteria as possible mechanisms yielding the formation of complex mushroom-shaped microcolonies.

In our model, we focus on the interplay between extracellular DNA production and bacterial motility. We show that a model where bacterial migration is stopped due to adherence to self-produced extracellular DNA can produce complex patterns of interconnected microcolonies.

Methods

We describe the model through the ODD protocol (Overview, Design concepts and Details) [START_REF] Grimm | A standard protocol for describing individual-based and agent-based models[END_REF], which was designed to facilitate descriptions of IBMs in ecology and social science [START_REF] Polhill | Using the ODD protocol for describing three agent-based social simulation models of land-use change[END_REF].

Purpose

The model we propose is intended to qualitatively explore the role of surfacemotility reduction by self-produced macromolecules in biofilm pattern formation.

We address this question at a rather abstract level, and the model results are not compared to specific experimental results.

State variables and scales

The model is a two-dimensional representation of a biofilm system and comprises two entities: bacterial cells and their local environment (patches) (Fig. 2). Bacterial cells are represented as discs. They are characterized by the state variables: continuous position (x,y), individual mass (m), individual diffusion factor (D f ) and individual substrate uptake rate (r)(table 1). A Patch is ∆l characterized by the state variables: discrete position (i,j), substrate concentration (s), extracellular substance concentration (p) and a substrate uptake rate (r s ) corresponding to the sum of uptake rates of the individual cells contained within the patch (table 1).

Scales

We simulated the first day of the biofilm development. We discretized time with a constant time step, denoted ∆t = 1s. A spatial patch has a size of 5 × 5 µm and the whole domain contains 400 × 400 patches (spatial domain side l = 2000 µm).

Bacterial cells have variable diameters (~2 µm) depending on their masses.

Process overview and scheduling

In the IBM, bacteria move stochastically along a 2-D surface while consuming substrate, growing and reproducing by binary fission and excreting a product.

Bacteria movement is then slowed down through interactions with the excreted product yielding different patterns of microcolonies.

At each time step the following processes are performed sequentially: Shoving: bacteria division and motion may produce cells overlaps. In this case cells are displaced using an algorithm proposed by [START_REF] Kreft | Individual-based modelling of biofilms[END_REF] that mimics a shoving process.

Substrate and excreted product mass balances: we modeled the substrate and the excreted product dynamics using two diffusion-reaction equations discretized on the lattice formed by the patches. The reaction terms in these equations are calculated using the substrate uptake rate (r s )

previously calculated for each patch.

Design Concepts

Emergence: the IBM is designed such that the spatial pattern of bacteria and product distribution emerge from local interactions.

Sensing: in our model, a bacterium senses the substrate and product concentration within the patch corresponding to its position. The substrate concentration affects the growth rate of the bacterium whereas the product concentration affects the motility of the bacteria.

Stochasticity: bacteria motility and the positioning of the daughter cells after a division event are the only stochastic processes that we considered in the IBM.

Observation: at each time step the state variables for bacteria and patches are recorded.

Submodels

Bacteria growth

We calculate the individual substrate uptake rate (r(t)) of a cell located in (x, y)

and having a mass m using the following Monod-like kinetic equation:

( ) ( ) ( ) s k t j i s t j i s m t r + ⋅ = , , , , max µ (1)
Where s(i, j, t) is the substrate concentration at patch (i, j) with i = floor(x/∆l) and j = floor(y/∆l) and µ max and k s are Monod kinetic parameters. The growth rate of the bacterium is given by:

( ) t r Y dt dm b ⋅ = (2)
Where Y b is the biomass yield (expressed in mass of bacteria per mass of consumed substrate). The time derivative is discretized using an Euler explicit scheme and the new mass of the cell is calculated by:

( ) ( ) ( ) t r Y t t m t t m b ⋅ ⋅ ∆ + = ∆ + (3)
With ∆t the time step. The product excretion rate r p (t) of the considered individual cell is given by:

( ) ( ) t r Y t r p p ⋅ = (4)

Substrate Uptake rates

For each patch we calculate a substrate uptake rate r s (t,i,j) by summing the individual substrate uptake rates of the cells contained within the patch:

( ) ∑ ∆ = k k s t r l j i t r ) ( 1 , , 2 (5)
where k is the number of cells in patch i,j.

Bacteria division

If the mass of a focal individual (a mother cell) becomes greater than twice the initial mass of an individual(m ≥ 2m 0 ) it divides into two daughter cells each with a mass m/2. The first daughter cell takes the position of the mother cell while the second daughter cell is placed randomly at a distance d (distance between the centers of both cells) corresponding to the diameter of the daughter cells (both daughter cells have the same diameter).

Bacteria motility

The motility of the cells is modeled as a Brownian motion process with an apparent diffusion factor (D f ) which is specific to each individual cell. For a given bacterium located at x, y at time t, the position of the bacterium at the instant t + ∆t is given by:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 , 0 , 2 1 , 0 , 2 N t y x D t y t t y N t y x D t x t t x f f ⋅ ∆ + = ∆ + ⋅ ∆ + = ∆ + (6)
where N(0, 1) draws a number from a centered normal distribution of standard deviation 1 (generated using the Mersenne Twister pseudo-random number generator). As we assumed that bacteria motility was reduced by the excreted product, we calculate diffusion factor (D f ) as a decreasing function of the excreted product concentration in the corresponding patch (i, j). We use the following function:

( )

t j i p D D f f , , 1 1 max ⋅ + ⋅ = β (7)
where D fmax is the maximum diffusion factor of the bacterium, β is a binding affinity factor and p(i, j, t) is the product concentration at patch (i,j). The parameter β rules the sensitivity of D f to the variation of p(i, j, t) as shown by We use periodic boundary conditions for the individuals such that an individual that exits the domain from one side appears from the opposite side.

Bacteria shoving

Bacteria shove each other when they overlap. If a bacterium with radius a is overlapped with n neighboring cells, it is displaced with a shoving vector d calculated using the following equation adapted from [START_REF] Kreft | Individual-based modelling of biofilms[END_REF]:

k n k k k u d a a d ∑ = - + = : 1 2 ( 8
)
a k is the radius of the neighbor cell k, d k is the Euclidean distance from the center of the bacterium to the center the k th neighboring cell and u k is a vector directed from the center of neighbor bacterium k towards the center of the bacterium and having a unitary norm. We consider periodic boundary conditions for the search of overlapping neighbors and for the repositioning of the shoved individuals.

Substrate and excreted product mass balance equation

The distribution of the substrate is the solution of the following continuous diffusion-reaction equation:

s s r s D t s - ∇ = ∂ ∂ 2 (9)
with periodic boundary conditions: [START_REF] Ferrer | Individual-based Modelling : An Essential Tool for Microbiology[END_REF] and having as initial conditions:

( ) ( ) y l x t s y x t s , , , 0 , = = = (10) ( ) ) , , ( 0 ,
( ) 0 , , 0 s y x t s = = (12)
The excreted product dynamic is also given by a diffusion-reaction equation: [START_REF] Grimm | Ten years of individual-based modelling in ecology: what have we learned and what could we learn in the future?[END_REF] with periodic boundary conditions:

s p p r Y p D t p + ∇ = ∂ ∂ 2
( ) ( ) y l x t p y x t p , , , 0 , = = = (14) ( ) ( ) l y x t p y x t p = = = , , 0 , , (15)
and having as initial conditions

( ) 0 , , 0 = = y x t p (16)
Y p in (eq. 13) is the product yield expressed in mass of excreted product per mass of consumed substrate. We discretize the substrate and the product mass balance equations with respect to space on the lattice formed with the patches using a four-point scheme. The reaction term has already been calculated for each patch (eq. 5). We discretize the time derivative term in the substrate and the product mass balance equations using an implicit scheme for the diffusion term. Note that the reaction term is calculated on the basis of an explicit scheme. The obtained discretized system is a sparse linear system that we solve using the conjugate gradient method [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF]. This gives the new substrate and product concentrations in each patch.

IBM parameters

Unless explicitly specified, we use the parameters values in table 2 for the individual-based model.

Initialization

We initialize all simulations with: a uniform initial substrate concentration s 0 a uniformly null concentration of excreted product N 0 bacterial cells drawn at random in the domain each with the same initial mass m 0 .

Initial conditions are detailed in table 3.

Model implementation

We implemented the model using the Java programming language and the Mason framework. Mason is a discrete event multiagent simulation library code developed at the George Mason University for implementing multi-agent models [START_REF] Luke | Mason : A Java multi-agent simulation toolkit[END_REF].

Results

We simulate the patterns formed by the spatial distribution of bacteria and the excreted product displayed after 9 and 18 hours for different values of β (the parameter determining the impact of the excreted product on bacteria motility) and for the case of immotile bacteria. Two extreme cases can be identified. In the first case, bacteria are motile and their motility is independent from the excreted product (β = 0) while in the second case bacteria are immotile D fmax = 0. Figures 4 and5 show the patterns obtained in these limiting cases. Motile cells disperse over the spatial domain (Fig. 4) in contrast with immotile cells that form isolated microcolonies (Fig. 5). These results are consistent with previous theoretical [START_REF] Picioreanu | Microbial mobility involvement in biofilm structure formation -a 3d modelling study[END_REF] and experimental (Klausen et al., 2003b) studies of P.

aeruginosa showing that motile cells tend to form flat biofilms while immotile bacteria form round shaped microcolonies. We also simulated intermediate cases where the bacteria motility depends on the local concentration of the excreted product. Figure 6 and7 show examples of spatial pattern obtained after 9 and 18 hours for a small and a large value of the parameter β respectively. As the system evolves over time, the patterns of microcolony formation under the two parameterizations begin to diverge. A large value of β results in microcolonies that are more spatially discrete with higher densities of individuals within them (Figure 7), while the smaller value of β results in a pattern of microcolonies that are more amorphous in size and shape, are more connected with each other, and densities of bacteria within them are lower (Figure 6).

In our simulations, microcolony formation is initiated by the local accumulation of the product excreted by the cells along their Brownian trajectories. The product excretion rate is maximal (~Y p µ max , see Equations 1 and 6) at the beginning of the simulation when the substrate (S>> k s in Equation 1) is abundant and accumulates due to its low diffusion factor. For large values of β bacteria are rapidly entrapped within the locations containing the excreted product and their daughter cells tend to accumulate locally yielding dense and discrete microcolonies. In the opposite for small values of β, bacteria and their daughter cells tend to disperse and the yielded microcolonies have amorphous shape. Analogously, simulation with different values of the product excretion ratio Y p (Equation 6), at a constant value of β, between the two extreme cases of Y p =0 (no product excretion) and Y p =1.0 (no growth, all the substrate is released back on the form of product) yields patterns that vary respectively from uniform distribution of the bacteria to the formation of isolated, round-shaped microcolonies (data not shown). This suggests that the binding affinity factor (β) and the product excretion rate may have a significant impact on the patterns of spatial distribution of the bacteria. In our model, the substrate concentration impacts the of product excretion. An initial low substrate level yields low rates of product excretion of the individuals which may not be sufficient to reduce the motility of the bacteria and the formation of microcolonies. In our simulation the initial substrate level is relatively high and only 90% of this initial stock is consumed by the end of the simulation. We investigate how the interconnections between the microcolonies form in the case of large value of binding affinity parameter (β=500). Figure 8 shows a zoomed view of the formation of the interconnection between two neighboring microcolonies. The interconnection seems to be created by bacteria which go from a microcolony to a neighboring one, and which progressively accumulate excreted product on the path. The resultant pattern is formed with dense and discrete microcolonies interconnected with relatively thin strands of bacteria and excreted product.

Discussion

In Individual-based models are appealing to microbiologists because of the emphasis on the individual cell as the fundamental unit and the richness of their dynamic.

However when they are too complex they often become difficult to analyze which limits their practical payoff [START_REF] Grimm | Individual-based Modeling and Ecology[END_REF]. Therefore, we have attempted to keep our model as simple as possible and inspired by the phenomenological results mentioned in the introduction, but rich enough to produce patterns of interconnected microcolonies. [START_REF] Johnson | Microcolony and biofilm formation as a survival strategy for bacteria[END_REF] proposed a model based on the assumption of the "direct" interaction between individual cells through attractive (and repulsive) forces arguing that forces between cells can be a proxy for the behavior we expect to see due to chemotaxis in response to chemicals released by other cells [START_REF] Johnson | Microcolony and biofilm formation as a survival strategy for bacteria[END_REF]. The assumption simplifies the model as the dynamic of the excreted product is not considered explicitly. This is a reasonable assumption when the dynamics of the chemical being produced is fast (high diffusion) compared to the dynamic of the bacteria (motility and growth) [START_REF] Lee | Non-local concepts and models in biology[END_REF]. Our model is based on the assumption that the excreted macromolecules like DNA and exopolymers diffuse at low rates and persist on the path of the motile bacteria. This seems to play a role in the formation of the interconnected microcolonies as microcolonies yielded by Johnson's model through attractive/repulsive forces seem not to be interconnected [START_REF] Johnson | Microcolony and biofilm formation as a survival strategy for bacteria[END_REF].

Our study included, however, some limitations that should be acknowledged for future research in this area. Examples, inherent to individual-based modeling approach, include the question of how to measure the "goodness" of an IBM. [START_REF] Grimm | Individual-based Modeling and Ecology[END_REF] suggested that testing an IBM against multiple observed patterns is a powerful way to assess the IBM capacity to capture system's essential characteristics. In this work we assessed the "goodness" of our IBM through comparison of an observed pattern (connected microcolonies in fig. 1) to a simulated one (fig. 8). Despite the similarities between the experimental and simulated patterns (thin path of bacteria and product connecting neighboring microcolonies) there are still differences. For example cells surrounded the thin path in figure 1 while they are contained in the area with high amount of product in figure 8. Therefore, additional patterns should be identified and compared to the simulated ones in order to strengthen or falsify the hypothetical implication of self-excreted product in motility reduction and formation of interconnected microcolonies. Another limitation of our model is the determination of parameters like the binding affinity factor (β) and cell diffusion (D f ). These parameters may be difficult to obtain from experimental observations and their impact should be assessed through sensitivity analysis.

Finally, our work can be extended beyond its basic premise and valuable research remains to be conducted. For instance the model can be extended to investigate the effect of other macromolecules like exopolymers on the bacteria motility and microcolonies formation. Another potential extension is to include additional processes like bacteria detachment and substrate feeding which allow investigating the effect of the interactions between cell motility and excreted macromolecules on the architecture of mature biofilms.

Fig. 1

 1 Fig. 1 Horizontal confocal laser scanning microscope section acquired in a 2-dayold DDAO-stained biofilm formed by Gfp-tagged P. aeruginosa PAO1. The images show the green fluorescent bacteria (A), the red fluorescent extracellular DNA (B), and an overlay of the two (C). Reproduced from Molecular Microbiology 59: 1114-1128 with permission from John Wiley & Sons, Inc.

Fig. 2

 2 Fig. 2 Schematic representation of the model entities: bacterial cells (discs) and patches (squares with a side ∆l). A bacterium with continuous coordinates x,y is contained in the patch i=floor(x/∆l,) j=floor(y/∆l).

  Bacteria growth: for each individual cell we (1) calculate and individual substrate uptake rate (r) which depends on the mass of the cell and the local substrate concentration.[START_REF] Allesen-Holm | A characterization of DNA release in Pseudomonas aeruginosa cultures and biofilms[END_REF] We use the calculated uptake rate to update the mass of the cell.Substrate uptake rates for patches: given the uptake rates of each individual cell, we calculate an uptake rate (r s ) for each patch by taking the sum of the uptake rates of the individual cells contained within it.Division: for each individual cell we compare the cell mass to a critical value. If the cell mass is higher than the critical value than the cell is divided into two daughter cells. One of the daughter cells takes the position of the mother cell while the second is placed at random around the mother cell position at a distance (d) corresponding to the diameter of the daughter cell.Surface motility: we model bacteria translocation as a Brownian process using an diffusion factor (D f ) proper to each individual. For each individual cell, including newly formed cells, we calculate a diffusion factor (D f ) using a decreasing function of the local product concentration and then move the individual cell accordingly.

Fig. 3

 3 figure 3.

Fig. 4

 4 Fig. 4 Overall and zoomed view of the spatial pattern formed by the bacteria (black dots) and the excreted macromolecules (gray scale) for motile bacteria for β= 0 at time 9 and 18 hours

Fig. 5 Fig. 6 Fig. 7

 567 Fig. 5 Overall and zoomed view of the spatial pattern formed by the bacteria (black dots) and the excreted macromolecules (gray scale) for immotile bacteria at time 9 and 18 hours

Fig. 8

 8 Fig. 8 Zoomed view of interconnection formation between two microcolonies for a binding affinity β = 500 at time 9 and 18 hours. Bacteria are represented with green dots and the excreted product concentration with gray scale.

  this work we abstract experimental observations on P. aeruginosa biofilm development (Whitchurch et al., 2002, Klausen et al., 2003a; 2003b; Allesen-Holm et al., 2006; Barken et al., 2008) in an individual-based model and investigate how bacteria motility reduction due to a self-produced substance yields different spatial patterns during the early stages of biofilm development.Our simulation results suggest that self-produced substance-mediated motility reduction does play a role in microcolony formation. Furthermore, in some cases, these simulated microcolonies build interconnections, similar to interconnected microcolonies observed during the early stages of P. aeruginosa biofilm development in flow chambers[START_REF] Allesen-Holm | A characterization of DNA release in Pseudomonas aeruginosa cultures and biofilms[END_REF], and also in biofilms formed by Pseudomonas species in marine environments[START_REF] Dalton | Substratum-induced morphological changes in a marine bacterium and their relevance to biofilm structure[END_REF][START_REF] Dalton | Diversity in surface colonization behavior in marine bacteria[END_REF].Several authors already explored the mechanisms yielding microcolonies in biofilms[START_REF] Alpkvist | Three-dimensional biofilm model with individual cells and continuum EPS matrix[END_REF][START_REF] Picioreanu | Microbial mobility involvement in biofilm structure formation -a 3d modelling study[END_REF]. However, the pattern of interconnected microcolonies cannot be obtained with these usual mechanisms: immotile bacteria form isolated microcolonies and constantly motile bacteria form flat biofilms. Based on experimental data and computer simulations we suggest a mechanism that could be responsible for the observed patterns. Our model shows that microcolonies may result from bacteria motility reduction by self-produced macromolecules. The analysis of the simulation results suggests that cells on the edge of a microcolony occasionally detach and undergo a surface-associated motility until being captured by a neighboring microcolony. The path of the migrating cell is marked by the excreted macromolecules and is progressively reinforced by other migrating cells. This results in the formation of an interconnection between the neighboring microcolonies. However, more investigation is necessary to strengthen or falsify this hypothesis. This investigation could be important for a better understanding of biofilm functions. Indeed, it is well accepted that the presence of different subpopulations in microcolonies can favor the survival of one or more subpopulations under adverse conditions[START_REF] Parsek | Pattern formation in Pseudomonas aeruginosa biofilms[END_REF] Tolker-Nielsen, 2008, Johnson, 2008). However, recent observations suggest that in some cases the matrix of extracellular DNA (and possibly other products) may guide migration of the cells between microcolonies[START_REF] Lu | Exopolysaccharide biosynthesis genes required for social motility in Myxococcus xanthus[END_REF][START_REF] Barken | Roles of type IV pili, flagellum mediated motility and extracellular DNA in the formation of mature multicellular structures in Pseudomonas aeruginosa biofilms[END_REF] favoring a type of collaboration between distinct subpopulations within a biofilm. This may be related to studies in landscape ecology where this pattern of interconnected habitats is recognized as particularly important(Burel and Baudry, 2003).

Table 1

 1 State variables of the individual-based model

		Variables	Description
	Bacterium state	x, y	Continuous position of the center of the bacterium
		m	Mass
		D f	Diffusion factor
		r	Individual substrate uptake rate
	Patch state	i,j	Discrete position
		s	Substrate concentration
		p	Excreted product concentration
			Reaction rate (Sum of
		r s	substrate uptake rates of the individual cells contained
			within the patch)

Table 2

 2

	Individual-based model parameters. (*): assumed, (**): adapted from [30], (***) adapted
	from [22]
	Parameter Description	Units	Value
	l	Domain side	µm	2000 (*)
	∆l	Spatial step	µm	5 (*)
	D s	Substrate diffusion constant m 2 /s	10 -10 (**)
	D p	Product diffusion constant	m 2 /s	10 -16 (*)
	D fmax	Maximum diffusion factor for the bacteria	m 2 /s	10 -12 (*)
	µ max	Maximum growth rate	kg substrate /(kg biomass s) 10 -4 (**)
	k s	Affinity constant	kg substrate /m 2	0.01 (**)
	Y p	Product yield	(kg product /kg substrate )	0.1 (*)
	Y b	Biomass yield	(kg biomass /kg substrate )	0.9 (*)
	ρ b	Density of a bacterial cell	(kg biomass /m 2 )	290 10 -6 (***)
	β	Binding affinity factor	(m 2 /kg product )	variable

Table 3

 3 Initial conditions of the individual-based model

	Parameter Description	Units	Value
	s 0	Initial substrate concentration in all patches	kg/m 2	10.0
	p 0	Initial product concentration in all patches	kg/m 2	0.0
	N 0	Initial number of bacterial cells	#	100
	m 0	Initial mass of a bacterium	kg	10 -15

Acknowledgements

The research was partly supported by Patres Project (NEST 043264).


Guillemin Jean-Philippe
email: jp.guillemin@giat-industries.fr

Brunet Luc

Bonnefoy Olivier

Thomas Gérard

A Flow Time Model for Melt-Cast Insensitive Explosive Process

Keywords: Melt-Cast Insensitive Explosive, Rheology of Concentrated Suspensions, Flow Time Model, Maximum Packing Density, Quemada's Model for Viscosity, De Larrard's Model

These parameters will be used as a basis to establish our model. First, a relationship between the viscous and flow time terms starting from Navier-Stokes equation is needed. Secondly, evaluations of some characteristics of the granular phase in concentrated suspensions will be obtained and inserted in our flow time equations. De Larrard's model for the calculation of the maximum packing density f m [4] will be considered. The study of the dependence of the casting time upon weak variations of the formulation and upon the choice of raw materials will make possible to compare the experimental results and those given by our model. Melt casting energetic formulations are established for use in a cylindrical tank of a total volume of 5700 cm 3 . The experiments for flow time determinations are all performed in this container.

Since the publication of Einstein analysing the viscosity of dilute suspensions of rigid spheres in a viscous liquid, numerous equations have been developed to try and extend Einstein's formula to suspensions of higher concentrations [5]. The various resulting formulas differ considerably from each other. In this study, three relationships have been selected to express the viscosity of dispersions of spherical particles as a function of the viscosity of the interstitial fluid , the volume fraction of solids f and the maximum packing density of the solids f 0 m m . Others equations are given in the literature, but the three retained hereafter seem to be the most encountered, discussed and reliable, and this is the reason why they have been chosen in this study. These equations are given in Tab. 1.

Introduction

The purpose of this article is to propose a predictive model of the flow time necessary for emptying a reactor filled with a concentrated suspension through a nozzle situated at the bottom. Similarly to our industrial process, the suspension is made out of insensitive energetic materials and flows under gravity. Experimental results are compared with three viscosity models largely used in the field of concentrated suspensions (Quemada [START_REF] Quemada | Rheology of Concentratred Disperse Systems and Minimum Energy Dissipation Principle, I : Viscosity-Concentration Relationship[END_REF], Krieger and Dougherty [START_REF] Krieger | A mechanism for non-Newtonien flow in suspensions of rigid spheres[END_REF] and Mooney [START_REF] Mooney | The Viscosity of a Concentrated Suspension of Spherical Particles[END_REF]). These models give the viscosity of a concentrated suspension as a function of the viscosity m of the pure liquid, the volume fraction f of the solid included in the paste, and the maximum packing density f 0 m .

Table 1. Rheological equations

Authors

Rheological equations

Krieger-Dougherty [2]

Quemada [START_REF] Quemada | Rheology of Concentratred Disperse Systems and Minimum Energy Dissipation Principle, I : Viscosity-Concentration Relationship[END_REF] Mooney [START_REF] Mooney | The Viscosity of a Concentrated Suspension of Spherical Particles[END_REF] 2.5

0 1 m m f f m m f - ae ö = - ç ÷ è ø 2 0 1 m f m m f - ae ö = - ç ÷ è ø 0 2.5 exp 1 m f m m f f ae ö ç ÷ ç ÷ = ç ÷ - ç ÷ è ø

Relationship between Viscosity and Flow Time

The Bernoulli relationship can be written: Because of the high solid content and the low flow velocity, it is assumed here that Reynolds number is not too important, although greater than unity. The inertial term of the Navier-Stokes equation is supposed to be small enough with respect to the viscous term. The momentum equation takes the following simplified form: The mass balance gives:

2 2 1 1 2 2 r r = v d v d (5)
With v d , so that Eq. ( 4) could be written now:

1 / = -h d . . v p v t r m ¶ = -Ñ + D + ¶ r r r r . g r (1)
.( )

+ ae ö = ç ÷ è ø dh a h b dt c (6)
As the suspension flows during a relatively long time (over one minute), the state can be regarded as quasi-stationary and one can write at any time: Where .

m r D = Ñ - r r r v p g ( 2
)
The integration of Eq. (6) between t=0 and , where the liquid heights are respectively and leads to Eq.7:

t = t 0 h ( ) t h
Eq. (2) implies that, for a given pressure gradient, the flow rate of a fluid is inversely proportional to its viscosity. This can be applied to the two fluids considered, the suspension with a viscosity and the interstitial fluid with a viscosity Therefore, for a given volume, the flow time t is directly proportional to the viscosity:

m 0 m 0 ( ) 4. t t + = + - a h b h b c (7)
Att , the tank is empty and . Finally, is given by Eq. ( 8)

0 0 0 ( ) 0 t = h t 0 0 t m t m = (3) 4 1 1 2 0 0 2 ( ) ( ) 2 1 . . t r r 1 2 é ù ae ö ae ö - - = - + - ê ú ç ÷ ç ÷ ê ú è ø è ø ë û d PP h d g g PP g (8)
Wheret represents the flow time for the interstitial fluid. 0 With Eq. ( 3), Eq. ( 8) and one of the rheological equations given in Tab. 1, three relationships, Eq. (9), Eq. (10) and Eq. ( 11), allow us to calculate the flow time for our energetic concentrated suspensions. These relationships depend on the solid volume fraction f , the flow for the interstitial fluidt and the maximum packing densityf . 0 m

Calculation of t by Bernoulli's approach

The interstitial fluid is a mixture of trinitrotoluene TNT and additives, presenting viscosities of 11 and 48.6 mPa.s respectively. Such values are rather small, and the interstitial fluid can be assimilated to a perfect fluid: this allows us to use the Bernoulli relationship in the tank described on Fig. 1.

2.5 0 1 m m f f t t f - ae ö = - ç ÷ è ø (9) 2 0 1 m f t t f - ae ö = - ç ÷ è ø (10) 0 2.5 exp 1 m f t t f f ae ö ç ÷ ç = ç ÷ - ç ÷ è ø ÷ (11)

Calculation of the Maximum Packing Density

Many studies aim at improving predictions of viscosity by refining the methods of calculation of the maximum packing density [6-7-8]. Many strong similarities appear between concrete and explosive made by casting process as the great number of granulometric scales, and the range of the component morphologies. In this kind of approach, the model developed by De Larrard seems very attractive [9]. This model predicts the maximum packing density of a polydisperse mix, from three parameters: the particle size distribution of the mix, its true density and the experimental packing density of the solid species. The use of software is required to determine the maximum packing density. In this study, the RENE-LCPC software developed by De Larrard and T. Sedran [10] has been selected. The details of the algorithm will not be described here but the interested reader may refer to some of the associated publications [11-12-13].

This model deals with grain mixtures in which linear combinations of packing densities allow to predict the packing density of a mixture of monosized particles d i (d 1 <d 2 <…<d n ) from:

1 1 1 1 (1 ) b f b a a - = = + = -- - å å i i i N i i j i i j j j i a b j f i (12) 1 min f < < = m i N ( 1 3 )
3 Raw Materials and Experimental Set-up

Raw Materials

The compounds used in this study are listed in Tab. 2.

Rheology of the Liquid Phase

The liquid phase is composed of TNT and a fusible additive A. The Newtonian viscosity of the additive is measured by a viscometer presenting a"Couette" geometry (Rheomat 30, Contrave). A thermo-regulated bath is used to control the temperature. Results show that the additive has a Newtonian behaviour and its Newtonian viscosity is equal to 48.6 mPa.s at 85°C. The TNT Newtonian viscosity is calculated with Eq. (14) established by [14] on temperature as follows:

m TNT * exp m m ae = ç è ø TNT Q T ö ÷ ( 14
)
Where T is the absolute temperature in Kelvin and and .

* 7

5.41 10 .

m - = × Pa s 3570 = Q K
With these parameter values, at 85°C.

11

. m = TNT mPa s

Solid Phase Characterization

The solid phase is made out of two species: NTO and aluminium. Four batches of NTO, labelled NTO 1, NTO 2, NTO 3, NTO 4 and two batches of aluminium, Al 1, Al 2 are used. They differ by their packing densities, grain size distributions and morphologies.

Where f i is the packing density of class; a i is the volume fraction of belonging to class i; b i is the residual packing density i.e. when the class i is alone and fully packed. To compute the packing density of the overall mixture, one considers that the bulk volume of the class i fills the porous space around the coarser grains; moreover, the volume of finer classes inserted in the voids of class i must be added. Two interaction effects must be taken into account in this calculation: the wall effect, a ij , exerted by the coarser particles and the loosening effect, b ij , exerted by the finer particles. Finally, f m is given by Eq. ( 13).

The experimental packing density C is calculated from the true and bulk densities. The true density is measured with a helium pycnometer from Micromeritics and the true bulk density with a volumenometer. C is defined by Eq. (15):

r t r b r r = t b C ( 1 5 )
The particle size distribution is determined by sifting under vibrations NTO particles, and for aluminium particles with a laser granulometer (Malvern Mastersizer 2000). Grain morphologies are determined by scanning electron microscopy. Average sizes of NTO particles reach 350-400 µm and 13 µm for aluminium particles. SEM analysis shows spherical morphology for Al 2.

Experimental Set-Up

Nine insensitive explosive formulations are elaborated with different volume fractions of raw materials. These formulations are given in the left hand side of Tab. 4. The flow time of 5700 cm 3 of each suspension has been measured and the texture of these energetic pastes has been observed. For the energetic formulations exhibiting a high apparent viscosity, the measurements of the flow time has been made several times to get significant results. The measurement precision decreases when explosive formulations are highly concentrated. In this case, we repeated twice the measurement of the flow time.

As in our industrial process, suspensions flow through a nozzle situated at the bottom of a tank under gravity. The dimensions of the cylindrical tanks are represented on Fig. 4. The flow time and texture of each formulation is reported in the right hand side of Tab. 4 and depicted on Fig. 5. The last formulation (Formulation 9) represents the interstitial fluid, composed of TNT and additive only.

Comparison between Experimental and Modelling Results

The maximum packing densityf and the ratio between the maximum packing density and the solid volume fraction f are given in the left hand side of Tab. 5. m Theoretical flow times are calculated Eq. (9), Eq. (10) and Eq. (11). Final results are listed in the right hand side of Tab. 5. Results from Eq. ( 10), which are a direct application of Quemeda's model, are in agreement with our experimental results.

Figure 6. Comparison between theoretical and experimental results

Discussion

We assumed that the interstitial fluid can be considered as inviscid; this seems to be quite realistic since the experimental flow time (11 s) for 5700 cm 3 corresponds to the theoretical time (13 s) calculated by Eq. (8). Moreover, the fair agreement between our experimental results for concentrated suspensions and the theoretical predictions from Quemada's model seems to support our second, intuitive, assumption that the loss of energy due to viscous effects remains negligible. The ratio of the maximum packing density to the solid volume fraction (left hand side of Tab 5) gives an interesting information to be compared with the texture of insensitive energetic pastes (Fig 5). When it is close to 1 (Formulations 3 and 4), i.e. near the jamming of the structure, the texture is very pasty. On the opposite side, when f /f m decreases (Formulations 1, 2 and 5) or far from 1, energetic pastes present a more liquid behaviour.

As shown on Fig. 6, theoretical values of the flow time resulting from Eq. (9) and Eq. (11) show significant deviations from experimental results. We find that Eq. (10) derived from Quemada's model reflects our experimental results more nicely.

Overall, this result completely ignores the non-Newtonian nature of these energetic suspensions. Typically, energetic suspensions exhibit yield stresses, wall slip and sensitivity to shearing stresses. For more details of the viscoplastic behaviour of concentrated energetic materials, the reader will be able to refer to last work of Kalyon [15-16].

Conclusions

For suspensions of different formulations, various experiments of casting under gravity have been carried out and the time required by a given volume to pass through a nozzle has been measured.

Considering that the suspension flow proceeds under conditions where the inertial term of the Navier-Stokes equation is negligible compared to the viscous term, we find the suspension flow time to be directly proportional to its viscosity. Consequently, the cast time can be calculated from the suspension viscosity. Three models, frequently used in rheology, make possible to calculate this viscosity as a function of the interstitial fluid Newtonian viscosity and the ratio f/f m where f m is calculated with De Larrard's model. The relatively low viscosity of the interstitial fluid allows us to consider perfect fluids and apply the Bernoulli relationship. Comparison between theoretical values and experimental results shows that Quemada's model gives a satisfactory modelling. In the configuration studied, the flow time can be estimated by the following relationship:

Figure 1 .

 1 Figure 1. Tank geometry Where V [m 3 ] is the suspension volume; and [m/s] are the fluid velocities; and [Pa] are the top and bottom pressures; and [m] are the diameters of sections 1 and 2; [m] is the fluid height; g [m/s²] is the acceleration of gravity equal to 9.81 m/s².1 v

Figure 2 .Figure 3 .

 23 Figure 2. Al and NTO particle size distributions

Figure 4 .

 4 Figure 4. Left: Industrial process; Right, scheme of the tank with some key dimensions

Table 2 .

 2 Studied compounds

	Compounds	Name	True Density Melting Point (g/cm 3 ) (°C)	Physical state at exp. Conditions
	TNT	2,4,6-trinitrotoluene	1.65	81	liquid
	A	Additive	1	83	liquid
	Al	Aluminium	2.7	660	solid
	NTO 3-nitro-1,2,4-triazol-5-one 1.92	279	solid

Table 3 .

 3 NTO and Al experimental packing densities

	Batches	r t (g/cm	3 )	r b (g/cm	3 )	C
	NTO 1	1.92		0.99		0.52
	NTO 2	1.92		0.97		0.51
	NTO 3	1.92		0.75		0.39
	NTO 4	1.92		0.79		0.41
	Al 1	2.7		1.1		0.41
	Al 2	2.7		1.2		0.44

Table 4 .

 4 Composition and flow time of each insensitive energetic formulation

	Formulation	Al	Batche	NTO	TNT	Volume Fraction f i A	Al	NTO	t (s)
	1	1		2	0.33	0.13	0.09	0.45	72
	2	1		1	0.33	0.13	0.09	0.45	168
	3	1		3	0.33	0.13	0.09	0.45	720
	4	1		3	0.33	0.13	0.09	0.45	1200
	5	2		2	0.33	0.13	0.09	0.45	64
	6	1		4	0.35	0.14	0.07	0.44	330
	7	1		4	0.35	0.14	0.07	0.44	303
	8	2		4	0.35	0.14	0.07	0.44	220
	9	N / A		N / A	0 . 7 2	0 . 2 8		0	0	1 1
	Figure 5. Texture of insensitive energetic pastes

Table 5 .

 5 Maximum packing density of each formulation and comparison between theoretical (ttheo) and experimental results (texp)

	Formulation f m f / f m	Quémada Krieger-Dougherty Mooney t theo (s)	t exp (s)
		1	0,79 0,68		128		124		912	72
		2	0,79 0,68		128		124		912	168
		3	0,61 0,89		972		348		1645790	720
		4	0,61 0,89		972		348		1645790	1200
		5	0,79 0,68		121		121		815	64
		6	0,63 0,81		353		174		10334	330
		7	0,63 0,81		353		174		10334	303
		8	0,63 0,81		353		174		10334	220
		9	0	0		1 3		1 3		1 3	1 1
				R²		0,942		0,914		0,910
		2000
								Eq. (9)
		1800						Eq. (10) Eq. (11)
								Experimental results
		1600
		1400
		1200
	t [ s ]	1000
		800
		600
		400
		200
		0.6 0	0.65	0.7	0.75	0.8	0.85	0.9	0.95	1
							f / f m

Acknowledgments

The authors thank L. Chouris as the Person in Charge of the realization of the tests.

Bulk density (g/cm 3 )