A biophysical threshold for biofilm formation

  1. Jenna A Moore-Ott
  2. Selena Chiu
  3. Daniel B Amchin
  4. Tapomoy Bhattacharjee
  5. Sujit S Datta  Is a corresponding author
  1. Department of Chemical and Biological Engineering, Princeton University, United States
  2. Andlinger Center for Energy and the Environment, Princeton University, United States

Abstract

Bacteria are ubiquitous in our daily lives, either as motile planktonic cells or as immobilized surface-attached biofilms. These different phenotypic states play key roles in agriculture, environment, industry, and medicine; hence, it is critically important to be able to predict the conditions under which bacteria transition from one state to the other. Unfortunately, these transitions depend on a dizzyingly complex array of factors that are determined by the intrinsic properties of the individual cells as well as those of their surrounding environments, and are thus challenging to describe. To address this issue, here, we develop a generally-applicable biophysical model of the interplay between motility-mediated dispersal and biofilm formation under positive quorum sensing control. Using this model, we establish a universal rule predicting how the onset and extent of biofilm formation depend collectively on cell concentration and motility, nutrient diffusion and consumption, chemotactic sensing, and autoinducer production. Our work thus provides a key step toward quantitatively predicting and controlling biofilm formation in diverse and complex settings.

Editor's evaluation

In this work, the authors develop a continuum description of biofilm formation from initially planktonic cells. The coupled partial differential equations that encode the dynamics of the cell populations, nutrients and autoinducers contain many parameters, but it is shown that only two dimensionless combinations of them are needed to understand the threshold for biofilm formation. This work should be of broad interest to a wide range of researchers in biophysics and cell biology.

https://doi.org/10.7554/eLife.76380.sa0

Introduction

Dating back to their discovery by van Leeuwenhoek over three centuries ago, it has been known that bacteria typically exist in one of two phenotypic states: either as motile, planktonic cells that self-propel using e.g., flagella or pili (“animalcules … moving among one another”; Van Leewenhoeck, 1677), or as immobilized, surface-attached biofilms (“little white matter … in the scurf of the teeth”; Leewenhoeck, 1684). These different states have critical functional implications for processes in agriculture, environment, industry, and medicine. For example, motility-mediated dispersal of planktonic cells enables populations to escape from harmful conditions and colonize new terrain (Adler, 1966a; Adler, 1966b; Saragosti et al., 2011; Fu et al., 2018; Cremer et al., 2019; Bhattacharjee et al., 2021)—underlying infection progression, drug delivery to hard-to-reach spots in the body, food spoilage, interactions with plant roots in agriculture, and bioremediation of environmental contaminants (Balzan et al., 2007; Chaban et al., 2015; Datta et al., 2016; Harman et al., 2012; Ribet and Cossart, 2015; Siitonen and Nurminen, 1992; Lux et al., 2001; O’Neil and Marquis, 2006; Gill and Penney, 1977; Shirai et al., 2017; Thornlow et al., 2015; Toley and Forbes, 2012; Dechesne et al., 2010; Souza et al., 2015; Turnbull et al., 2001; Watt et al., 2006; Babalola, 2010; Adadevoh et al., 2016; Adadevoh et al., 2018; Ford and Harvey, 2007; Wang et al., 2008; Reddy and Ford, 1996; Martínez-Calvo et al., 2021). In addition, the formation of immobilized biofilms can initiate antibiotic-resistant infections, foul biomedical devices and industrial equipment, or conversely, help sequester and remove contaminants in dirty water (Davey and O’toole, 2000; Hall-Stoodley et al., 2004; Mah et al., 2003; O’Toole and Stewart, 2005; Fux et al., 2005; Nicolella et al., 2000; Donlan and Costerton, 2002; Davies et al., 1998). Hence, extensive research has focused on understanding bacterial behavior in either the planktonic or biofilm state.

For example, studies of planktonic cells have provided important insights into bacterial motility—which can be either undirected (Berg, 2018; Berg, 2004; Bhattacharjee and Datta, 2019a; Bhattacharjee and Datta, 2019b) or directed in response to e.g., a chemical gradient via chemotaxis (Adler, 1966b; Adler, 1966a; Saragosti et al., 2011; Fu et al., 2018; Cremer et al., 2019; Bhattacharjee et al., 2021; Keller and Segel, 1971; Odell and Keller, 1976; Keller and Odell, 1975; Lauffenburger, 1991; Seyrich et al., 2019; Croze et al., 2011; Amchin et al., 2022). These processes are now known to be regulated not just by intrinsic cellular properties, such as swimming kinematics and the amplitude and frequency of cell body reorientations, but also by the properties of their environment, such as cellular concentration, chemical/nutrient conditions, and confinement by surrounding obstacles (Berg, 2018; Berg, 2004; Bhattacharjee and Datta, 2019a; Bhattacharjee and Datta, 2019b; Adler, 1966b; Adler, 1966a; Saragosti et al., 2011; Fu et al., 2018; Cremer et al., 2019; Bhattacharjee et al., 2021; Keller and Segel, 1971; Odell and Keller, 1976; Keller and Odell, 1975; Lauffenburger, 1991; Seyrich et al., 2019; Croze et al., 2011; Amchin et al., 2022). Thus, the manner in which planktonic bacteria disperse can strongly vary between different species and environmental conditions.

Similarly, studies of biofilms under defined laboratory conditions have also provided key insights—such as by revealing the pivotal role of intercellular chemical signaling in biofilm formation (Nadell et al., 2008; Bassler and Losick, 2006; Davey and O’toole, 2000; Hall-Stoodley et al., 2004). In this process, termed quorum sensing, individual cells produce, secrete, and sense freely diffusible autoinducer molecules, thereby enabling different bacteria to coordinate their behavior (Davies et al., 1998; Sakuragi and Kolter, 2007; Bassler and Losick, 2006; Miller and Bassler, 2001; Herzberg et al., 2006; Laganenka et al., 2018; McLean et al., 1997; Paul et al., 2009). For example, in many cases, quorum sensing positively controls biofilm formation (Herzberg et al., 2006; Laganenka et al., 2018; Davies et al., 1998; Sakuragi and Kolter, 2007; McLean et al., 1997; González Barrios et al., 2006; Yarwood et al., 2004; Koutsoudis et al., 2006; Waters and Bassler, 2005; Parsek and Greenberg, 2005; Jayaraman and Wood, 2008; Hentzer et al., 2003; Kirisits and Parsek, 2006): autoinducer accumulation above a threshold concentration upregulates the expression of genes involved in biofilm formation, ultimately driving a transition from the planktonic to the biofilm state (Nadell et al., 2008). Again, however, the cellular factors that control this transition, such as the autoinducer production rate, diffusivity, and threshold concentration, can strongly vary between different species and environmental conditions.

Because planktonic dispersal and biofilm formation both depend on a dizzyingly complex array of factors, these distinct processes are typically studied in isolation. Thus, while each is well understood on its own, quantitative prediction of the conditions under which a population of planktonic bacteria transitions to the biofilm state—or instead, continues to disperse away and remains in the planktonic state—remains challenging. Here, we address this challenge by developing a mathematical model that describes essential features of motility-mediated dispersal of planktonic cells and autoinducer-mediated biofilm formation together. Using numerical simulations of this model, we systematically examine the influence of cellular concentration, motility, and chemotactic sensing; nutrient availability, diffusion, and consumption; and autoinducer production, diffusion, and accumulation on biofilm formation. Guided by these results, we establish a potentially-universal biophysical threshold that unifies the influence of all these factors in predicting the onset and extent of biofilm formation across different species and environmental conditions. Our work therefore provides a theoretical foundation for the prediction and control of biofilm formation in diverse and complex settings, and yields new quantitative predictions to guide future experiments.

Results

Development of the governing equations

As an illustrative example, and to connect our model to recent experiments of bacterial dispersal (Bhattacharjee et al., 2021), we consider a rectilinear geometry with a starting inoculum of planktonic cells at a maximal concentration b1,0 and of width x0. In general, the continuum variable b(x,t) describes the number concentration of bacteria, where x is the position coordinate and t is time, and the subscripts {1,2} represent planktonic or biofilm-associated cells, respectively. Following previous work (Lauffenburger, 1991; Keller and Segel, 1971; Adler, 1966a; Croze et al., 2011; Fu et al., 2018; Bhattacharjee et al., 2021), we consider a sole diffusible nutrient that also acts as the chemoattractant, with a number concentration represented by the continuum variable c(x,t) with diffusivity Dc. Initially, nutrient is replete throughout the system at a constant concentration c0. The bacteria then consume the nutrient at a rate b1κ1g(c), where κ1 is the maximum consumption rate per cell and the Michaelis-Menten function g(c)cc+cchar quantifies the nutrient dependence of consumption relative to the characteristic concentration cchar (Croze et al., 2011; Monod, 1949; Cremer et al., 2019; Woodward et al., 1995; Shehata and Marr, 1971; Schellenberg and Furlong, 1977; Cremer et al., 2016).

As time progresses, the bacteria thereby establish a local nutrient gradient that they respond to via chemotaxis (Figure 1A). In particular, planktonic cells disperse through two processes: undirected active diffusion with a constant diffusivity D1 (Berg, 2018), and directed chemotaxis with a drift velocity vcχ1log(1+c/c-1+c/c+) that quantifies the ability of the bacteria to sense and respond to the local nutrient gradient (Keller and Segel, 1970; Keller and Segel, 1971; Odell and Keller, 1976; Keller and Odell, 1975) with characteristic bounds c and c+ (Cremer et al., 2019; Sourjik and Wingreen, 2012; Shimizu et al., 2010; Tu et al., 2008; Kalinin et al., 2009; Shoval et al., 2010; Lazova et al., 2011; Celani et al., 2011; Fu et al., 2018; Dufour et al., 2014; Yang et al., 2015; Cai et al., 2016; Chen and Jin, 2011) and a chemotactic coefficient χ1. The planktonic cells also proliferate at a rate bγ1g(c), where γ1 is the maximal proliferation rate per cell. Finally, as the planktonic bacteria consume nutrients, they produce and secrete a diffusible autoinducer, with a number concentration represented by a(x,t) and with diffusivity Da, at a maximal rate k1 per cell. Motivated by some previous work (Hense et al., 2012; Hense and Schuster, 2015; Kirisits et al., 2007; Bollinger et al., 2001; Duan and Surette, 2007; Mellbye and Schuster, 2014; De Kievit et al., 2001; Pérez-Osorio et al., 2010), we take this process (hereafter referred to as ‘production’ for brevity) to also be nutrient-dependent via the same Michaelis-Menten function g(c) for the results presented in the main text, but we also consider the alternate case of ‘protected’ nutrient-independent production in the supplementary materials. Following previous work (Koerber et al., 2002; Ward et al., 2001; Ward et al., 2003), we also model natural degradation of autoinducer as a first-order process with a rate constant λ.

Figure 1 with 2 supplements see all
Competition between motility-mediated dispersal and autoinducer-mediated biofilm formation.

(A) Schematic of chemotactic dispersal: planktonic bacteria (green) consume nutrient (purple) and establish a local gradient that they, in turn, direct their motion in response to. (B) Schematic of positive quorum sensing-controlled biofilm formation: accumulation of produced autoinducer (red) above a threshold concentration causes cells to transition to the biofilm state (blue). (C) Results of an example simulation of Equations 1–4 showing the dynamics of the nutrient, planktonic cells, autoinducer, and biofilm cells from top to bottom, quantified by the normalized concentrations c/c0, b1/b1,0, a/a*, b2/b1,0, respectively; c0, b1,0, and a* represent the initial nutrient concentration, initial bacterial concentration, and autoinducer threshold for biofilm formation, respectively. The position coordinate is represented by the normalized position x/x0, where x0 is the width of the initial cellular inoculum. Different shades indicate different time points as listed. The inoculum initially centered about the origin consumes nutrient (purple), establishing a gradient that drives outward dispersal by chemotaxis (outward moving green curves); the cells also produce autoinducer (red) concomitantly. At t13 h, sufficient autoinducer has been produced to trigger biofilm formation at the origin; at even longer times (t16 h), nutrient depletion limits autoinducer production at this position. However, accumulation of autoinducer by the dispersing planktonic cells triggers partial biofilm formation at x/x06 as well. This competition between dispersal and biofilm formation leads to a final biofilm fraction of f=21% at the final time of t=20 h. An animated form of this figure is shown in Video 1. The values of the simulation parameters are given in Supplementary file 2.

As autoinducer is produced, it binds to receptors on the surfaces of the planktonic cells with a second-order rate constant α, as established previously (Koerber et al., 2002; Ward et al., 2001; Ward et al., 2003). Motivated by experiments on diverse bacteria, including the prominent and well-studied species Escherichia coli, Pseudomonas putida, and Pseudomonas aeruginosa (Davies et al., 1998; Sakuragi and Kolter, 2007; Bassler and Losick, 2006; Miller and Bassler, 2001; Herzberg et al., 2006; Laganenka et al., 2018; McLean et al., 1997; Paul et al., 2009; González Barrios et al., 2006; Yarwood et al., 2004; Koutsoudis et al., 2006; Waters and Bassler, 2005; Parsek and Greenberg, 2005; Jayaraman and Wood, 2008; Hentzer et al., 2003; Kirisits and Parsek, 2006), we assume that planktonic cells transition to the biofilm state at a rate τ1 when the local autoinducer concentration exceeds a threshold value a* (Figure 1B). Because our focus is on this transition, we assume that it is irreversible, and that cells in the biofilm lose motility. However, they still continue to consume nutrient, proliferate, and produce autoinducer with maximal rates κ2, γ2, and k2 per cell, respectively; additional behaviors such as subsequent production of extracellular polymeric substances or transitioning back to the planktonic state can be incorporated as future extensions to this model.

Hence, while planktonic cells can disperse via active diffusion and chemotaxis, their dispersal is hindered—and biofilm formation is instead promoted—when autoinducer accumulates sufficiently, as schematized in Figure 1A–B. The central goal of this paper is to examine the processes underlying this competition between dispersal and biofilm formation. Our model is thus summarized as:

(1) Planktonic: b1t=D12b1(b1vc)Motility+b1γ1g(c)Proliferation b1τ1H(aa)Biofilm formation
(2) Biofilm: b2t=b2γ2g(c)Proliferation+b1τ1H(aa)Biofilm formation
(3) Nutrient: ct=Dc2cDiffusion(b1κ1+b2κ2)g(c)Consumption
(4) Autoinducer: at=Da2aDiffusion+(b1k1+b2k2)g(c)Productiona(λ+αb1)Loss

where is the Heaviside step function describing the transition from the planktonic to biofilm state. To explore the competition between motility-mediated dispersal and autoinducer-mediated biofilm formation, we then numerically solve this system of coupled equations using values of all parameters—which are either intrinsic descriptors of cellular physiology or are solely/additionally influenced by the local environment—that are derived from experiments (Supplementary file 1). Further details are provided in the Materials and methods. Additional simulations indicate that the results obtained are not appreciably influenced by variations in the exact nature of how our model treats the arrest in planktonic cell motility while transitioning to the biofilm state (Figure 1—figure supplement 1) or the initial inoculum shape (Figure 1—figure supplement 2).

Representative numerical simulations

The results of a prototypical example are shown in Figure 1C and Video 1. Consumption by the planktonic cells (green curves) rapidly establishes a steep nutrient gradient (purple) at the leading edge of the inoculum. This gradient forces the planktonic cells to then move outward via chemotaxis. In particular, they self-organize into a coherent front that expands from the initial inoculum and continually propagates, sustained by continued consumption of the surrounding nutrient—consistent with the findings of previous studies of planktonic bacteria (Bhattacharjee et al., 2021). In this case, however, the cells also concomitantly produce autoinducer that accumulates into a growing plume (red). In some locations, the autoinducer eventually exceeds the threshold a*, thus driving the formation of an immobilized biofilm (blue). Hence, at long times, f=21% of the overall population is biofilm-associated, while the remaining 1-f=79% continues to disperse in the planktonic state.

Video 1
Animated form of Figure 1C: Results of an example simulation of Equations 1–4 showing the dynamics of the nutrient, planktonic cells, autoinducer, and biofilm cells from top to bottom, quantified by the normalized concentrations c/c0, b1/b1,0, a/a*, and b2/b1,0, respectively; c0, b1,0, and a* represent the initial nutrient concentration, initial bacterial concentration, and autoinducer threshold for biofilm formation, respectively.

The position coordinate is represented by the normalized position x/x0, where x0 is the width of the initial cellular inoculum. The inoculum initially centered about the origin consumes nutrient (purple), establishing a gradient that drives outward dispersal by chemotaxis (outward moving green curves); the cells also produce autoinducer (red) concomitantly. At t13 h, sufficient autoinducer has been produced to trigger biofilm formation at the origin; at even longer times (t16 h), nutrient depletion limits autoinducer production at this position. However, accumulation of autoinducer by the dispersing planktonic cells triggers partial biofilm formation at x/x04 as well. This competition between dispersal and biofilm formation leads to a final biofilm fraction of f=21% at the final time of t=20 h. The values of the simulation parameters are given in Supplementary file 2. The video displays the profiles every 30 min, to retain a manageable file size; however, the temporal step size in the actual simulations is 0.1 s.

Because the processes underlying motility-mediated dispersal and autoinducer-mediated biofilm formation are highly species- and environment-dependent, the values of the parameters in Equations 1–4 can span broad ranges—giving rise to different emergent behaviors under different conditions. Our simulations provide a way to examine how these behaviors depend on cellular concentration and motility, quantified by {b1,0,D1,χ1,c-,c+}, nutrient availability and consumption, quantified by {Dc,c0,κ1,κ2,cchar}, cellular proliferation, quantified by {γ1,γ2}, and autoinducer production, availability, and sensing, quantified by {Da,k1,k2,λ,α,τ,a*}. For example, implementing the same simulation as in Figure 1C, but for cells with faster nutrient consumption, yields a population that completely disperses in the planktonic state (the fraction of the population in the biofilm state at the final time of t=20 h is f=0%, as shown in Figure 2 and Video 2). Conversely, when cells consume nutrient slower, a larger fraction of the population forms an immobilized biofilm (f=52%, Figure 2—figure supplement 1 and Video 3).

Figure 2 with 1 supplement see all
Faster nutrient consumption limits autoinducer production, leading to complete dispersal.

Retion as in Figure 1C, but for planktonic cells with faster nutrient consumption (larger κ1). Panels and colors show the same quantities as in Figure 1C. The inoculum initially centered about the origin consumes nutrient (purple), establishing a gradient that drives outward dispersal by chemotaxis (outward moving green curves); the cells also produce autoinducer (red) concomitantly. However, nutrient is depleted at this position more rapidly, limiting autoinducer production; as a result, the population continues to disperse in the planktonic state and the final biofilm fraction is f=0%. An animated form of this figure is shown in Video 2. The values of the simulation parameters are given in Supplementary file 2.

Video 2
Animated form of Figure 2: Results of the same simulation as in Video 1, but for planktonic cells with faster nutrient consumption (larger κ1).

Panels and colors show the same quantities as in Video 1. The inoculum initially centered about the origin consumes nutrient (purple), establishing a gradient that drives outward dispersal by chemotaxis (outward moving green curves); the cells also produce autoinducer (red) concomitantly. However, nutrient is depleted at this position more rapidly, limiting autoinducer production; as a result, the population continues to disperse in the planktonic state and the final biofilm fraction is f=0%. The values of the simulation parameters are given in Supplementary file 2. The video displays the profiles every 30 min, to retain a manageable file size; however, the temporal step size in the actual simulations is 0.1 s.

Video 3
Animated form of Figure 2—figure supplement 1: Results of the same simulation as in Video 1, but for planktonic cells with slower nutrient consumption (smaller κ1).

Panels and colors show the same quantities as in Video 1. The inoculum initially centered about the origin slowly consumes nutrient (purple), establishing a slight gradient that allows partial planktonic dispersal (green curves moving outward); the cells also produce autoinducer (red) concomitantly. Because nutrient is consumed slowly, autoinducer production is not limited, resulting in partial biofilm formation (blue). Autoinducer has sufficiently accumulated above the threshold after t14 h, which causes a population of biofilm cells to form at the origin (x/x00). After 20 h, the biofilm population continues to grow, and additionally, autoinducer concentration exceeds the threshold concentration at x/x010. Thus, we see a second population of biofilm cells form, centered at x/x010. The slower nutrient consumption results in a greater final biofilm fraction than in Video 1—here, f=52%. The values of the simulation parameters are given in Supplementary file 2. The video displays the profiles every 30 min, to retain a manageable file size; however, the temporal step size in the actual simulations is 0.1 s.

Given that the competition between motility-mediated dispersal and autoinducer-mediated biofilm formation depends sensitively on such a bewildering array of cellular and environmental factors, we ask whether these dependencies can be captured by simple, generalizable, biophysical rules. Nondimensionalization of Equations 1–4 yields characteristic quantities and dimensionless groups that can parameterize these dependencies, as detailed in Appendix 1; however, given the large number of such groups, we seek an even simpler representation of the underlying processes that could unify the influence of all these different factors. To do so, we examine the fundamental processes underlying biofilm formation in our model.

Availability of nutrient for autoinducer production

When autoinducer production is nutrient-dependent, we expect that a necessary condition for biofilm formation is that enough nutrient is available for sufficient autoinducer to be produced to eventually exceed the threshold a*. To quantify this condition, we estimate two time scales: τd, the time taken by the population of planktonic cells to deplete all the available nutrient locally, and τa, the time at which produced autoinducer reaches the threshold for biofilm formation. While τd and τa can be directly obtained in each simulation, we seek a more generally-applicable analytical expression for both, solely using parameters that act as inputs to the model. In particular, for simplicity, we consider nutrient consumption and autoinducer production, both occurring at their maximal rates κ1 and k1, respectively, by an exponentially-growing population of planktonic cells that are uniformly distributed in a well-mixed and fixed domain. Integrating Equations 3 and 4 then yields (Appendix 2)

(5) τd=γ1-1ln(1+β~1,0)
(6) τa=γ1-1ln[1-ζ~1,0-1ln(1-η~)].

Three key dimensionless quantities, denoted by the tilde (~) notation, emerge from this calculation. The first, β~1,0γ1/(b1,0κ1/c0), describes the yield of new cells produced as the population consumes nutrient—quantified by the rates of cellular proliferation and nutrient consumption, γ1 and b1,0κ1/c0, respectively (Amchin et al., 2022). The second, η~αa/k1, describes the competition between autoinducer loss and production, quantified by their respective rates αa* and k1, at the single-cell scale. The third, ζ~1,0αb0/γ1, describes the loss of autoinducer due to cell-surface binding as the population continues to grow, quantified by the population-scale rates of autoinducer loss and cellular proliferation, αb0 and γ1, respectively; for simplicity, this quantity neglects natural degradation of autoinducer, given that the degradation rate is relatively small, with λαb0.

The ratio between Equations 5 and 6 then defines a nutrient availability parameter, D~τd/τa. When D~ is large, produced autoinducer rapidly reaches the threshold for biofilm formation before the available nutrient is depleted; by contrast, when D~ is small, nutrient depletion limits autoinducer production. Hence, we hypothesize that D~D~ specifies a necessary condition for biofilm formation, where D~ is a threshold value of order unity. The simulations shown in Figures 1C and 2 and Figure 2—figure supplement 1 enable us to directly test this hypothesis. Consistent with our expectation, the simulation in Figure 1C is characterized by D~=0.33, near the expected threshold for biofilm formation; as a result, f=21%. When consumption is faster as in Figure 2 (D~=0.033), the available nutrient is rapidly depleted; thus, cells disperse away before sufficient autoinducer is produced to initiate biofilm formation, and f=0%. Conversely, when nutrient consumption is slow as in Figure 2—figure supplement 1 (D~=3.1), nutrient continues to be available for autoinducer production, eventually driving biofilm formation, with a larger fraction f=52%.

Taken together, these results support our hypothesis that D~D~1 is a necessary condition for biofilm formation. It is not, however, a sufficient condition: repeating the simulation of Figure 1C but for faster-moving cells yields a population that rapidly disperses without forming a biofilm at all (f=0%, Figure 3A and Video 4)—despite having the same value of D~=0.33. Thus, our mathematical description of the conditions that determine biofilm formation is, as yet, incomplete.

Figure 3 with 1 supplement see all
Enhanced motility enables cells to disperse before sufficient autoinducer accumulates, leading to complete dispersal.

(A) Results of the same simulation as in Figure 1C, but for faster-moving planktonic cells (larger D1 and χ1). Panels and colors show the same quantities as in Figure 1C. The inoculum initially centered about the origin consumes nutrient (purple), establishing a gradient that drives outward dispersal by chemotaxis (outward moving green curves); the cells also produce autoinducer (red) concomitantly. More rapid dispersal enables the planktonic cells to ‘outrun’ the growing autoinducer plume, as shown by the extended and magnified view in (B). As a result, the population continues to disperse in the planktonic state and the final biofilm fraction is f=0%. An animated form of this figure is shown in Video 4. The values of the simulation parameters are given in Supplementary file 2.

Video 4
Animated form of Figure 3: Results of the same simulation as in Video 1, but for faster-moving planktonic cells (larger D1 and χ1).

Panels and colors show the same quantities as in Video 1. The inoculum initially centered about the origin consumes nutrient (purple), establishing a gradient that drives outward dispersal by chemotaxis (outward moving green curves); the cells also produce autoinducer (red) concomitantly. More rapid dispersal enables the planktonic cells to ‘outrun’ the growing autoinducer plume, as shown by the extended and magnified view in (B). As a result, the population continues to disperse in the planktonic state and the final biofilm fraction is f=0%. The values of the simulation parameters are given in Supplementary file 2. The video displays the profiles every 30 min, to retain a manageable file size; however, the temporal step size in the actual simulations is 0.1 s.

Competition between motility-mediated dispersal and autoinducer accumulation

The results shown in Figure 3 indicate that the ability of planktonic bacteria to move, which is not incorporated in the nutrient consumption parameter D~, also plays a key role in regulating whether a biofilm forms. Indeed, close inspection of Figure 3A hints at another necessary condition for biofilm formation: as shown by the magnified view in Figure 3B (e.g., at t=4 h), the leading edge of the dispersing planktonic cells extends beyond the plume of produced autoinducer. Therefore, we expect that even when sufficient nutrient is available for autoinducer production (D~D~1), autoinducer production must be rapid enough to reach the threshold for biofilm formation before cells have dispersed away. To quantify this condition, we estimate the the time τc at which the motile planktonic cells begin to ‘outrun’ the growing autoinducer plume. Specifically, we quantify the dynamics of the leading edge positions of the chemotactic front of planktonic cells and the autoinducer plume, x1,edge(t) and xa,edge(t), respectively. The front position x1,edge(t) is known to depend on cellular motility, nutrient diffusion, and nutrient consumption in a non-trivial manner (Berg, 2004; Cremer et al., 2019; Fu et al., 2018; Amchin et al., 2022), and we are not aware of a way to compute this quantity a priori from input parameters; instead, we extract this sole quantity from each simulation by identifying the largest value of x at which b110-4b1,0. While the plume position xa,edge(t) can also be directly obtained in each simulation, we again develop a more generally applicable analytical expression by assuming that the autoinducer continually diffuses from the initial inoculum: xa,edge(t)=x0+2Dat. Then, τc can be directly determined as the time at which x1,edge(t) begins to exceed xa,edge(t).

The ratio between τc thereby determined and τa, the time required for produced autoinducer to reach the threshold for biofilm formation (Equation 6), then defines a cellular dispersal parameter, J~τc/τa. When J~ is large, autoinducer accumulation is sufficiently rapid to drive biofilm formation; by contrast, when J~ is small, the planktonic cells rapidly disperse without forming a biofilm. Hence, we hypothesize that J~J~ specifies another necessary condition for biofilm formation, where J~ is, again, a threshold value of order unity. The simulations shown in Figures 1C and 3A enable us to directly test this hypothesis. Consistent with our expectation, the simulations in Figure 1C and Figure 2—figure supplement 1 are characterized by J~=1.6, near the expected threshold for biofilm formation; as a result, f>0 in both cases. Furthermore, implementing the same simulation as Figure 1C (with the same D~=0.33) but for slower-moving cells, characterized by a larger J~=120, yields a population that forms an even larger biofilm fraction f=82% (Figure 3—figure supplement 1 and Video 5). Conversely, when cellular dispersal is faster as in Figure 3, characterized by a smaller J~=0.1, the cells disperse away before sufficient autoinducer is produced to initiate biofilm formation, and f=0%. Taken together, these results support our hypothesis that J~J~1 is another necessary condition for biofilm formation.

Video 5
Animated form of Figure 3—figure supplement 1: Results of the same simulation as in Video 1, but for slower-moving planktonic cells (smaller D1 and χ1).

Panels and colors show the same quantities as in Video 1. The inoculum initially centered about the origin consumes nutrient (purple), establishing a slight gradient—however, because the motility parameters are diminished, the planktonic population (green) remains around the origin. The planktonic cells produce autoinducer (red) concomitantly, and after 1 h, the autoinducer concentration exceeds the threshold concentration. Thus, some of the planktonic cells transition to biofilm cells, centered at the origin. Both the biofilm cells and planktonic cells continue to grow, produce autoinducer, and consume nutrient; the planktonic cells do not disperse due to their diminished motility, resulting in a larger fraction of biofilm cells (f=82%) than in Video 1. The values of the simulation parameters are given in . The video displays the profiles every 30 min, to retain a manageable file size; however, the temporal step size in the actual simulations is 0.1 s.

A universal biophysical threshold for biofilm formation

Thus far, we have shown that the two conditions D~D~ and J~J~ are both necessary for biofilm formation. Is the combination of both sufficient to fully specify the conditions required for biofilm formation? To test this possibility, we implement 10,983 numerical simulations of Equations 1–4 exploring the full physiological ranges of the input parameters that describe cellular, nutrient, and autoinducer properties for diverse bacterial species/strains and environmental conditions (Supplementary file 1). For each simulation, we compute D~, J~, and f. Remarkably, despite the extensive variability in the values of the underlying parameters, all the results cluster between two states parameterized by D~ and J~, as shown in Figure 4A: motility-mediated dispersal without biofilm formation (f=0%, green points) when either D~<D~ or J~<J~, and biofilm formation without dispersal (f=100%, blue points) when both D~>D~ and J~>J~. Many different combinations of the input parameters yield the same (D~, J~); yet, no matter the input values of these parameters, which vary over broad ranges for different cells and environmental conditions, (D~, J~) uniquely specify the resulting biofilm fraction f for all points, as shown in Figure 4B–C—indicating that these two dimensionless parameters reasonably encompass all the factors determining biofilm formation within our model. We observe some exceptions at the boundary between these two states, likely because the simplifying assumptions underlying the derivation of the D~ and J~ parameters begin to break down. Nevertheless, the boundary between both states, summarized by the relation D~/D~+J~/J~1 with D~ and J~ both 1 (black curve), thus specifies a universal biophysical threshold for biofilm formation.

Figure 4 with 2 supplements see all
The two states of complete dispersal by planktonic cells (green) and complete formation of a biofilm (blue) can be universally described by three dimensionless parameters.

(A) State diagram showing the fraction of biofilm formed, f, at the final time (t=20 h) for different values of the nutrient availability and cellular dispersal parameters, D~ and J~, respectively. The state diagram summarizes the results of 10,983 simulations of Equations 1–4 exploring the full range of parameter values describing different bacterial species/strains and different environmental conditions (Supplementary file 1). Each point represents the mean value of f obtained from multiple simulations with different parameter values, but with similar D~ and J~ (identical within each bin defined by the spacing between points). (B) represents the same data, but each point represents the standard deviation of the values of f obtained from the same simulations. Despite the vastly differing conditions explored in each simulation, they cluster into the two states of planktonic dispersal (green) and biofilm formation (blue) when parameterized by D~ and J~. The boundary between the two states can be described by the relation D~/D~+J~/J~1, as shown by the black line; this relation combines the transition between the two states that occurs at both D~1 and J~1. Away from this boundary, all simulations for the same D~ and J~ collapse to have the same biofilm fraction f, as shown by the points in (B) and examples (i)–(iii) and (v) in (C)—confirming the universality of our parameterization. Near the boundary, we observe some slight differences between simulations, as shown in (B) and examples (iv) and (vi) in (C). The values of the simulation parameters for the examples in (C) are given in Source Data file 1. The data in (AC) correspond to a fixed value of the third dimensionless parameter S~=50, which describes the case of biofilm cells that produce autoinducer rapidly; repeating these simulations for the opposite case of slow autoinducer production by biofilm cells (S~=1/50) yields the state diagram shown in (D), but for 14,351 simulation runs; again, (E) shows the standard deviation of the corresponding values of f. As shown by (DE), while the transition between the two states (black line) is unaffected by the change in S~, the transition to complete biofilm formation is more gradual. Together, the three parameters D~, J~, and S~ provide a full description of the onset and extent of biofilm formation across vastly different conditions.

Discussion

The transition from the planktonic to biofilm state is known to depend on a large array of factors that describe cellular concentration, motility, and proliferation; nutrient availability and consumption; and autoinducer production, availability, and sensing—all of which can vary considerably for different strains/species of bacteria and environmental conditions. Therefore, quantitative prediction of the onset of biofilm formation is challenging. The biophysical model presented here provides a key step toward addressing this challenge. In particular, for the illustrative case we consider—in which cells can either disperse through active motility, retaining them in the planktonic state, or form an immobilized biofilm when exposed to sufficient autoinducer—we have shown that the onset of biofilm formation is uniquely specified by a biophysical threshold set by the two dimensionless parameters D~ (quantifying nutrient availability) and J~ (quantifying bacterial dispersal). Importantly, within the formulation of our model, this threshold is universal: many different combinations of cellular and environmental factors are described by the same (D~,J~), and thus, yield the same onset of biofilm formation. Therefore, given a bacterial strain and set of environmental conditions, extensions of our model could help provide a way to predict whether a biofilm will form a priori. Indeed, because the factors that define D~ and J~ can be directly measured, our work now provides quantitative principles and predictions (as summarized in Figure 4) to guide future experiments.

For generality, our model also incorporates proliferation, nutrient consumption, and autoinducer production by cells after they have transitioned to the biofilm state. Hence, within our model, biofilm-produced autoinducer could also drive surrounding planktonic cells to transition to the biofilm state. In this case, we expect that the long-time fraction of the population in the biofilm state, f, will also depend on nutrient depletion and autoinducer production by the growing biofilm. Indeed, performing a similar calculation as that underlying the nutrient availability parameter, D~, yields a third dimensionless parameter, S~τd,2/τa,2; here, τd,2 and τa,2 describe the times at which biofilm cells have depleted all the available nutrient and produced enough autoinducer to reach the threshold for biofilm formation, respectively (Appendix 2). Thus, we hypothesize that, while the onset of biofilm formation is specified by (D~,J~), the final extent of biofilm that has formed will also be described by S~. The results shown in Figures 14 have a fixed S~=50, which describes the case of a biofilm that produces autoinducer rapidly; repeating these simulations for the opposite case of slow autoinducer production by biofilm cells, with S~=1/50, yields the state diagram shown in Figure 4D. In agreement with our hypothesis, while the transition to the biofilm state (black line) is not appreciably altered by the change in S~, the transition to complete biofilm formation (f=1) is more gradual in this case (compare Figure 4A,B,D,E). Moreover, we note that our analysis thus far has focused on the case in which autoinducer production is nutrient-dependent; however, this process may sometimes be nutrient-independent (Narla et al., 2021). In this case, we expect that our overall analysis still applies, but with the onset of biofilm formation specified by only the dispersal parameter J~—as confirmed in Figure 4—figure supplement 1.

Possible extensions of our work

The transition from the planktonic to biofilm state is highly complex and, in many cases, has features that are unique to different species of bacteria. Nevertheless, our model provides a minimal description that can capture many of the essential features of biofilm formation more generally—thereby providing a foundation for future extensions of our work, some of which are described below.

  1. For simplicity, our model considers only one spatial dimension; however, fascinating new effects may arise in higher-dimensional implementations of our model. For example, in our prior work modeling the collective migration of planktonic bacteria in the absence of quorum sensing-mediated biofilm formation, we found that variations in the shape of the cellular front orthogonal to the main propagation direction ‘smooth out’ over time (Alert et al., 2022; Bhattacharjee et al., 2022). In particular, cells at outward-bulging parts of the front are exposed to more nutrient, which diminishes their ability to respond to the nutrient gradient via chemotaxis and thus slows them down. As a result, the migrating front eventually smooths to a flat shape whose subsequent dynamics can then be described using just one spatial dimension, just as in our treatment here. However, we expect that this behavior could be altered in interesting new ways when the cells can additionally produce and sense autoinducer and thereby transition to the biofilm state, as is the case here. In this case, we speculate that because cells at outward-bulging parts of the front are exposed to more nutrients and have a weaker chemotactic response, autoinducer production and accumulation will be more rapid relative to cellular dispersal. That is, at these parts of the front, τa and τc will be shorter and longer, respectively, causing the dispersal parameter J~ to be larger locally. Thus, our model would predict biofilm formation to occur first at these parts of the front, potentially also influencing subsequent dispersal and biofilm formation at other locations along the front. Therefore, while our conclusions here could be the same locally at different parts of the front, the global behavior of the population could be different—potentially giving rise to e.g., spatially-heterogeneous biofilm formation.

  2. As an illustrative example, our model considers the case in which cells produce a single autoinducer; however, some quorum sensing systems utilize multiple autoinducers (Miller and Bassler, 2001; Miller et al., 2002; Pesci et al., 1997), which could be described using additional field variables and equations similar to Equation 4. Moreover, while we take the nutrient to be the sole chemoattractant, in some cases, autoinducers can also act as chemoattractants (Laganenka et al., 2016), which could also be described in our framework by e.g., introducing autoinducer-dependent chemotaxis in the drift velocity in Equation 1.

  3. Our model considers positive quorum sensing control in which planktonic cells transition to the biofilm state in a step-like fashion when the local autoinducer concentration exceeds a threshold value. That is, when planktonic cells encounter sufficiently concentrated autoinducer, the diffusivity and chemotactic coefficient transition in a step-like fashion from the constant values D1 and χ1, respectively, to zero after the time duration τ, for simplicity. In real systems, the change in cellular motility may not be as temporally abrupt. Future work could address a more gradual loss of motility in our theoretical framework by, for example, considering a cellular diffusivity and chemotactic coefficient that gradually transition from their planktonic values to zero over a non-zero time scale. Given that the same cells would be transitioning from the motile planktonic to immotile biofilm state—but in this case with the introduction of a time-varying diffusivity and chemotactic coefficient—we expect that the long-time biofilm fraction f will be similar, and only the spatial profile of the biofilm population may be altered. Hence, we expect that our main findings summarized in Figure 4 will be unaffected by such a change. Indeed, performing the same representative simulation shown in Figure 1C, but with both motility parameters D1 and χ1 smoothly transitioning to zero in time, shows nearly identical results (Figure 1—figure supplement 1)—confirming our expectation that the temporal nature of the arrest in motility does not appreciably influence our model results and conclusions.

  4. While we take the transition to the biofilm state as being irreversible, this is often not the case (Barraud et al., 2006; Kaplan, 2010; Abdel-Aziz, 2014). Longer-time transitions back to the planktonic state could be described using additional terms similar to the last terms of Equations 1; 2, but with the opposite sign. Similar modifications could be made to describe other species of bacteria (e.g., Vibrio cholerae) that utilize the opposite case of negative quorum sensing control, in which biofilm cells instead transition to the planktonic state when the autoinducer accumulates above a threshold value (Hammer and Bassler, 2003; Bridges and Bassler, 2019).

  5. Biofilms are often formed by multiple different microbial species, whereas our model describes biofilm formation by a single species, for simplicity. Nevertheless, we expect that our theoretical framework can be extended by following reasoning similar to that described in this paper, but with the introduction of additional equations and variables in the governing Equations 1–4 to describe the distinct cell and chemical types, as appropriate. For example, if the different species i consume and respond to distinct nutrients ci, and secrete and respond to distinct autoinducers ai, each species could be described in isolation using our same governing Equations 1–4, but now extended to incorporate the distinct variables ci, ai, b1,i, and b2,i. Then, directly following our approach, each species would be described by its own dimensionless parameters D~i and J~i, with D~/D~i+J~/J~i1 again specifying the threshold for biofilm formation for each. We hypothesize that the composition of the final two-species biofilm community would then be given by the combination of each single-species biofilm. Alternatively, in the case that the different species consume and respond to the same nutrient c, and secrete and respond to the same autoinducer a, our Equations 1–4 could again be extended to consider the cellular parameters specific to each species i. In this approach, however, biofilm formation by each of the species cannot be described in isolation, because they are coupled through the nutrient and autoinducer dynamics. Instead, the calculations of the characteristic time scales τd, τa, and τc would need to be extended, following our approach, to now reflect contributions from all the different species. We hypothesize that the overall multi-species community would then be described by one set of governing dimensionless parameters (D~,J~), and D~/D~+J~/J~1 would again specify a universal biophysical threshold for the onset of biofilm formation for the overall community—but the composition of the final multi-species biofilm that results above this threshold may not be uniquely specified by (D~,J~).

  6. Biofilm formation may be regulated by other, non-quorum sensing-based, processes not considered in our model. For example, the intracellular accumulation of secondary signaling molecules such as cyclic di-GMP can also regulate biofilm formation (Valentini and Filloux, 2016; Simm et al., 2004; Jenal et al., 2017; Römling et al., 2013; Hengge, 2009; Krasteva et al., 2010; Baraquet and Harwood, 2013; Trampari et al., 2015; Davis et al., 2013; Boehm et al., 2010; Russell et al., 2013). In some cases, this process may be controlled by quorum sensing (Waters et al., 2008) and thus could be described by our model, while in others, it is controlled by other cues such as e.g., contact with surfaces, which would need to additionally be incorporated into our theoretical framework.

  7. Finally, we note that our model is deterministic—describing the cellular processes of motility, nutrient consumption, proliferation, and autoinducer production, availability, and sensing using the parameters {D1,χ1,c-,c+,κ1,κ2,cchar,γ1,γ2,k1,k2,α,τ,a*}, respectively. Each of these is taken to be single-valued in each of our simulations. However, these parameters can have a distribution of values arising from e.g., inherent cell-to-cell variability. Because these values define the governing D~ and J~ that specify the threshold for biofilm formation in our model, we expect that variability in the parameter values would broaden the planktonic-to-biofilm transition predicted by our model. That is, we expect the transition specified by the black curve in Figure 4A to be smeared out, similar to what is seen in Figure 4D, though due to a fundamentally different reason—with biofilm formation arising in some cases at lower (D~,J~) than predicted by the black curve. Indeed, similar behavior was recently observed in a distinct model of biofilm formation on flat surfaces (Sinclair et al., 2022). Exploring the influence of such variations by using a more probabilistic approach in our theoretical framework will thus be a useful direction for future research.

Materials and methods

To numerically solve the continuum model described by Equations 1–4, we follow the experimentally validated approach used in our previous work (Bhattacharjee et al., 2021; Amchin et al., 2022). Specifically, we use an Adams-Bashforth-Moulton predictor-corrector method in which the order of the predictor and corrector are 3 and 2, respectively. Because the predictor-corrector method requires past time points to inform future steps, the starting time points must be found with another method; we choose the Shanks starter of order 6 as described previously (Rodabaugh and Wesson, 1965; Shanks, 1966). For the first and second derivatives in space, we use finite difference equations with central difference forms in rectilinear coordinates. The temporal and spatial resolution of the simulations are δt=0.1 s and δx=20 μm, respectively; furthermore, we constrain our analysis to simulations for which the peak of the overall bacteria population moves slower than δx/δt. Repeating representative simulations with different spatial and temporal resolution indicates that even finer discretization does not appreciably alter the results (Figure 4—figure supplement 2). Thus, our choice of discretization is sufficiently finely-resolved such that the results in the numerical simulations are not appreciably influenced by discretization. Furthermore, performing the same representative simulation shown in Figure 1C, but with the shape of the initial inoculum changed from a Gaussian profile to a step function with the same maximum cellular concentration and width, shows nearly identical results (Figure 1—figure supplement 2)—suggesting that our results are robust to variations in this initial condition chosen. Further probing the mathematical structure of our biophysical model to examine additional influences of initial conditions and explore the possibility of oscillatory solutions, closed orbits, or singularities would be a fascinating direction for future work.

To connect the simulations to our previous experiments (Bhattacharjee et al., 2021), we choose a total extent of 1.75×104 μm for the size of the entire simulated system, with no-flux conditions for the field variables b1, b2, c, and a applied to both boundaries at x=0 and 1.75×104 μm. As in the experiments, we initialize each simulation with a starting inoculum of planktonic cells with a Gaussian profile defined by the maximum concentration b1,0 at x=0µm and a full width at half maximum of 100 μm. Nutrient is initially uniform at a fixed concentration c0, and the autoinducer and biofilm concentrations are initially zero, throughout. Furthermore, following previous work (Amchin et al., 2022; Dell’Arciprete et al., 2018; Volfson et al., 2008; Farrell et al., 2013; Klapper and Dockery, 2002; Head, 2013), we also incorporate jammed growth expansion of the population in which growing cells push outward on their neighbors when the total concentration of bacteria is large enough. In particular, whenever the total concentration of bacteria (planktonic and biofilm) exceeds the jamming limit of 0.95 cells μm3 at a location xi, the excess cell concentration is removed from xi and added to the neighboring location, xi+δx, where δx represents the spatial resolution of the simulation, retaining the same ratio of planktonic to biofilm cells in the new location. We repeat this process for every location in the simulated space for each time step.

We run each simulation for a total simulated duration of tsim=20 h. At this final time, we use the simulation data to directly compute fb2dxb2dx+b1dx, the total fraction of the population in the biofilm state. We also compute the values of the dimensionless parameters D~, J~, and S~ using the equations presented in the main text. We note that the autoinducer production time τa (Equation 6) is only finite for η~αa/k1<1; when η~1, the rate of autoinducer loss exceeds that of autoinducer production, and thus the time required to reach the threshold for biofilm formation diverges. Because both D~ and J~ are defined as τd/τa and τc/τa, respectively, for simulations with η~1, we represent them on the state diagrams in Figure 4 and Figure 4—figure supplement 1 at (D,J)=(102,103), the smallest values shown on the diagrams. All of these simulations have f=0, as expected. Furthermore, to ensure that tsim is sufficiently long, we (i) only perform simulations with τa and τa,2 smaller than tsim, and (ii) do not include simulations with f=0 but τc=τsim, for which sufficient time has not elapsed for planktonic cells to chemotactically disperse.

Appendix 1

Nondimensionalizing the governing equations

The governing equations Equations 1–4 are described by six variables: those describing the concentrations of planktonic bacteria (b1), biofilm bacteria (b2), nutrient (c), autoinducer molecules (a), as well as the one-dimensional space (x), and time (t) coordinates. Additional constants for our equations are highlighted in Supplementary file 1, with initial conditions b1(t=0)=b1,0, c(t=0)=c0, and x0 as the width of the initial planktonic inoculum. We define the dimensionless variables b~1b1B1, b~2b2B2, c~cC, a~aA, x~xX, and t~tT, where the tilde (~) notation indicates a dimensionless quantity and the dimensional quantities B1, B2, C, A, X, and T are not specified a priori. Thus, in nondimensional form, Equations 1–4 can be represented as:

(S1) Planktonic: b~1t~=D1(X2/T)~2b~1χ1(X2/T)~(b~1~log(1+c~/c~1+c~/c~+))+ b~1(γ1T)g(c~)b~1(τ1T)H(a~a/A)
(S2) Biofilm: b~2t~=(γ2T)b~2g(c~)+(B1/B2)(τ1T)b~1H(a~a/A)
(S3) Nutrient: c~t~=Dc(X2/T)~2c~((B1/C)κ1Tb~1+(B2/C)κ2Tb~2)g(c~)
(S4) X2/D1,X2/chi1,X2/Dc,X2/Daγ11,γ21,τ,cob1, 0k1,cob1, 0k2,ab1, 0k1,ab1, 0k2,1,(αβ1,0)1

where g(c~)c~c~+c~char. Given that the characteristic autoinducer concentration a* arises in the argument of the Heaviside step function in Equations S1 and S2, we choose A=a*. Moreover, given that the planktonic cells have a characteristic concentration b1,0 defined by the initial inoculum, we choose B1=b1,0. The fraction of the population in the biofilm state is defined as f=b2/(b2+b1); thus, to ensure that f~=f for simplicity, we also choose B2=B1=b1,0. Finally, given that the nutrient has a characteristic concentration c0 defined by the initial saturation, we choose C=c0. With these choices of characteristic quantities, multiple length and time scales emerge as possible choices for X and T, respectively:

 Length scale:TD1, TDa, TDc, TDa

 Time scale:X2/D1,X2/chi1,X2/Dc,X2/Da,γ11,γ21,τ,cob1, 0k1,cob1, 0k2,ab1, 0k1,ab1, 0k2,1,(αβ1,0)1

 Each such choice will lead to the emergence of many different dimensionless groups characterizing this problem. Nevertheless, all these different groupings are accounted for in the dimensionless parameters D~, J~, and S~ described in the main text, with the exception of quantities involving the nutrient diffusivity Dc, planktonic-to-biofilm transition rate τ-1, and the natural autoinducer degradation rate λ, which have corresponding time scales that are much smaller than the other time scales of the systems considered here and are neglected from our analysis for simplicity.

Appendix 2

Derivation of the dimensionless parameters D~ and S~

We first estimate the time τd taken for cells to deplete available nutrient through consumption. To do so, for simplicity, we consider a population of planktonic cells exponentially growing at the maximal rate γ1, uniformly distributed in a well-mixed and fixed domain (i.e., neglecting motility-mediated spreading), and consuming nutrient at the maximal rate κ1. Thus, dcdt=-κ1b1,0etγ1; integrating this equation from t=0 (with c=c0) to t=τd (with c=0) yields Equation 5 of the main text.

We use a similar approach to estimate the time τa taken for produced autoinducer to reach the threshold for biofilm formation a*. In particular, we consider the same population of planktonic cells secreting autoinducer at the maximal rate k1. We neglect natural degradation of autoinducer, given that the degradation rate is relatively small compared to binding to the cell surface receptors with a second-order rate constant α, that is, λαb0. The rate of autoinducer production and loss are then given by b1,0etγ1×k1 and b1,0etγ1×αa, respectively, ultimately yielding dadt=b1,0etγ1(k1-αa). Integrating this equation from t=0 (with a=0) to t=τa (with a=a*) then yields Equation 6 of the main text. Notably, this analytical solution for the time scale τa is only defined for η~αa/k1<1; when η~1, the rate of autoinducer loss exceeds that of autoinducer production and secretion, and thus the time required to reach the threshold for biofilm formation diverges. Finally, the ratio of τd and τa thus derived yields the nutrient availability parameter D~ as described in the main text.

Thus far, we have only considered nutrient consumption by planktonic bacteria. However, cellular proliferation, autoinducer production, and nutrient consumption can also occur for cells after they have transitioned to the biofilm state, causing biofilm-produced autoinducer to also drive surrounding planktonic cells to transition to the biofilm state. Hence, we repeat the same calculations for τa and τd as described above, but now for a population of cells in the biofilm state (still with the characteristic concentration b1,0 defined in our model), exponentially growing at the maximal rate γ2 and consuming nutrient at the maximal rate κ2. In this case, dcdt=-κ2b1,0etγ2, and integrating this equation from t=0 (with c=c0) to t=τd,2 (with c=0) yields τd,2=γ2-1ln(1+β~2,0), where β~2,0γ2/(b1,0κ2/c0) describes the yield of new biofilm cells produced as the population consumes nutrient. For the calculation of autoinducer production, we adopt a similar approach as that described above to calculate τa, but now assuming that the biofilm surface receptors are saturated (i.e., neglecting autoinducer loss). As a result, dadt=b1,0k2etγ2. Integrating this equation from t=0 (with a=0) to t=τa,2 (with a=a*) finally yields τa,2=γ2-1ln(1+θ~2,0), where θ~2,0γ2b1,0k2/a*. The ratio of τd,2 and τa,2 thus derived then yields the parameter S~ as described in the main text.

Data availability

All data generated or analyzed during this study are included in the manuscript and supporting file; source data files have been provided for all figures.

References

  1. Report
    1. Rodabaugh D
    2. Wesson JR
    (1965)
    On the Efficient Use of Predictor-Corrector Methods in the Numerical Solution of Differential Equations. NASA-TN-D-2946
    WashingtonWASHINGTON, D.C: National Aeronautics and Space AdministrationNATIONAL AERONAUTICS AND SPACE ADMINISTRATION.

Decision letter

  1. Raymond E Goldstein
    Reviewing Editor; University of Cambridge, United Kingdom
  2. Aleksandra M Walczak
    Senior Editor; CNRS LPENS, France

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Decision letter after peer review:

Thank you for submitting your article "A biophysical threshold for biofilm formation" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Aleksandra Walczak as the Senior Editor. The reviewers have opted to remain anonymous. We regret the lengthy delay in furnishing this report.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) Given that the continuum equations are by definition evolution dynamics for the mean values of the various quantities, are there situations in which a more probabilistic approach (e.g. by a Fokker-Planck equation) ought to be considered?

2) The problem studied here is simplified to one spatial dimension. How might the conclusions change if there were variations orthogonal to the main front propagation direction? Could this lead to qualitatively different conclusions?

3) Please clarify the use of a single effective diffusivity parameter for a system that is slowly being arrested. That is, is the diffusion coefficient a constant in equation 1? The effective diffusivity of cells changes dramatically from the planktonic to the immobilized state. How is this issue resolved in the model?

4) Biofilms are usually formed by different types of micro-organisms, and I am struggling to see whether this single species model is general enough to capture that reality. We do not mean to ask for the authors to perform more simulations, but we would like to hear know their thoughts on this issue.

5) While we believe that the results are compelling and quite useful, the work can be strengthened. For example, biofilms are usually heterogeneous (i.e., different microbial populations) and many of the parameters in the model are homogeneous. It would be interesting to understand how even a second cell species would alter (or not) the results, in particular the claim that onset of biofilm formation can be uniquely defined by the parameters D and J.

6) Are these results sensitive to the chosen initial condition? That is, it would be good to explore the dynamics of the biophysical model (perhaps that has been done already, but we did not see it in the manuscript). For example, are there oscillatory solutions or closed orbits, or even singularities that arise from the solution of the equations?Reviewer #1 (Recommendations for the authors):

The problem of biofilm formation is considered here in the context of a well-defined geometry (rectilinear) in which there are spatio-temporally evolving populations of planktonic and surface-bound cells in the biofilm, consuming nutrients, dividing and responding to secreted autoinducers. While there are many individual parameters that describe the various processes involved, including diffusion constants, motilities, growth, consumption, and transitions rates, the main result of the paper is that only two dimensionless combinations of parameters determine the threshold for biofilm formation. These are a ratio of timescales for nutrient depletion by planktonic cells and that for autoinducer to reach threshold defines one parameter, while the second is a ratio of times associated with the moving front of transitioning cells and the autoinducer timescale. Investigating a very wide range of the individual parameters in the problem shows that these two composite ones serve as a reliable predictor of biofilm formation.

This is an excellent paper that provides real insight into a notoriously difficult problem in biology. While the individual assumptions of the model can certainly be tweaked or generalised, the overall framework, and in particular the simplified geometry, should have a great impact on the field.Reviewer #2 (Recommendations for the authors):

The authors investigate the onset of biofilm formation (from a planktonic cell state) using a set of coupled, nonlinear differential equations. The equations model the two main cell states, planktonic and biofilm; two auxiliary equations for nutrient consumption and autoinducer production are introduced to model the transition from planktonic to biofilm state. Analysis of the biophysical model shows that there are two main dimensionless parameters that govern the onset of biofilm formation (in the model), namely a nutrient (D) and an autoinducer (J) availability parameters. An important result is that the biophysical threshold for biofilm formation seems quite robust across many different sets of model parameters suggesting, perhaps, a universal threshold. A strength of the formulation is that these dimensionless parameters can be experimentally measured, and therefore the presented biophysical model can be tested in the laboratory. Overall, the results from the model support the authors' claims.

https://doi.org/10.7554/eLife.76380.sa1

Author response

Essential revisions:

1) Given that the continuum equations are by definition evolution dynamics for the mean values of the various quantities, are there situations in which a more probabilistic approach (e.g. by a Fokker-Planck equation) ought to be considered?

The Reviewers and Editors raise an excellent point, and we thank them for this constructive suggestion. Our (deterministic) framework provides a first step toward modeling the process of biofilm formation. However, we fully agree that a probabilistic approach that can additionally describe e.g., stochastic variations in cellular behaviors would be a valuable extension of our work. Indeed, motivated by this insightful suggestion, we have started to explore how to incorporate such an approach into our model.

For example: our model describes the cellular processes of motility, nutrient consumption, proliferation, and autoinducer production, availability, and sensing using the parameters {D1,x1,c_,c+,k,cchar,γ,k,α,τ,a*}, respectively—each of which is taken to be single-valued in each of our simulations. However, these parameters can have a distribution of values arising from e.g., inherent cell-to-cell variability in the underlying processes they describe. Because these values define the governing dimensionless parameters D~ and J~ that specify the threshold for biofilm formation in our model, we expect that variability in the parameter values would broaden the planktonic-to-biofilm transition predicted by our model. That is, we expect the transition specified by the black curve in Figure 4A to be smeared out, similar to what is seen in Figure 4D, though due to a fundamentally different reason—with biofilm formation arising in some cases at lower (D~, J~) than predicted by the black curve. Indeed, similar behavior was recently observed in a distinct model of biofilm formation on flat surfaces [Ref. 118 newly added to our manuscript]. Exploring the influence of such variations within our theoretical framework will thus be a useful direction for future research.

We appreciate the Reviewers and Editors for their insightful comment. In the revised manuscript, we have now explicitly discussed the possible ways in which a more probabilistic approach could be integrated with our model, and our expectations for how our results might be altered in such an approach, as described above. We anticipate that this extended discussion will help to motivate further work that builds on the present manuscript.

2) The problem studied here is simplified to one spatial dimension. How might the conclusions change if there were variations orthogonal to the main front propagation direction? Could this lead to qualitatively different conclusions?

We thank the Reviewers and Editors for this thoughtful question; such orthogonal variations could give rise to rich new effects that will be interesting to study in future work. Indeed, we examined a similar question in our prior work modeling the collective migration of planktonic bacteria in the absence of quorum sensing-mediated biofilm formation. In that case, we found that variations orthogonal to the main front propagation direction “smooth out” as the front propagates, due to corresponding variations in the chemotactic response of the cells; cells at outward-bulging parts of the front are exposed to more nutrients, which diminishes their ability to respond to the nutrient gradient via chemotaxis and thus slows them down [Refs. 96 and 97 of our manuscript]. As a result, the migrating front eventually smooths to a flat shape whose subsequent dynamics can then be described using just one spatial dimension, just as in our treatment in the present manuscript.

However, we expect that this behavior could be altered in interesting new ways when the cells can additionally produce and sense autoinducer and thereby transition to the biofilm state, as is the case in the present manuscript. In this case, we speculate that because cells at outward-bulging parts of the front are exposed to more nutrients and have a weaker chemotactic response, autoinducer production and accumulation will be more rapid relative to cellular dispersal. That is, at these parts of the front, τα and τc will be shorter and longer, respectively, causing the dispersal parameter J~ to be larger locally. Thus, our model would predict biofilm formation to occur first at these parts of the front, potentially also influencing subsequent dispersal and biofilm formation at other locations along the front. Therefore, while the conclusions of our present manuscript could be the same locally at different parts of the front, the global behavior of the population could indeed be different—potentially giving rise to e.g., spatially-heterogeneous biofilm formation.

In the revised manuscript, we have now explicitly discussed possible ways in which variations orthogonal to the main front propagation direction could lead to fascinating new behaviors, as described above. We anticipate that this extended discussion will help to motivate further work exploring this direction of inquiry, and thank the Reviewers and Editors for encouraging us to think about these possibilities.

3) Please clarify the use of a single effective diffusivity parameter for a system that is slowly being arrested. That is, is the diffusion coefficient a constant in equation 1? The effective diffusivity of cells changes dramatically from the planktonic to the immobilized state. How is this issue resolved in the model?

We apologize for not adequately clarifying how the cellular diffusivity changes from the planktonic to the immobilized state in the previous version of the manuscript, and appreciate having the chance to do so here. In our model, the diffusivity of cells in the planktonic state (in Equation 1) is indeed a constant D1. When these cells encounter sufficiently concentrated autoinducer (a ≥ a), they transition to the immotile biofilm state after a time lag ~τ.

The Reviewers and Editors astutely point out that in real systems, the change in cellular diffusivity (and chemotactic coefficient) may not be as temporally abrupt—although we are not aware of specific characterization of these dynamics, which would be useful for future experiments to probe. Thus, in the absence of such characterization, for simplicity, we treat this transition as being step-like from a single constant value of diffusivity D1 and a single constant value of the chemotactic coefficient χ1 in the planktonic state to zero diffusivity and chemotactic coefficient in the biofilm state after the time lag ~τ. Future work could, for example, incorporate a temporally-varying diffusivity (and chemotactic coefficient) into our model that slowly transitions from the planktonic values D1 and χ1 to the biofilm value of zero over a non-zero time scale. Given that the same cells would be transitioning from the motile planktonic to immotile biofilm state—but in this case with the introduction of a time-varying diffusivity and chemotactic coefficient—we expect that the long-time biofilm fraction f will be similar, and only the spatial profile of the biofilm population may be altered. Hence, we expect that our main findings summarized in Figure 4 will be unaffected. To test this expectation, we have now run a new version of the identical representative simulation shown in Figure 1C, but with both motility parameters D1 and χ1 smoothly transitioning to zero as described above. We observe nearly-identical results for both cases, confirming our expectation that the temporal nature of the arrest in motility does not appreciably influence our model results and conclusions.

We thank the Reviewers and Editors for asking this interesting question. In the revised manuscript, we have now explicitly discussed possible ways in which these dynamics could be incorporated into our model, as described above, along with the new results presented above. Further investigating these dynamics will be an interesting direction to explore in future work.

4) Biofilms are usually formed by different types of micro-organisms, and I am struggling to see whether this single species model is general enough to capture that reality. We do not mean to ask for the authors to perform more simulations, but we would like to hear know their thoughts on this issue.

5) While we believe that the results are compelling and quite useful, the work can be strengthened. For example, biofilms are usually heterogeneous (i.e., different microbial populations) and many of the parameters in the model are homogeneous. It would be interesting to understand how even a second cell species would alter (or not) the results, in particular the claim that onset of biofilm formation can be uniquely defined by the parameters D and J.

We are responding to both questions 4 and 5 together, because they are closely related. We are grateful to the Reviewers and Editors for encouraging us to clarify possible ways our model could be extended to consider variability in cellular behaviors in a population, as well as the inclusion of multiple species in a biofilm community. Indeed, we believe our theoretical framework provides a useful foundation for modeling these important features of real-life systems.

As the Reviewers and Editors astutely pointed out, biofilms are often formed by multiple different microbial species, which we do not consider in this paper. Instead, as a first step, our work describes biofilm formation by a single species, for simplicity. Nevertheless, we expect that our theoretical framework can be extended by following reasoning similar to that described in our paper, but with the introduction of additional equations and variables in the governing Equations (1)— (4) to describe the distinct cell and chemical types. For example, in the case of biofilm formation by two different species, A and B, under positive quorum sensing control:

i) In the case that the different species consume and respond to distinct nutrients ci (where i ∈ {A,B} indexes each species), and secrete and respond to distinct autoinducers ai, each species can be described in isolation using our same governing Equations (1)—(4), but now extended to incorporate the distinct variables ci, a ,b1,i, and b2,i-. Then, directly following the approach described in our paper, each species would be described by its own dimensionless parameters Di~  Td, i/Ta, i and Ji~  Tc, i/Ta, i, with D~/D~ +J~/J~   1 again specifying the threshold for biofilm formation for each. We hypothesize that the composition of the final two-species biofilm community would then be given by the combination of each single-species biofilm.

ii) In the case that the different species consume and respond to the same nutrient c, and secrete and respond to the same autoinducer +, our governing Equations (1)—(4) could be extended to consider the cellular parameters specific to each species i. In this approach, however, biofilm formation by each of the two species cannot be described in isolation, because they are coupled through the nutrient and autoinducer dynamics. For example, the calculation of the nutrient depletion time scale τd would need to be extended—following the same derivation as in our paper—to now reflect the aggregate nutrient consumption by both species. Similarly, the time scales τa and τc would be extended, following our paper, to now reflect contributions from both species. Then, just as in our paper, the overall two-species community would be described by one set of governing dimensionless parameters (D~,J~) with D~/D~ +J~/J~~1 specifying the threshold for biofilm formation for the overall twospecies community. We hypothesize that while this relation would specify a universal biophysical threshold for the onset of biofilm formation by the entire community, the composition of the final two-species biofilm that results above this threshold may not be uniquely specified by (D~,J~).

The above serve as two possible examples of how our model could be extended for the case of two-species communities. Similar extensions can be made as needed for e.g., communities with even more species, or communities with species that transition to/from the biofilm state in different ways (e.g., reversible biofilm formation, biofilm formation under negative quorum sensing control or regulated by other non-quorum sensing-based mechanisms).

We thank the Reviewers and Editors for encouraging us to think about the possible ways our work could be extended to describe additional complexities of biofilm formation in real-life settings. In the revised manuscript, we have now explicitly discussed the simplifications made in our model, its limitations, and the possible ways in which it could be extended to describe different cell types and mechanisms of chemotaxis/biofilm formation. We anticipate that our model will provide a foundation for other researchers to theoretically describe these and other complexities.

6) Are these results sensitive to the chosen initial condition? That is, it would be good to explore the dynamics of the biophysical model (perhaps that has been done already, but we did not see it in the manuscript). For example, are there oscillatory solutions or closed orbits, or even singularities that arise from the solution of the equations?

We thank the Reviewers and Editors for these interesting questions, and appreciate having the chance to clarify the details of our model and simulations. Because the goal of this present manuscript is to present the biophysical model and examine the predictions it makes for the onset of biofilm formation, our work thus far used numerical solutions of the full system of equations. So, we have not yet fully probed the mathematical structure of the governing equations to examine the possibility of e.g., oscillatory solutions, closed orbits, or singularities—which would be fascinating to explore in the future. We expect that oscillatory solutions/closed orbits are not likely to arise, given that our model takes biofilm formation to be irreversible; however, such complex long-time dynamics may emerge if the transition to biofilm formation were taken to be reversible. More broadly, exploring our system of equations within the framework of dynamical systems theory will be an interesting direction for future work.

However, to assess the sensitivity of our results to the chosen initial condition, we performed a new simulation identical to that shown in Figure 1C, but with the shape of the initial inoculum changed from a Gaussian profile to a step function (with the same maximum cellular concentration and width). We observe nearly-identical results for both cases, indicating that the results are robust to variations in this initial condition chosen. The influence of other changes in the initial conditions (e.g., maximal bacterial concentration, nutrient concentration) is already described by our model through the time scales τd,τa, and τc that define the governing dimensionless parameters (D~,J~).

https://doi.org/10.7554/eLife.76380.sa2

Article and author information

Author details

  1. Jenna A Moore-Ott

    Department of Chemical and Biological Engineering, Princeton University, Princeton, United States
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6832-0658
  2. Selena Chiu

    Department of Chemical and Biological Engineering, Princeton University, Princeton, United States
    Contribution
    Data curation, Investigation, Software
    Competing interests
    No competing interests declared
  3. Daniel B Amchin

    Department of Chemical and Biological Engineering, Princeton University, Princeton, United States
    Contribution
    Methodology, Software
    Competing interests
    No competing interests declared
  4. Tapomoy Bhattacharjee

    Andlinger Center for Energy and the Environment, Princeton University, Princeton, United States
    Contribution
    Methodology
    Competing interests
    No competing interests declared
  5. Sujit S Datta

    Department of Chemical and Biological Engineering, Princeton University, Princeton, United States
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Visualization, Writing – original draft, Writing – review and editing
    For correspondence
    ssdatta@princeton.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2400-1561

Funding

National Science Foundation (CBET-1941716)

  • Sujit S Datta

National Science Foundation (EF-2124863)

  • Sujit S Datta

National Science Foundation (DMR-2011750)

  • Sujit S Datta

Pew Charitable Trusts (Pew Biomedical Scholars Program)

  • Sujit S Datta

National Science Foundation (DGE-1656466)

  • Jenna A Moore-Ott

Princeton University (Eric and Wendy Schmidt Transformative Technology Fund)

  • Sujit S Datta

Princeton University (Princeton Catalysis Initiative)

  • Sujit S Datta

Princeton University (Reiner G. Stoll Undergraduate Summer Fellowship)

  • Selena Chiu

Princeton University (Princeton University Library Open Access Fund)

  • Jenna A Moore-Ott

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

It is a pleasure to acknowledge R Kōnane Bay, Vernita Gordon, Anna M Hancock, Andrej Košmrlj, Alejandro Martínez-Calvo, JTN Moore-Ott, PG-A Moore-Ott, NA Moore-Ott, Howard Stone, Carolina Trenado-Yuste, and Ned Wingreen for stimulating discussions. This work was supported by NSF grants CBET-1941716 and EF-2124863, the Eric and Wendy Schmidt Transformative Technology Fund at Princeton, the Princeton Catalysis Initiative, a Reiner G Stoll Undergraduate Summer Fellowship (to SC), in part by funding from the Princeton Center for Complex Materials, a Materials Research Science and Engineering Center supported by NSF grant DMR-2011750, and the Pew Charitable Trusts through the Pew Biomedical Scholars Program. This material is also based upon work supported by the National Science Foundation Graduate Research Fellowship Program (to JAM-O) under Grant No. DGE-1656466. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. This publication was supported by the Princeton University Library Open Access Fund.

Senior Editor

  1. Aleksandra M Walczak, CNRS LPENS, France

Reviewing Editor

  1. Raymond E Goldstein, University of Cambridge, United Kingdom

Publication history

  1. Preprint posted: December 5, 2021 (view preprint)
  2. Received: December 14, 2021
  3. Accepted: June 1, 2022
  4. Accepted Manuscript published: June 1, 2022 (version 1)
  5. Version of Record published: July 21, 2022 (version 2)

Copyright

© 2022, Moore-Ott et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,331
    Page views
  • 543
    Downloads
  • 5
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Jenna A Moore-Ott
  2. Selena Chiu
  3. Daniel B Amchin
  4. Tapomoy Bhattacharjee
  5. Sujit S Datta
(2022)
A biophysical threshold for biofilm formation
eLife 11:e76380.
https://doi.org/10.7554/eLife.76380
  1. Further reading

Further reading

    1. Computational and Systems Biology
    2. Physics of Living Systems
    Nil Z Gurel, Koustubh B Sudarshan ... Olujimi A Ajijola
    Research Article

    Stellate ganglia within the intrathoracic cardiac control system receive and integrate central, peripheral, and cardiopulmonary information to produce postganglionic cardiac sympathetic inputs. Pathological anatomical and structural remodeling occurs within the neurons of the stellate ganglion (SG) in the setting of heart failure. A large proportion of SG neurons function as interneurons whose networking capabilities are largely unknown. Current therapies are limited to targeting sympathetic activity at the cardiac level or surgical interventions such as stellectomy, to treat heart failure. Future therapies that target the stellate ganglion will require understanding of their networking capabilities to modify any pathological remodeling. We observe SG networking by examining cofluctuation and specificity of SG networked activity to cardiac cycle phases. We investigate network processing of cardiopulmonary transduction by SG neuronal populations in porcine with chronic pacing-induced heart failure and control subjects during extended in-vivo extracellular microelectrode recordings. We find that information processing and cardiac control in chronic heart failure by the SG, relative to controls, exhibits: i) more frequent, short-lived, high magnitude cofluctuations, ii) greater variation in neural specificity to cardiac cycles, and iii) neural network activity and cardiac control linkage that depends on disease state and cofluctuation magnitude.

    1. Computational and Systems Biology
    2. Neuroscience
    Sergio Oscar Verduzco-Flores, Erik De Schutter
    Research Article Updated

    How dynamic interactions between nervous system regions in mammals performs online motor control remains an unsolved problem. In this paper, we show that feedback control is a simple, yet powerful way to understand the neural dynamics of sensorimotor control. We make our case using a minimal model comprising spinal cord, sensory and motor cortex, coupled by long connections that are plastic. It succeeds in learning how to perform reaching movements of a planar arm with 6 muscles in several directions from scratch. The model satisfies biological plausibility constraints, like neural implementation, transmission delays, local synaptic learning and continuous online learning. Using differential Hebbian plasticity the model can go from motor babbling to reaching arbitrary targets in less than 10 min of in silico time. Moreover, independently of the learning mechanism, properly configured feedback control has many emergent properties: neural populations in motor cortex show directional tuning and oscillatory dynamics, the spinal cord creates convergent force fields that add linearly, and movements are ataxic (as in a motor system without a cerebellum).