Extending the minimal model of metabolic oscillations in Bacillus subtilis biofilms

Biofilms are composed of microorganisms attached to a solid surface or floating on top of a liquid surface. They pose challenges in the field of medicine but can also have useful applications in industry. Regulation of biofilm growth is complex and still largely elusive. Oscillations are thought to be advantageous for biofilms to cope with nutrient starvation and chemical attacks. Recently, a minimal mathematical model has been employed to describe the oscillations in Bacillus subtilis biofilms. In this paper, we investigate four different modifications to that minimal model in order to better understand the oscillations in biofilms. Our first modification is towards making a gradient of metabolites from the center of the biofilm to the periphery. We find that it does not improve the model and is therefore, unnecessary. We then use realistic Michaelis-Menten kinetics to replace the highly simple mass-action kinetics for one of the reactions. Further, we use reversible reactions to mimic the diffusion in biofilms. As the final modification, we check the combined effect of using Michaelis-Menten kinetics and reversible reactions on the model behavior. We find that these two modifications alone or in combination improve the description of the biological scenario.

and glutamate, and two more ODEs describe the dynamics of housekeeping proteins and the enzyme glutamate dehydrogenase (GDH). In our approach, we focus solely on the two metabolites: glutamate and ammonia, and thus describe the same process in a minimal model 11 . We did this by employing a previously known minimal oscillator involving three ODEs, proposed by Wilhelm and Heinrich, aimed to describe chemical reactions 12 . The three ODEs in variables X, Y, and Z in the minimal oscillator corresponded, in our adapted model, to glutamate in the periphery (G p ), interior (G i ) and ammonia (A) respectively and were found to adequately describe the biofilm oscillations 11 . The system of equations of the basic model (BM) is as follows: The parameter names and values are listed in Table 1. Glutamate is the focal amino acid and supplied in the experimental setup. Hence, its external concentration is represented by a constant term (G E ). Its concentration within the cells in the periphery of the biofilm is dynamic and is represented by the variable G p . It is plausible to assume a positive feedback of G p on itself (k 1 G E G p ) because as more glutamate is taken up, more proteins are produced, including membrane proteins, which help take up more glutamate. It then diffuses to the interior cells (with rate k 4 G p ) and, there, it is represented by G i . We distinguish the same metabolite using two different variables based on their location in the biofilm since these regions use glutamate for different purposes. Interior cells produce ammonia using glutamate (rate k 5 G i ) since that produced by the periphery is lost to the environment and is a waste of nitrogen. Peripheral cells use the ammonia diffusing out from the interior in addition to the glutamate supplied in the microfluidics chamber to produce biomass (k 2 AG p ) 11 . However, the loss of ammonia due to diffusion to the environment (k 3 A) is much greater than its consumption to produce biomass; hence the term k 2 AG p is neglected in equation BM2.
This model describes all these processes as irreversible, which is quite unrealistic from a biological point of view. Wilhelm and coworkers 13 did present a reversible version of their oscillator. However, they considered all reactions to be reversible and assumed that the reverse reactions are considerably (10×) slower than the forward reactions, which may not be applicable for all biological scenarios. In particular, the diffusion constant is usually the same for either direction of diffusion.
In this study, we develop four sequels of the basic model to describe the oscillations in the system much more realistically, without compromising too much on minimalism. We first introduce a gradient of metabolites from the center of the biofilm to the periphery using six ODEs (Model 6ODE). Subsequently, we use realistic Michaelis-Menten kinetics instead of mass-action kinetics in addition to including some reversible reactions (Model RMM). Finally, we check the effect of using Michaelis-Menten kinetics (Model MMK) or reversible reactions (Model R) individually in two distinct versions on the model behavior.

Results and Discussion
Towards a continuous diffusion gradient within the biofilm. Having discrete interior and peripheral regions in a biofilm is quite adequate in the minimalist philosophy, but in a realistic scenario, there is no rigid boundary separating the interior from the periphery. Instead, there is a gradient of substances, the concentrations of which vary with the distance from the biofilm center. This is an important feature to test since it has not been considered by the BM 11 . Bocci et. al. have already modeled these gradients using PDEs 10 .
In our model, we restrict ourselves to ODEs since PDEs are computationally expensive. This can be done by adding several layers (based on ODEs) to the biofilm. As a first step in this regard, we consider one more layer, represented by G m , at the interface between G i and G p (see Fig. 1a,b). Furthermore, the BM only uses one variable (eq. BM2) to describe the dynamics of ammonia with the assumption that diffusion does not play a role in its dynamics because ammonia diffuses fast. We aim to test this assumption by using three separate ODEs for ammonia in the three layers of the model biofilm, namely, A i , A m, and A p . We further subdivide this model into www.nature.com/scientificreports www.nature.com/scientificreports/ two categories based on the function of the middle layer. Since the layer is added between the interior and periphery, it can either act simply as a transitional layer (referred to as s6ODE model, Fig. 1a) or it could have a more complex function (referred to as c6ODE model, Fig. 1b) as a layer that contributes to biomass formation (like the periphery) and also ammonia production (like the interior). serves not only as a physical but also a functional transitional layer, since it contributes to ammonia production (k 5 G m ) and to biomass production (k 2 A m G m ) in addition to being connected through diffusion as in s6ODE. For values of parameters, see Table 1.
These two possibilities are modeled and analyzed separately. The results are shown in Fig. 2. For the bifurcation plot of the s6ODE model, see Fig. S1.
Simulating either of these model variants leads to results that do not show qualitative differences from the three-variable version 11 , except that the oscillations are more spike-like (Fig. 2a). This is because the middle layer causes a delay in the diffusion of glutamate and ammonia. A comparison of the bifurcation plots (Figs. 2b and S1) and the analogous Fig. 2b(inset) shows that there is no drastic effect on the Hopf bifurcation point nor the amplitude. Thus, we concluded that it is not necessary to go on adding further layers to the s6ODE or c6ODE variant until a smooth gradient is created, and it is adequate to use a three-variable system instead of an n-variable system. This is in support of our assumption that an ODE model instead of a PDE model is sufficient.

Michaelis-Menten kinetics along with reversible reactions.
In order to further improve the model by using biologically more realistic kinetics, we employ Michaelis-Menten kinetics to describe the uptake of glutamate from the environment, which is mediated by transport proteins embedded in the membrane of the cell 14 . Additionally, inspired by the reversible version 13 of the BM 12 , we also make the following reactions reversible: The consumption of A and G p for protein (biomass) is considered reversible. This is biologically meaningful because biomass is subject to a permanent turnover, involving degradation to amino acids and ammonia. In the above reaction equation, we consider A on both sides for technical reasons. This formally implies that  www.nature.com/scientificreports www.nature.com/scientificreports/ the stoichiometric coefficient of ammonia is zero, which is in line with the approximative assumption that the ammonia balance is almost exclusively affected by diffusion rather than by biomass turnover. Additionally, the diffusion of glutamate between the interior and periphery and the production of ammonia from glutamate is also reversible. The diffusion of ammonia to the surroundings is a rapid irreversible process 8 . Additionally, the uptake of glutamate from the environment is also irreversible since the process is mediated by transporters on the bacterial membrane 14 . Thus, the reversible Michaelis-Menten version of BM 11 , called model RMM is governed by the following equations: We then examine the sensitivity of the bifurcation to the parameters of the model. To simplify this, we assume Note that for bifurcation analyses with respect to all parameters other than r, we use the plausible equality k −2 = k 2 , k −4 = k 4 and k −5 = k 5 . Bifurcation plots for the parameters G E , V max , K m, and r are shown in Fig. 3a-d. We used a wide scan range for each parameter because we do not have experimental values for all parameters. Once the parameters have been measured more accurately, their values can easily be included in the model.
For bifurcation parameter G E , we once again found the supercritical Hopf bifurcation just like in the BM. Interestingly, when G E was further increased, we observed another bifurcation which was laterally flipped, supercritical and conjoined to the first bifurcation, so that the diagram resembles a 'bubble' (Fig. 3a). Such a bubble-like pair of Hopf bifurcations was also observed in several models of calcium oscillations 15 . The supercritical bifurcations indicate that there are certain values of G E where the oscillations vanish with a smoothly decreasing amplitude. This implies that the oscillations will vanish with high availability of glutamate (nutrients), this was also found to be true experimentally 4 . The supercritical nature of the second bifurcation seems more biologically realistic than a switch-like transition (as would occur in a subcritical Hopf bifurcation) since the oscillatory process cannot suddenly halt. A similar effect is also observed when V max is used as a bifurcation parameter instead of G E (Fig. 3b).
A second Hopf bifurcation also occurs upon variation of the Michaelis-Menten constant, K m (Fig. 3c). When K m is very low, the saturation range is reached early, so that the kinetics becomes practically independent of G p . Therefore, no undamped oscillations can arise in that case. When K m is very high, the slope of the kinetics is low, so that the feedback is not strong enough to enable undamped oscillations.
As can be seen in Fig. 3d, oscillations occur also if r is in the range of 2 h −1 , similar to those of the forward reactions (see Table 1), which implies that oscillations can be observed when the rate constants of forward and reverse reactions are approximately equal. In our biological scenario, equal rates of forward and backward reactions are quite suitable.
Michaelis-menten kinetics -fortifying the model to make it biologically plausible. Now we make two modifications to investigate the individual effects of either Michaelis-Menten kinetics or reversibility of reactions on this bubble-like bifurcation and determine which one of them are the cause of such a bifurcation.
We modify equation BM1 describing G p as follows: The other equations (MMK2 and MMK3) remain unchanged from the BM (eqs BM2 and BM3) 11 . We found the same bubble-like pair of Hopf bifurcations that was seen for model RMM. This bifurcation could be explained as follows: If G E gets very high, the glutamate input increases considerably, so that G p increases up to the saturation range. Then, the kinetics becomes practically independent of G p , so that the essential properties of the minimal basic model are no longer granted and undamped oscillations can no longer occur (Fig. 4). While usually, the amplitude increases gradually at supercritical Hopf bifurcations, in our system, it increases quite abruptly in one of the bifurcations (Fig. 4a, right part of the bubble). A drawback of this model is the complete lack of reversible reactions which is unrealistic especially in the case of diffusion. the effect of reversible reactions. In the final modification, we investigate the effect of reversible reactions alone. The model is similar to model RMM except that it uses mass action kinetics instead of Michaelis-Menten. Thus the equation RMM1 from above becomes: And the equations RMM 2,3 are as above: www.nature.com/scientificreports www.nature.com/scientificreports/ In (a-c), we can observe a peculiar bubble-like pair of Hopf bifurcations. In the case of K m , the bubble is laterally inverted with respect to that of G E and V max indicating that increasing K m implies a steep increase in amplitude at the left bifurcation whereas a gradual decrease at the right bifurcation. For parameter values refer to Table 1. www.nature.com/scientificreports www.nature.com/scientificreports/ Similar to model RMM, r is the rate constant of reverse reactions for the bifurcation analysis with respect to these reactions and k −2 = k −4 = k −5 = r. In the bifurcation analysis for the other parameters, we assume k −2 = k 2 , k −4 = k 4 and k −5 = k 5 . This is the reversible version of BM 11 , called model R. Note that Eqs. R1 and R2 are identical to equations RMM2 and RMM3, respectively.
The system shows no change in the oscillatory behavior in comparison to the BM, in particular, bifurcation analysis with respect to G E does not show the bubble-like Hopf bifurcation (see Fig. 5a). This indicated that Michaelis-Menten kinetics is the cause for the bubble-like Hopf bifurcation, in line with the explanation for that given above. However, when the reaction for the uptake of glutamate is made reversible, the bubble-like bifurcation can be seen once again (see Fig. S2). This is a new observation and was not anticipated by Wilhelm and Heinrich in their reversible model 13 where they consider all reactions to be reversible with a rate constant of backward reaction 0.1 times that of the forward reaction. The drawback of model R, however, is that it tends to a steady state when all the backward rate constants are equal to the forward rate constants (data not shown). This can be addressed and oscillations can still be obtained if the value of k 1 is increased from 0.34 (mM* h) −1 to 0.88 (mM* h) −1 (Fig. 5b).

Quasi-steady-state approximation for the model R.
Since ammonia is a small molecule, it diffuses fast. Thus, the question arises whether a quasi-steady-state approximation 16 can be applied. In the following, we analyze the special case where reaction 3 would be very fast. We could then assume that the variable A attains a quasi-steady state: From Eq. (5), we derive a reduced system: The system shows two steady states:  www.nature.com/scientificreports www.nature.com/scientificreports/ which is the trivial steady-state (TSS) which is always positive. The Jacobian matrix for the NTSS reads:  The eigenvalues are which can be achieved with k 1 G E < k 4 . If the root value is smaller than the b value, then both eigenvalues are negative, so that the trivial steady state is a stable node. Otherwise, one eigenvalue is negative and the other one positive. The steady-state then is unstable, it is a saddle point. Furthermore, the case where b ac 4 2 < leads to complex eigenvalues and thus the steady-state is a stable focus (damped oscillations).
This implies that undamped oscillations cannot occur when ammonia diffuses so fast that it can be eliminated as a variable by the quasi-steady-state approximation. This also follows from general results by Hanusse 17,18 saying a system involving two variables and only linear and bilinear terms cannot give rise to limit-cycle oscillations. As oscillations were observed in B. subtilis biofilms, we do not consider extremely high diffusion coefficients of ammonia in this study.

Materials and methods
Simulation. For computer simulations, we used the software COPASI versions 4.16 and 4.24 19 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.
Except for k 2 , k 3, and k 4 , the rate constants were adopted from the publication of Wilhelm and Heinrich 12 and were rescaled such that the period of oscillations matched the one in the experimental work by Liu et al. 4 . The rate constant of diffusion of ammonia, k 3 , is tuned to be twice as high as that of glutamate, k 4 . This is because the diffusion coefficient for ammonia 20 is about 1.6E-05 cm 2 s −1 while that for glutamate 21,22 is about 8E-06 cm 2 s −1 . The values of K m and V max are arbitrary because we just want to analyse the qualitative effect of enzyme kinetics, without a specific enzyme in mind.

conclusion
Here, we have analysed various model variants to describe biofilm oscillations. We used ODE systems, as is often done in mathematical biology 16,23 . Bocci et. al. 10 used PDEs to describe biofilm oscillations, which is a computationally expensive approach. Our results above show that ODEs yield satisfactory results as well. ODE systems must involve at least two variables to describe limit-cycle oscillations 15,16,24 . When only linear and bilinear kinetics is used, at least three variables are needed 17,18 . Note that two-dimensional oscillator models such as the Brusselator 25 and the Somogyi-Stucki model 26 involve higher non-linearities. Therefore, our basic model is three-dimensional.
We have tested different modifications of the basic model 11 and discussed salient features and disadvantages for each. All the models analysed here involve a positive feedback, notably of peripheral glutamate on its own uptake. In many biological systems, a positive feedback is the cause of oscillations 15,16,27 .
Adding additional layers makes the model more spatial. It may be possible to have more experimental inputs to the model. For example, the diffusion rates of ammonia and glutamate could be obtained from experiments to obtain more realistic amplitudes and wavelengths of oscillation. However, this does not qualitatively change the model. The bifurcation obtained is a supercritical Hopf bifurcation just like in the previous study 11 . We conclude that the modification of adding a middle region (models s6ODE and c6ODE) does not improve the model in any way. It highlights that ODEs describing just two layers are sufficient to model this biological phenomenon and that similar observations can also be obtained, if spatial effects are ignored.
The next modification, RMM, investigates the effect of Michaelis-Menten kinetics and reversibility. It has reversible reactions and the system shows undamped oscillations when the equilibrium constants of the reversible reactions are unity, which is realistic. It also uses Michaelis-Menten kinetics which is more robust to describe the process of nutrient uptake. Thus this modification is well suited to describe the periodic halting of growth of the Bacillus subtilis biofilm and has a greater scope to include experimental parameters like K m and V max . We also obtain a peculiar bubble-like pair of Hopf bifurcations for this system. Such bifurcations were also obtained on models for cellular calcium oscillations 15 . The range beyond the second Hopf bifurcation is then called "overstimulation", meaning that at high agonist levels, no oscillations occur 26,27 .
The RMM model predicts that if the nutrient supply gets high enough, oscillations will no longer continue (see Fig. 3a). In other words, metabolic oscillations are a mechanism to mitigate the nutrient limitation and therefore, can only be observed when nutrients are scarce. This observation, although quite intuitive, was not reported in the models of biofilm oscillations shown earlier 4,10,11,28 . This prediction could be validated in experiments, and knowing this threshold could help us control and modulate biofilm growth.
A similar bifurcation was also observed when V max was used as a bifurcation parameter instead of G E . A second Hopf bifurcation also occurs upon variation of the Michaelis-Menten constant, K m (Fig. 3c). When K m is very low, the system reaches the saturation range early, so that the kinetics becomes practically independent of G p . Thus, the sole bilinear term is lost and the remainder linear system cannot show a limit cycle 16 . When K m is very high, the slope of the kinetics is low, so that the feedback is not strong enough to enable undamped oscillations. The model also retains its minimality with respect to the number of variables and the number of bilinear terms used which makes it quite easy to analyze.
In order to investigate whether the cause of such bifurcation is the inclusion of reversibility or Michaelis-Menten kinetics, we analysed the two modifications separately in the further models. We also discussed the inspiration behind these individual modifications in the respective paragraphs.
The third extension, MMK, is achieved by using Michaelis-Menten kinetics for the uptake reaction. Uptake of glutamate has been shown to be mediated by a proton/glutamate symport protein in B. subtilis 14 . Here we neglect the diffusion of glutamate from the environment to the biofilm periphery but focus mainly on the uptake of glutamate by the cells in the periphery. We found that at one of the Hopf bifurcations, the amplitude increases very steeply. This could be considered physiologically advantageous because a pronounced division of labour is already reached close to the bifurcation. Analogously, in calcium oscillations, this phenomenon is beneficial for signal recognition 29 . For the system studied here, this modification led to a revelation that oscillations vanish after a high enough value of the bifurcation parameter. In terms of G E , it suggests that beyond a threshold value, the oscillations will no longer persist. This seems quite a likely case since oscillations arise in order to allow the interior of the biofilm obtain a steady input of glutamate. When glutamate is supplied in large excess, the interior will always have a steady supply of it and thus the oscillations will no longer be necessary. From a mathematical point of view, we observe that the saturation kinetics viz. Michaelis-Menten kinetics are crucial for obtaining such a bubble-like Hopf bifurcation. The major contribution of this model is that it identifies the cause of the bubble-like Hopf bifurcations.
The model R has already been inspired by a model established earlier 13 . We adopted and improved it by applying reversibility to only those reactions which are really reversible. For instance, the loss of ammonia is a highly irreversible reaction since ammonia is a small molecule that diffuses very rapidly once it is generated. From a physical point of view, all diffusion processes are reversible. However, when the concentration differences are large, the reverse process can be neglected. This was also done in models of calcium oscillations with respect to the calcium concentrations in the cytosol (low) and endoplasmic reticulum and external space (high) 26,27 . As the difference in glutamate concentrations between the peripheral and interior regions is less pronounced, we here extended the BM so that glutamate diffusion is reversible.
In model R, we found the same bubble-like bifurcation that was seen for the Michaelis-Menten model only when the reaction for the uptake of glutamate (given by the term k 1 G E G p ) is made reversible (see Fig. S2). This was not anticipated by Wilhelm and coworkers 13 , who also proposed a fully reversible version of their model, and is thus an interesting new observation. Biologically it implies that when there are reversible reactions involved, there will be a stabilizing effect in the system that will prevent the amplitude to increase indefinitely with the bifurcation parameter. However, the uptake of glutamate is an active process and, therefore, can be modeled to be irreversible. In that case, only a single Hopf bifurcation is observed (Fig. 5a). The drawback is that the system tends to a steady-state (i.e. oscillations vanish) when all the backward rate constants are equal to the forward rate constants, which means that the equilibrium constants of the reversible reactions are unity. This can be resolved by increasing k 1 indicating that a higher glutamate uptake rate is required for oscillations to persist at greater reversibility. Thus, this modification helps to correct the glutamate uptake rate k 1 . It also helps identify the reaction that causes the bubble-like Hopf bifurcation. www.nature.com/scientificreports www.nature.com/scientificreports/ This work highlights the importance of simple ODE based models to describe the observations of Liu et al. 4 . Each modification presented here reveals a new piece of the puzzle which, when put together, will help us see the broader picture that is biofilm oscillations. The extended models presented here are still simpler and easier to handle than the Liu model (which involves 6 ODEs) and describe the experimental observations equally well. The aim of both the Liu model and ours is to describe biofilm oscillations qualitatively or at best semi-quantitatively. It is difficult to say whether the extended models describe the observations better than the basic model because no comprehensive parameter scan had been performed in experiment. As the extended models involve more parameters, it is likely that they can be fitted better to data obtained in the future.
The transition between oscillatory and stationary regimes can be gradual as in the case of model RMM or abrupt as in the case of model MMK which makes the BM flexible to describe two contrasting scenarios. An additional experimental investigation will provide more clarity about the exact nature of oscillations in the biofilm of Bacillus subtilis.
Ethical approval and informed consent. Ethics approval was not required for this study and not applicable

Data availability
No data are associated with this article.