Differential equation-based minimal model describing metabolic oscillations in Bacillus subtilis biofilms

Biofilms offer an excellent example of ecological interaction among bacteria. Temporal and spatial oscillations in biofilms are an emerging topic. In this paper, we describe the metabolic oscillations in Bacillus subtilis biofilms by applying the smallest theoretical chemical reaction system showing Hopf bifurcation proposed by Wilhelm and Heinrich in 1995. The system involves three differential equations and a single bilinear term. We specifically select parameters that are suitable for the biological scenario of biofilm oscillations. We perform computer simulations and a detailed analysis of the system including bifurcation analysis and quasi-steady-state approximation. We also discuss the feedback structure of the system and the correspondence of the simulations to biological observations. Our theoretical work suggests potential scenarios about the oscillatory behaviour of biofilms and also serves as an application of a previously described chemical oscillator to a biological system.


Introduction
social dynamics come into play [2,3]. One of these is the division of labour [4,5]. The core of the biofilm growing on a solid surface shows a different metabolic state than the periphery. The periphery can freely access the nutrients from the surrounding environment. The interior, however, faces hindrance in obtaining a stable inflow of nutrients because the peripheral cells use up the nutrients that diffuse towards the interior. An experimental set-up to simulate that situation is provided by a microfluidics chamber [4].
An example of such a nutrient gradient is the production and diffusion of ammonia in the biofilm. Every cell in the biofilm has the ability to produce ammonia [4,6]. However, this small chemical compound is highly diffusive and therefore escapes into the environment as soon as it is produced by the cells in the periphery, thus leading to waste of nitrogen. In the interior, the ammonia produced by the cells diffuses out into the periphery. Thus, the interior cells monopolize ammonia production for the entire biofilm. Ammonia being an essential component of glutamine metabolism could be used to control the growth rate of the periphery by limiting its supply. The interplay between the inner and outer cells is required for glutamine synthesis and therefore the growth of the biofilm [4,6,7].
To understand biofilms more closely and make predictions based on empirical data, several models have been developed [4,[7][8][9][10][11][12]. Liu et al. [4] observed oscillations in the biofilm, which they explained by different metabolic roles performed by the different compartments in the biofilm. They also established a model based on six differential equations. They defined two regions: the interior and periphery. Each of the regions has a variable representing glutamate and another representing the concentration of housekeeping proteins like ribosomal proteins. Ammonia and the active form of the enzyme glutamate dehydrogenase are also variables of the model.
Since many biological oscillators have been described by less than six variables [13,14], a simpler model could be established for biofilm oscillations as well. Our ultimate aim was to develop a minimal model to describe the metabolic oscillations happening in a biofilm. Minimal models are the simplest way to describe a certain phenomenon with the least number of parameters [15] and this is in agreement with Occam's razor. For example, minimal models were established for glycolytic oscillations by Higgins [16] and Sel'kov [17] and for calcium oscillations by Somogyi & Stucki [18].
Here, we employ the smallest chemical reaction system showing a Hopf bifurcation [19], which was further analysed [20,21] and used to describe p53 oscillations [22]. At a Hopf bifurcation, damped oscillations turn into undamped oscillations [15,23]. In particular, Wilhelm & Heinrich [19,20] performed a thorough stability analysis of the model. We test to what extent the terms in this model match the processes in a biofilm system. In this analysis, we focus on the Hopf bifurcation, discuss the feedback structure and point out the correspondence of the simulations to biological observations.
In our model, we use three variables only: ammonia, and interior and peripheral glutamate. Besides the quest for minimality, a reason for not considering the concentrations of housekeeping proteins as variables is that they change on a longer time-scale than metabolites. A similarity to the larger Liu model [4] is that, among the various amino acids, we focus on the metabolism of glutamate since glutamate and ammonia are both involved in the production of various amino acids through trans-amination, which is then equated to growth.
To study the effect and possible benefit of oscillations, it is of interest to compute the average values of variables, as was done for several oscillators [24][25][26][27][28][29]. For linear differential equation systems showing oscillations (such as the system describing the harmonic pendulum), the average values equal the values at the marginally stable steady state. For nonlinear differential equation systems, the average values often differ from the values at the unstable steady state surrounded by the oscillations. However, there are some types of nonlinear systems for which equality holds, for example, Lotka-Volterra systems of any dimension [25]. The equality property has also been proved for some models of calcium oscillations [27,30] and the Higgins-Selkov oscillator [27]. Here, we probe the model employed for describing biofilm oscillations for the above-mentioned property.

The model
Based on the scenario described by Liu et al. [4], the biofilm was separated into two compartments-the interior and the periphery. Here, we use a minimalist approach and try to model biofilm oscillations with the simplest model possible. Accordingly, we use the smallest chemical reaction system with Hopf royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 190810 bifurcation [19]. The term chemical system mathematically means that only up to bilinear terms are involved. It turns out that this model matches the biological set-up.
The model includes five reactions (figure 1) and three species with variable concentrations. The general variables X, Y and Z from the Wilhelm and Heinrich model can be assigned for the biofilm system to: peripheral glutamate (G p ), ammonia (A) and internal glutamate (G i ), respectively. Based on mass-action kinetics, the reactions have been translated as follows into a system of ordinary differential equations (ODEs): Model assumptions and interpretation of terms in the model: The uptake of glutamate from the environment (G E ) by the periphery of the biofilm. G E is supplied in a large excess, hence considered constant. The uptake of glutamate (G p ) is dependent on itself because glutamate represents the total amino acid and thus protein concentration in the biofilm periphery and can be assumed, in rough approximation, to be proportional to the concentration of various transport proteins embedded in the cell membranes. The greater the concentration of these proteins, the higher is the glutamate uptake rate. Without this self-amplification of glutamate, the system would not oscillate by construction of the minimal model. k 4 G p : Diffusion of glutamate from the periphery of the biofilm into its interior. We do not consider self-amplification by G i in the main text. We analysed the effect of self-amplification of G i using the term k 4 G i G p (electronic supplementary material, figure S4). k 2 AG p : Consumption of glutamate and ammonia to produce biomass. As a simplification, we assumed that only the interior cells produce ammonia since that produced by the peripheral cells is rapidly lost to the environment. k 5 G i : Consumption of glutamate to produce ammonia.
k 3 A: Diffusion of ammonia into the surroundings. The loss of ammonia due to diffusion is much larger than that taken up by the periphery to produce biomass. Therefore, the term k 2 AG p does not appear in equation (2.1b).

Simulation
For computer simulations, we used the software COPASI v. 4.16 and 4.24 [31] and its LSODA deterministic solver. The simulations were double-checked using the Matlab ode15s (MathWorks) function. The figures of the simulations were produced using COPASI, and the three-dimensional (3D) phase plot was generated using the lines3D function of R plot3D library. The biomass plot was generated using the R function ggplot. Parameter values are given in table 1. They are obtained by rescaling the parameter values from the Wilhelm & Heinrich paper [19] such that the oscillation period observed in the experimental work by Liu et al. [4] is matched. The glutamate concentration in the environment was adopted from the Liu et al. paper [4]. We have chosen the rate constant of diffusion of ammonia, k 3 , to be twice as high as that of glutamate, k 4 . This is because the diffusion coefficient for ammonia [32] is about 1.6 × 10 −5 cm 2 s −1 , while that for glutamate [33,34] is about 8 × 10 −6 cm 2 s −1 . k 1 and k 2 had to be increased in order to obtain undamped oscillations and to match the same period. Overall, the parameters allow a good comparison to the results by Wilhelm & Heinrich [19], while also being realistic from a physicochemical point of view.
The predicted doubling time was calculated by averaging the relative increase in biomass at four consecutive time points of the maxima of ammonia concentration.

Steady states
The steady states of the system can be calculated analytically. This gives a trivial steady state (TSS) and a non-trivial steady state (NTSS) It is worth noting that the concentrations at the latter state are linear functions of G E . The TSS and NTSS are stable if k 1 G Ek 4 is negative or positive, respectively [19]. At the threshold, a transcritical bifurcation [35] occurs; that is, the two steady states interchange their stability. At a further threshold, k 1 G E = k 3 + k 4 + k 5 , the stable NTSS turns unstable in a Hopf bifurcation [19].

Time course shows oscillations
We run the time course calculation of the system (2.1a-c) for 25 simulation hours with 1000 steps each of size 0.025 h (1.5 min). The period for oscillations is about 126 min (2 h 6 min), which is in agreement with the experimental observations [4], because the parameters have been rescaled accordingly (see above). The amplitude of oscillations is observed to be 3.0 mmol l −1 for ammonia, 7.1 mmol l −1 for interior glutamate and 14.9 mmol l −1 for peripheral glutamate (figure 2). It can be seen that the three variables oscillate with phase shifts, i.e. asynchronously. In order to see the interdependence between the variables of our model, we plot the phase portrait of all three variables for various values of G E ( figure 3). G E is an appropriate bifurcation parameter because   As per the assumptions of our model, k 2 AG p is a proxy for the input for the synthesis of biomass from ammonia and glutamate and is thus related to the growth of the biofilm. Biomass production can be described by the following differential equation: where b is a conversion factor and was tuned to 0.1 ((mmol l −1 ) 2 h) −1 so that the doubling time is in agreement with the experimental values [36]. The number of cells can be converted to biomass by taking the volume of a typical B. subtils cell of about 0.85 ± 0.38 µm 3 [37] and the average density of 1 g ml −1 which results in 8.5 ± 0.38 × 10 −13 g cell −1 . The numerical solution of equation (3.3) for various initial values of G p is shown in figure 4. It can be seen that there is periodic retardation in growth. Figure 4 also displays the growth curve in the hypothetical case where G p and A subsisted at steady state (black curve).
The initial value of G p for the growth with constant growth rate (black monotonic curve) was chosen such that biomass is comparable to that for oscillating growth in the first 10 h. If the same initial values as for the growth with varying growth rate were chosen, biomass would grow to higher values right from the beginning. Thus, the numerical calculations suggest that oscillating growth for this system is not in favour of increasing growth rate. As can be seen from electronic supplementary material, figure S5, the steady-state growth rate overtakes the oscillating growth rate at about 10.5 h.
It is an important result that biofilm oscillations can be described by considering a few processes only, which are listed below equation (2.1). They include considerably less processes than the Liu model [4]. However, they do include the diffusion of ammonia to the surroundings, unlike that model. Thus, it is plausible to assume that these are the most relevant processes for the phenomenon of oscillations.  All other processes (such as the diffusion of glutamate from the interior of the biofilm to its periphery and from there to the environment) can be neglected.

Average concentrations and average growth rate
We integrate over one oscillation period, Now, we calculate the integral of dA=dt: Now, we calculate the integral of dG i =dt and derive, analogously Thus, the average concentrations are equal to the steady-state concentrations. The question arises whether the oscillations have an effect on the average of the bilinear term k 2 AG p . This is not immediately clear, although figure 4 suggests that the average growth rate is slower than that at the metabolic steady state. Note that the growth term bAG p B is trilinear.
To check whether hA G p i ¼ A ss : G p,ss , we integrate dG p =dt over one period, Since hG p i ¼ G p,ss and hAi ¼ A ss ¼ k k 2 , dividing by k 2 and T gives hA G p i ¼ A ss : G p,ss : Thus, the average of the bilinear term, which can be interpreted as the input to biomass, is indeed unaffected by oscillations, although the ammonia and peripheral glutamate levels oscillate asynchronously. For two-dimensional Lotka-Volterra systems, this property was shown earlier [24]. These values can also be calculated from the general formulae for the bifurcations [19]. The steady-state value of G p is a linear function of the bifurcation parameter G E , as shown in equation (3.2a). It can be seen that the Hopf bifurcation is supercritical; that is, the amplitude grows gradually starting from zero and the limit cycle is stable right from the beginning. Near the Hopf bifurcation, the obtained time course curve ( figure 2) is sinusoidal. For G E ≫ 24.41 mmol l −1 the oscillations get spike-like and are no longer sinusoidal. It is of interest to speculate about the physiological advantage of spike-like oscillations. This question has been discussed earlier in the context of calcium oscillations [26,38,39]. Whenever the kinetic effect of the oscillating variable (e.g. in activating a protein or in a biochemical conversion) is nonlinear and follows a convex function, the spikes contribute more than proportionately to the effect. Thus, spike-like oscillations can lead to a high average effect even at low average value of the variable. In order that oscillations really enable division of labour in the case of biofilms, it can be expected that they should not be sinusoidal. This deserves further studies. A biological explanation of the bifurcations is given in the Discussion.

Bifurcations
We checked the parameter sensitivity in two ways. First, we performed a bifurcation analysis for all parameters (except k 1 since it is equivalent to G E ), see figures 5 and 6 and electronic supplementary material, figures S1-S3. We see a steep increase in the oscillation amplitude with respect to G E and k 3 , a moderately steep increase for k 4 , whereas the other parameters show a very gradual increase in the vicinity of the bifurcation. Second, we applied local parameter sensitivity analysis to the steadystate concentrations, which are equal to the average values. This can be done in an analytical way by differentiating the steady-state values given in equation (3.2) consecutively with respect to all parameters. The resulting unscaled sensitivities for the parameter values from table 1 are given in electronic supplementary material, table S1a. Thereafter, we computed the scaled sensitivities by multiplying by the parameter and dividing by the concentration (electronic supplementary material, table S1b). The obtained values for all sensitivities were confirmed by numerical computation using COPASI [40].
The results show that A is not sensitive at all to k 3 and k 5 , nor is G p to k 5 . This is counterintuitive because increasing k 5 corresponds to over-expression of glutamate dehydrogenase, which produces ammonia. However, in our model, increasing k 5 leads to a decrease in G i , so that the term k 5 G i stays constant. It deserves further study whether a more realistic kinetics leads to non-zero sensitivity for these cases.
Scaled sensitivities are equal to unity if the parameter enters the formula for the steady value in a multiplicative way (that is, as a factor in the numerator) and equal to minus unity if it is a factor in the denominator. Examples are provided by the scaled sensitivities of G i and G p with respect to k 3 and of all concentrations with respect to k 2 , respectively. The highest scaled sensitivity, 1.24, is found for all concentrations with respect to k 1 . That value means that a 1% increase in a parameter value leads to a 1.24% increase in the concentration.

Quasi-steady-state approximation
To ascertain the cause of oscillations, which could be a negative feedback with a delay or a positive feedback, we can study a subsystem by eliminating a variable. This can be done by the quasi-steadystate approximation (QSSA) [15]. In our system, we see that G p exerts a positive feedback on itself which is linear and, thus, quite weak. For example, in the Higgins-Selkov oscillator involving two variables, the feedback is quadratic [16,17]. Moreover, the above system involves a negative feedback: G p is converted to G i , G i is converted to A and A promotes the degradation of G p (figure 1). In the Goodwin oscillator, which also consists of three variables, a negative feedback is the cause of oscillation [41,42]. Inspired by the observation in figure 5 that the oscillations vanish at high k 5 values, we apply the QSSA for G i . This corresponds to the special case where glutamate dehydrogenase is overexpressed. In that case, indeed quenching of oscillations was observed in experiment and also predicted by the Liu model [4]. The variable G i then attains a quasi-steady state This leads to Substituting this into the above model equations (2.1a-c) yields a simplified system and dA dt ¼ Àk 3 A þ k 4 G p : ð3:6bÞ It is of interest to analyse its dynamics. It follows from the general result by Hanusse [43,44] that it cannot give rise to a limit cycle because it involves two variables and only linear and bilinear terms. However, the question still remains whether it gives rise to a stable or unstable steady state, whether damped oscillations are possible, etc. System (3.6) shows two steady states which is the TSS, and ð3:8a,bÞ which is the NTSS (see equations (3. 2a,b)). The Jacobian matrix reads while for the TSS, it reads For matrices with such a triangular structure, the eigenvalues are given by the diagonal elements. In our case ð3:11Þ In any case, the eigenvalues are real, so that not even damped oscillations are possible. For k 1 G E < k 4 , both eigenvalues are negative, so that the TSS is a stable node. For k 1 G E > k 4 , one eigenvalue is negative and the other one positive. The steady state then is unstable, it is a saddle point. For the NTSS (3.8), the Jacobian matrix becomes The characteristic equation reads This has the solutions Now, we distinguish three cases: (a) For k 1 G E < k 4 , the term under the square root is positive, so that the root is real. Moreover, it is larger than k 3 /2. Thus, one eigenvalue is negative and the other one positive. The steady state then is unstable, it is a saddle point. (b) For 0 < k 1 G E − k 4 < k 3 /4, the root is again real. It is less than k 5 /2, though. Both eigenvalues are negative; the steady state is a stable node. (c) For k 1 G E − k 4 > k 3 /4, the root is imaginary. Both eigenvalues are complex numbers, with the same negative real part −k 3 /2. The steady state is a stable focus. This state is, thus, reached by damped oscillations.
From these calculations, the following conclusions can be drawn. At k 1 G E = k 4 , the two steady states of the simplified system (3.6) coincide, as in the complete system (2.1). Since the TSS and NTSS interchange their stability at that point, it is a transcritical bifurcation. There is a second transition point where the qualitative behaviour changes, at k 1 G E − k 4 = k 3 /4. This is another point than the Hopf bifurcation in the complete system, which is at k 1 G E − k 4 = k 3 + k 5 . At this transition in the simplified system, a stable node turns into a stable focus. Such a transition must also occur in the complete system between the transcritical bifurcation, where an unstable node turns into a stable node, and the Hopf bifurcation, where a stable focus gets unstable. It is difficult to find it exactly in the three-dimensional system. Beyond this transition, the simplified system shows damped oscillations. This implies that the positive feedback of G p on itself can be considered as a cause of oscillation, yet not as a cause of a limit cycle. The 'inflation' of the oscillation to a limit cycle in the royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 190810 complete system appears to be brought about by the negative feedback via G i . The loop via G i can, moreover, be interpreted as a delay.
In this subsection, we have considered the special case of high k 5 . This corresponds to a situation realized in experiments by Liu et al. where they overexpress the glutamate dehydrogenase leading to an excessive production of ammonia [4]. In that situation, indeed no oscillations were observed. We have proved analytically that the limit cycle disappears in that case. This generalizes the numerical finding shown in figure 6. Thus, to model limit cycle oscillations in biofilms by equation (2.1), we need the full three-dimensional system with values of k 5 that are not too high.
The QSSA for the rate constant k 3 is given in the electronic supplementary material.

Discussion
Here, we have used the smallest chemical system showing a Hopf bifurcation to model metabolic oscillations in B. subtilis biofilms. That model had been used earlier to describe p53 oscillations [22]. Here, we have applied the model to describe another biological phenomenon. We have specifically selected the parameter values to describe biofilm dynamics, which makes the model more relevant in the light of the biological observation. In our system, the diffusion of ammonia is critical for biofilm oscillations. All the terms in the model are linear, except k 2 AG p , which is bilinear. The model describes metabolic and diffusion processes as outlined above. As an output, the growth of the biofilm (consisting of incremental and halting phases) was also computed (figure 4). A major reason of the observed oscillations was demonstrated to be the division of labour between the central and peripheral zones of the biofilm. While the release of usable ammonia is mainly delimited to the former, the production of biomass and, thus, growth, is mainly delimited to the latter.
We have presented bifurcation diagrams, which clearly show supercritical Hopf bifurcations (figure 5 and 6, electronic supplementary material, figures S1-S3). Earlier, Wilhelm & Heinrich [19,20] had analysed that bifurcation and had presented one bifurcation diagram. Here, we have added some mathematical analysis. For example, we show the maxima and minima of oscillations, the knowledge of which gives us a quantitative insight into the biofilm dynamics. Moreover, we performed a QSSA and probed for the equality property of the average values. We analysed the Hopf bifurcation by changing not only the external glutamate G E but alternatively also all the rate constants except k 1 , since changing k 1 has the same effect as changing G E . Interestingly, a recent model [7] has shown a subcritical bifurcation in describing the behaviour of the stress levels in the biofilm periphery. However, they modelled the stress with a single delay differential equation and did not consider other molecular details, while we do not consider stress. Using a single differential equation meets the quest for minimality. However, a delay differential equation (meaning that the time derivative of a variable depends on that variable at a previous time point) is, from a mathematical point of view, very complex because it requires infinitely many initial values (from zero to the delay period, with a simplifying assumption being that they are all equal). Moreover, stability analysis is then considerably more complicated. Our model is complementary to their model. It is closer to Liu's original model [4] but much simpler because it involves only three rather than six variables, and thus only requires three initial values. We chose the parameters of the model such that they are in agreement with Liu's experimental results, namely the period and the amplitude of oscillations.
In our model, peripheral glutamate exerts a positive feedback on itself. Mathematically, this has the form of a bilinear term involving peripheral and external glutamate concentrations. At very low values of G E , the feedback is not strong enough to enable a positive steady state. The system then tends to the TSS. In that state, G p is zero, so that growth is impossible. Biologically, this can be interpreted that the biofilm is too small to be viable. This is in agreement with observations in the recent study from the Suel group [7]. At a certain threshold value of G E (5.88 mmol l −1 ), the NTSS turns stable in a transcritical bifurcation. Beyond that value, the feedback is strong enough to enable growth. At high values of G E , the feedback becomes so strong that an overshoot occurs: more glutamate is taken up than needed, so that the G p level transiently exceeds the steady-state value. Then, more peripheral glutamate is consumed for release of ammonia or for growth, so that the concentration decreases again. This leads to oscillations.
From a functional point of view, a steady state is quite appropriate [29]. Growth of the biofilm does not require oscillations. However, in this system, oscillations help in mitigating the chemical attack that challenges the biofilm [4]. This may have interesting clinical implications in view of treatment of biofilmforming bacteria by antibiotics. Furthermore, another study [9] indicates that oscillations in growth royalsocietypublishing.org/journal/rsos R. Soc. open sci. 7: 190810 actually help in sharing the nutrients among several biofilms more efficiently. However, not all biofilms show oscillations, indicating that it is not critical for biofilms. Our numerical calculations suggest that the average growth rate is lower when compared with growth at the metabolic steady state.
By contrast, the individual concentration variables (ammonia, peripheral glutamate etc. but not biomass) show the equality property that their average during oscillations equals the value at the unstable state, as usual for Lotka-Volterra systems [25]. Here, we have shown that the bilinear input term k 2 AG p exhibits this equality property as well. This may come as a surprise because the ammonia and peripheral glutamate levels oscillate with a phase shift.
In the paper by Liu et al. [4], the oscillations computed by their model have a sinusoidal shape. In our model, such a shape only occurs in the neighbourhood of the Hopf bifurcation. Further beyond it, the shape is more spike-like with the crests being sharper than the troughs.
The question arises whether the model used and analysed here is minimal. On the basis of ODE systems (without delays), at least two variables are needed to generate oscillations [13,14]. However, when only linear and bilinear terms are included, at least three variables are needed, as was proved by Hanusse by an analysis of the Jacobian matrix [43,44]. As shown by Wilhelm & Heinrich [19], such a model requires at least five reactions. Thus, the above model is minimal in terms of the number of variables (criterion with highest priority) and number of reactions, given the type of kinetics. The famous two-variable Brusselator model [45] and the Higgins-Selkov oscillator involve a term of degree three each [16,17]. If the number of reactions is granted the highest priority, the model may look different. Thus, it depends on the criteria what a minimal model is. Note that a delay differential equation [7] is, from the viewpoint of the number of initial values, of infinite dimension. While our model is not necessarily the simplest, it provides a trade-off between simplicity and adequacy to match the observed oscillation in biofilms.
As for any oscillatory system, it is interesting to elucidate the feedback structure. The term k 1 G E G p represents a positive feedback because peripheral glutamate stimulates its own uptake. This is because glutamate is a proxy for the concentration of various transport proteins embedded in the cell membranes. The higher the concentration of these proteins, the higher is the glutamate uptake rate. Since this positive feedback is the driving force for oscillations, at low values of k 1 G E , we observe a steady state rather than oscillations. In glycolytic and calcium oscillations, the cause of oscillations is also a positive feedback [13,14,17,46], while in a Goodwin oscillator, it is a negative feedback [41,42].
In addition to the positive feedback, there is also a negative feedback loop in the system (figure 1). As seen from the differential equations, peripheral glutamate positively influences internal glutamate, which positively affects ammonia, which then negatively influences peripheral glutamate. Thus, the overall effect is inhibitory. This feedback structure of the Wilhelm-Heinrich model has been highlighted earlier [22].
By applying the QSSA, we have proved analytically that the limit cycle disappears if glutamate is degraded very fast or ammonia diffuses very easily. As mentioned in the Results section, the former case corresponds to a situation realized in experiments by overexpressing the glutamate dehydrogenase [4]. In that situation, no oscillations were observed. By contrast, in the case where oscillations occur, a description by a simple mass-action system requires three variables. Analyses in this direction may be relevant for clinical interventions via inhibition or activation of bacterial enzymes or changing diffusivity in the biofilm.
The model analysed here has several pros and cons. In view of the mathematical analysis, its simplicity is certainly an advantage. In view of an adequate description of the biological and biochemical processes involved, the model may appear oversimplified. For example, describing growth by trilinear term is quite simplistic; usually, it is described by saturation kinetics (e.g. Michaelis-Menten kinetics). In addition, the assumption that glutamate uptake by the periphery is proportional to the glutamate concentration as explained above could be refined in future studies. Moreover, diffusion processes are usually reversible. In the above model, we neglected the backward processes in diffusion, which is justified if concentration differences are high.
Many theoretical and experimental studies have been published on glycolytic oscillations [13,16,17,47,48]. However, these oscillations only occur under very special or even artificial conditions. In living cells, metabolic oscillations are rare, while being quite frequent in signalling systems [13,14,49,50]. The lights of a car are a helpful analogy: the headlights illuminate the street in a permanent way; there is no point in oscillations. By contrast, the side indicators (as signalling device) flash; that is, they emit oscillating light. Interestingly, in the case of biofilms, metabolic oscillations could provide advantages. While the work described here is quite theoretical, we consider it to be an appropriate basis for refined and more sophisticated models of biofilm oscillations.