Perspectives in mathematical modelling for microbial ecology

Although mathematical modelling has reached a degree of maturity in the last decades, microbial ecology is still developing, albeit at a rapid pace thanks to new insights provided by modern molecular tools. However, whilst microbiologists have long enjoyed the perspectives that particular mathematical frameworks can provide, there remains a reluctance to fully embrace the potential of models, which appear too complex, esoteric or distant from the “real-world”. Nevertheless there is a strong case for pursuing the development of mathematical models to describe microbial behaviour and interactions, dynamically, spatially and across scales. Here we put forward perspectives on the current state of mathematical modelling in microbial ecology, looking back at the developments that have defined the synergies between the disciplines, and outline some of the existing challenges that motivate us to provide practical models in the hope that greater engagement with empiricists and practitioners in the microbiological domain may be achieved. We also indicate recent advances in modelling that have had impact in both the fundamental understanding of microbial ecology and its practical application in engineered biological systems. In this way, it is anticipated that interest can be garnered from across the microbiological spectrum resulting in a broader uptake of mathematical concepts in lecture theatres, laboratories and industrial systems.


Introduction
Mathematical modelling and its applications can be traced back to ancient civilisations as an attempt to understand and analyse the world around them. Microbial ecology, on the other hand, is a much more recent branch of science, having its roots in the medical setting, before Sergei Winogradsky took the theory into the environment with the discovery of chemosynthesis at the turn of the twentieth century. However modern microbial ecology, which considered not only microbe physiology but also their interactions with the environmental and other organisms, only took form about half a century ago through the study of anaerobic microbes in the rumen of cattle (Hungate, 1960).
The ecologist and mathematician Pierre-Franç ois Verhulst developed the logistic equation to describe the growth of populations, which included terms for intrinsic growth rate and a carrying capacity (Verhulst, 1838). In 1910, Alfred Lotka proposed a pair of equations, heavily influenced by the logistic equation, to describe autocatalytic chemical reactions (Lotka, 1910), which developed into the well-known predator-prey Lotka-Volterra model describing the dynamics of biological systems, forming the basis for ecological modelling. Nearly 30 years later, Jacques Monod looked to test the logistic equation, which had recently been rediscovered by Pearl and Reed, who used it to predict population growth in the United States (Pearl and Reed, 1920), and selected bacteria for the purpose, due to the doubling-time of these organisms being far shorter than the larger lifeforms originally the subject of the model (Monod, 1942).
Monod's equation for microbial growth (Eq. (1)) actually follows the same form as Michaelis-Menten's equation describing enzyme kinetics (Eq. (2)) (S) = max S K S + S (1) where the specific microbial growth rate = (dX/dt)/X (X is the concentration of biomass growing on substrate, S) is related to the rate of product formation dP/dt, the maximum specific growth rate max is equivalent to the enzyme concentration, E, multiplied by the turnover number, k 2 . The enzyme-substrate (ES) binding rate constant is denoted by k 1 , and K S is an inverse measure of the microbial substrate affinity, denoted as K m in Michaelis-Menten; the inverse measure of substrate affinity with the enzyme (not shown). Monod's model for growth is not applicable to all systems as it is based on assumptions (e.g., fixed cell composition, saturation kinetics) that are not universally relevant to all environments (or even phases of growth), and others have subsequently built upon the work of Monod in proposing new forms or expansions to the model (Moser, 1958;Contois, 1959;Andrews, 1968).
Whilst it is acknowledged that mechanistic approaches in microbial ecological modelling have helped make sense of empirical observations, or have tried to express some physical first principle rules governing behaviour, one must still consider the underlying objective of any modelling exercise, the outputs and precision required, and the effort needed to parameterise and simulate the model under the desired conditions. As the volume, resolution, quality and types of data continue to grow and expand, modelling must move apace and find ways to best incorporate the new knowledge accessible from this information. Indeed, stochastic models, models that incorporate metabolic fluxes and heuristic approaches, such as genetic algorithms, or even equation-free models (DeAngelis and Yurek, 2015;Ye et al., 2015) are becoming popular in light of this increase in data and information. Nevertheless, methods such as the latter are typically "black-box" and fail to allow for a deeper understanding of the underlying mechanics of a process, essentially restricting the advancement of knowledge in favour of prediction. Although this manuscript focuses primarily on deterministic modelling, we do not discount the contribution from the wider field of mathematical modelling, but stress the importance of understanding the mechanics of the problem to be addressed prior to model selection and application, as well as the real need to account for complexity given the available information, data and theorems we currently have.
A recent review of current thought in microbial ecology have described the importance of considering ecological theory that underpins much that has emerged through experimental endeavour, and how some of the challenges in the field are now being addressed through theory driven by tools and concepts created thanks to advances in molecular and genomic techniques (Prosser et al., 2007). Whilst we describe some of these aspects, this review presents a summary of significant developments related specifically to the use of mathematical modelling for microbial ecology, together with current perspectives and future challenges aimed at stimulating both researchers and practitioners working in the many related fields.

Current developments
We present here an overview of developments in mathematical modelling with examples related specifically to engineered biological systems, where microbial ecology and its control are key. However, we do acknowledge that other areas such as predictive modelling in food microbiology are also worthy of focus, but have decided on this field for reasons of brevity and focus. The emergence of novel insights in ecological theory insights guided by greater molecular scrutiny has meant that engineers and microbiologists alike have access to data and knowledge that shift microbial ecology from purely a passive discipline towards one that allows active engagement across scientific fields. Advances in empirical study of microbial systems now allow scientists to control, manipulate and predict theoretical concepts at the microbial scale with much more certainty and scope (Prosser et al., 2007).
Whilst practical experience and empiricism have led to development through observation, trial-and-error or by chance, many of the greatest advancements in engineering have come through the acquisition of a theoretical understanding guided by mathematics. Advances in microbial ecology methodologies over the last few decades have yielded the tools for qualifying and quantifying microbially driven systems, allowing theorists to better test hypotheses resulting in the potential for better engineering design and operation of biological processes. Furthermore, there has been a historical transition from microbiology of undifferentiated biomass or single strains, towards the concept of multi-species microbial ecosystems incorporating an ecological dimension, which has opened up the possibilities for both individual and community based mathematical modelling. Table 1 Examples of accepted models for different engineered biological processes.

Models Dimension Notes
Primary clarification Lessard and Beck (1988) 1 Simple 5-layer model Activated sludge Henze et al. (1987) 13 (9) Model of 8 biochemical processes Henze et al. (1999) 19 Model extensions include Bio-P Secondary clarification Takács et al. (1991) 1 10-layer model Anaerobic digestion Batstone et al. (2002) 32 Model of 19 biochemical processes Biofilms Picioreanu et al. (1998) 1 Cellular automata Noguera and Picioreanu (2004) 2 2D automata &ODE approaches Plant-wide Jeppsson et al. (2006) 114 Benchmark models Although many conventional engineered biological systems have a well-established history, the development and acceptance of mathematical models describing biological phenomena in these systems is much more recent. The activated sludge process for wastewater treatment, a century old in 2014, did not attract the interest of modellers until the late 1970s.
With the rapid advances in computing capabilities, the scope and extent of mathematical modelling broadened to allow for previously intractable systems to be analysed, simulated and developed with greater accuracy and speed. This burgeoning in silico capacity provided the support for the development of the first reference model for engineered biological systems, the Activated Sludge Model No. 1 (ASM1) (Henze et al., 1987), which built on the model structure first proposed by Dold et al. (1980). ASM1 and its extensions have allowed both theorists and practitioners to gain a more formal understanding of the biochemical mechanisms involved in wastewater treatment, such as carbon oxidation, nitrification and denitrification, through analytical studies (Nelson and Sidhu, 2009), computer simulation (Gernaey et al., 2000) and process prediction (Hu et al., 2003). These mechanistic models allow for the coupling of mass-balance derived equations describing biological growth combined with physical and chemical properties of the system. The relationships can then be utilised in process characterisation, design and control, but are sensitive to the assumptions built into the model.
As interest in the renewable energy sector increased at the turn of the twenty-first century due to a growing awareness of the limitations and impacts of energy intensive aerobic processes, anaerobic digestion became a fertile area for both mathematical modelling and microbial ecology, much in the same way that rising fossil fuel prices in the 1970s stimulated a first wave of anaerobic digestion research, primarily through the efforts of J.F. Andrews and his work on dynamic modelling, stability and control (Andrews, 1969(Andrews, , 1971Graef and Andrews, 1974).
Subsequent developments in anaerobic digestion modelling, and identification and characterisation of functional groups culminated in the publication of an advanced model describing the complex stoichiometry and kinetics of a standardised anaerobic process, the Anaerobic Digestion Model No. 1 (ADM1) (Batstone et al., 2002).
However, it has become apparent that ADM1 is limited and it has been postulated that with greater empirical understanding (e.g., more realistic parameterisation) and process extensions (e.g., inclusion of phosphorus related components), a better trade-off between model realism and complexity may be achieved (Jeppsson et al., 2013). Additionally, as insights into the fundamental mechanisms of microbial ecology within anaerobic digesters emerge, it is vital that this knowledge is integrated into models. However, finding the best ways to develop a workable 'virtual' AD process that incorporates this information remains an open challenge.
Aside from the suspended growth systems described, models for fixed or attached growth have also received a lot of attention, with applications ranging from healthcare to water treatment systems (Böl et al., 2013). A specific extension of the chemostat model used for suspended growth was analysed for stability by considering separate rate models for planktonic growth and floc formation (attachment rate) and separation (detachment rate) (Fekih-Salem et al., 2013). The specific form of the model considers that the growth rate and dilution terms are dependent on the density of biomass in the system (cf. Section 4.3). Biofilms or 'wall growth' models have received a large amount of attention from microbiologists and modellers alike. The combination of succession and spatialisation makes these systems very interesting to study, with a host of methods and concepts applied to understanding their characteristics, predicting their behaviour and managing their formation (Picioreanu et al., 1998;Noguera and Picioreanu, 2004;D'Acunto et al., 2015). Typically these spatial models are more complex than their compartmental Ordinary Differential Equation (ODE)-based counterparts, relying on Partial Differential Equations (PDEs) to describe the transport terms. Nevertheless, biofilm modelling is fairly mature and presents a necessary means to bridge the gap between ecological theory and engineering practice.
A summary of some of the most popular models for engineered biological systems is shown in Table 1.
The ability to effectively and robustly model engineered biological systems is fundamental in their progression from being poorly understood black box processes to facilities that can be operated, controlled and optimised despite significant levels of uncertainty in their inputs. A comprehensive understanding of the role of microbial communities, their characterisation and quantification, will lead to better system design and operation so long as the models can integrate the ever more abundant and resolved information generated by rapidly evolving microbial analysis techniques and technologies. The future of modelling of engineered biological systems can only be strengthened by the evolving techniques and tools for observing and characterising the form and function of their microbial ecology, complimented by the involvement of mathematicians in order to handle their inherently complex nature.
It should be noted here that over recent years there has been a step-change in the ability to observe, characterise and analyse biological systems brought about primarily by the developments in molecular tools such as next generation sequencing and quantification techniques. This has resulted in vast amounts of biological data becoming available at relatively low cost and with fast acquisition times. Although not the focus of this review, we feel it is important to highlight the impact that these tools are having in modelling of microbial ecology, and a summary is presented in supplemental material.

Using mechanistic models to understand biological phenomena
When microbial ecosystems are not at equilibrium, or evolve to different equilibria depending on the operational conditions and initial composition, the study of mechanistic models is often well suited to propose explanations and produce predictions of their evolution through time. An example of a fairly well-studied model is given by the classical chemostat model shown in Table 2A, where S and X are substrate and biomass concentrations, respectively, D is dilution rate, Y is the yield coefficient, and is the growth function for X on S. In this section, we give two illustrations of situations for which the dynamical outcome is not straightforward and can even be surprising.
Consideration of aggregation in chemostat mode operations. Biomass aggregation is a common phenomenon observed in continuous cultures, which results in the formation of flocs or biofilms. Modelling in a simple way the difference between "attached" and "free" biomass as two compartments (whose concentrations we denote respectively by X a and X f ) leads to a system of three differential equations, as shown in Table 2B.
The system contains parameters g a (·) and g d (·) that model the attachment-detachment process. This general formulation gathers various models already proposed in the literature, such as adaptive nutrient uptake (Tang et al., 1997), wall attachment (Pilyugin and Waltman, 1999), intestine model (Freter et al., 1983) or flocs (Haegeman and Rapaport, 2008), amongst others, that correspond to different choices of functions g a (·) and g d (·).
The main point in the analysis of this model is to assume that the attachment and detachment processes are intrinsically fast phenomena compared to bacterial growth (Thomas et al., 1999). Then, the model (Table 2B) can be reduced under the quasi-state approximation to a system of two equations only, with the same structure as the standard chemostat model (Table 2A), describing with good accuracy the dynamic of the total biomass X = X a + X f , and the substrate with new uptake and removal rates depending on S and X (Table 2F) expressed by: where p(S, X) is a certain function derived from the mathematical analysis. It can be recognised from the equations in Table 2F the classical chemostat model with a justification of density dependence for the growth function (S, X) (as considered by Lobry et al., 2005), and more surprisingly, the equivalent removal rate, D 1 (S, X), that depends on S and X.
The mathematical analysis of such a model has revealed some non-intuitive properties, for the single species case and its extension to multi-species cases: • multiplicity of positive equilibria even when the growth functions a (·) and f (·)are monotonic, with the possibility of the wash-out equilibrium to be attractive (Fekih-Salem et al., 2013), • possible coexistence of different species when one of them aggregates (Haegeman and Rapaport, 2008;Fekih-Salem et al., 2013, • appearance of stable periodic solutions for multiple species models (Fekih-Salem et al., 2015).
It is important to note that these phenomena cannot occur in the model without aggregation.
Consideration of fluctuating environments. In the classical chemostat model with multiple species (cf. Table 2C) it is well known that one cannot have more than one species present at steadystate (Hsu et al., 1977) (a property that is usually called Competitive Exclusion Principle (CEP)). Since this (mathematical) result contrasts with reality, there has been a concerted effort to understand model adjustments that can result in coexistence. Several possible mechanisms have been proposed to explain coexistence, such as intra-specific interactions or "crowding effects" (Lobry et al., 2005). One such modification could be the consideration of a variable dilution rate t → D(t) replacing the typical constant value. It has been recognised for over 30 years that a periodic time varying dilution rate could lead to the coexistence of several species . More recently, sophisticated time-dilution rate relationships have been considered that can model more realistically the seasonality of an environment: piecewise constant, "almost periodic" dilution rate and Piecewise Deterministic Markov Processes Lobry, 2013;Benaï m and Lobry, 2014).
Nevertheless, one has to be aware that during the process some species can reach such small concentrations that the validity of the deterministic population model is questionable (see further the atto-fox problem).

Understanding extinction through stochastic models
As is well understood the modelling of population dynamics with deterministic differential equations is an approximation of more complex probabilistic birth and death processes. One of the main drawbacks of modelling with differential equations is the fact that we represent the size of populations by continuous variables, which prevents them becoming exactly equal to zero. To understand this consider the trivial differential equation which represents the evolution of the size of a population when the growth rate depends on time. With initial condition x(0) = 1, it can be easily checked that x(t) decreases until t = 10 and then increases for t > 10 and is equal to 1 for t = 20. Thus, if x(t) represents a population, it is apparently persistent. But, assume that x = 1 represents 10 4 individuals; since, as x(10) = e −50 ≈ 1.93 × 10 −22 we see that x(10) represents approximately 1.93 × 10 −18 individuals, which is absurd. It was noticed by Mollison (1991) that in the work by Murray et al. (1986), who developed a model of the propagation of rabies in Great Britain, a variable representing the number of foxes per square kilometre could decrease (in the model!) down to 10 −18 , and the expression "atto-fox" problem was subsequently coined for this phenomenon. Campillo and Lobry showed that the very basic Rosenzweig-MacArthur model presents an atto-fox problem for quite reasonable values of the parameters even if one unit of resource represents 10 6 individuals (Campillo and Lobry, 2012). But, conversely, a drawback of stochastic models is that, due to unlikely events of negligible probability (say <10 −9 ), most models predict extinction with probability one, which is not qualitatively accurate in describing reality since, if it is "mathematically true" that extinction will occur, it will occur in a time that can be greater than the age of the universe. More precisely, consider the logistic model Here a represents the exponential growth rate of the population (when the population is small, there is no competition for the resources) and K represents a limiting term due to the resources. A natural stochastic process associated to this is the birth and death process with density-dependence parameters (Kurtz, 1970). That is, we have three parameters, b, d, c. Each individual gives birth at constant rate b and dies at rate d. Moreover, for each couple of individuals, one bacterium of this couple dies at rate c (a death by competition). If b and d are much larger than c (b, d c) (cf. Theorem 17 in Méléard and Villemonais, 2012), then the random dynamics is close to the deterministic model with parameter a = (b − d) and Although the long time behaviour of the approximating differential system is well known: the population goes to a positive limit, it is not the same for the random system, where the population tends to extinction with probability one. This can seem inconsistent, but it is not. Indeed the extinction time, T 0 may be very large and the population size can remain around a form of equilibrium over a long period of time. Thus, although the population goes to extinction, we do not see it in simulation and during experiments. The equilibrium before extinction is called, in probability theory, a quasi-stationary distribution (or the Yaglom limit) (Méléard and Villemonais, 2012). An extensive bibliography on this subject can be found in a manuscript by Pollett (2015). If a random process attains a quasi-stationary equilibrium then the extinction time T 0 has the property to be exponentially distributed and, as a consequence, there is no memory, and the expectancy lifetime then remains constant through time (cf. Méléard and Villemonais, 2012, Section 2.1.5). This phenomenon was observed, for instance, with experiments in Figure 2 in Aalen and Gjessing (2001) and Figure 1 from Carey et al. (1992).
Following very recent work (Chazottes et al., 2014), one can see that the logistic birth and death process attains its quasi-stationary equilibrium at a time around K ln(K) unit of time, rests around its equilibrium during an exponentially long time interval in K unit of time, and then becomes extinct. The quasi-stationary equilibrium is found close to the deterministic limit.
Let us now return to the chemostat. Adding randomness to the chemostat model probably dates from studies carried out in 1979 (Crump and O'Young, 1979;Stephanopoulos et al., 1979). In contrast with Stephanopoulos et al. (1979), where the stochasticity was introduced according to an ad-hoc approach, a more recent study (Campillo et al., 2010) proposes a family of models where the randomness emerges from the microscopic dynamics. Various stochastic models are developed such as SDEs (stochastic differential equations), pure jump process and hybrid process (processes whose components are deterministic with a jump component). In particular, they recover the model already proposed by Crump and O'Young (1979) and give an overview on the literature on the subject. More recently, other studies have developed an individual-based approach to justify some model with partial differential equations and the standard chemostat equation in large populations (Table 2A) (Campillo and Fritsch, 2014;Fritsch et al., 2015). Other individual based models were introduced in order to understand the genetic evolution of the population (Champagnat et al., 2014). For these kind of models, the challenge is to obtain, as for the logistic model (Chazottes et al., 2014), the different time-scales of the process. Namely, a time-scale where the process behaviour attains its quasi-equilibrium (that is satisfactorily close to the deterministic equilibirum) and an estimation of the extinction time. Preliminary results in this area are the subject of current research.
For instance, consider the Crump-Young model (Crump and O'Young, 1979), which is viewed as a benchmark example with results also applicable to other chemostat models. The model is described briefly as follows. The nutrient concentration evolves continuously but is dependent on the population size, which is a birth and death process (pure jump process) with coefficients depending on time through the nutrient concentration. The random fluctuations in this model are only due to the individual births and deaths of bacteria; see for instance Fig. 1, as one example. More precisely, we keep the notation for S, X, D, S in and we denote the birth rate per individual by b and the death rate per individual by d. These two new parameters depend on the nutrient rate S (and could depend on X in a ratio-dependence type model) and is close to the main deterministic model when b = /Y and d = D. Despite its simplicity, the above results on the logistic model do not apply to the chemostat model. That is, even if this model is simpler than the individual-based model, one can not estimate the time of extinction or the quasi-equilibrium as in (Chazottes et al., 2014).
However, an algorithm amenable for simulation is described by Campillo et al. (2010). A trajectory of this process is given in Fig. 1; namely a possible realisation of this process.
Nevertheless, even if the quasi-stationary distributions and the extinction rates are not well understood in the chemostat model, they are amenable to simulation. Indeed, there exists some numerical method to approximate these quantities. One can use the particle method described in Section 6 of work by Méléard and Villemonais (2012, Section 6). This method was introduced by Burdzy et al. (2000) and Moral and Miclo (2000) and has generated an enormous amount of research (Asselah et al., 2011;Cloez and Thai, 2013;Groisman and Jonckheere, 2013;Villemonais, 2011). Another method was recently studied and is based on stochastic approximation, as described in Benaïm and Cloez (2014) and Blanchet et al. (2014).

Generalisation and complexity
In this section, we use the term "complex" in a very general sense. It is used either by its classical meaning, that is the capacity of the model to reproduce properties that emerge from interactions within its basic entities, but also to designate complicated models we have described previously that are not necessarily complex in the normal sense. In previous sections, models of microbial ecosystems aim at integrating "all" the knowledge available at a given time in order to develop virtual biosystems. However, whilst this approach is well suited for simulation purposes (notably because they limit the number of costly and timely experiments), it does not necessarily allow for a better understanding of fundamental principles in biological processes. An alternative way to proceed is to develop simpler models from basic concepts, by reducing such complicated models, or directly from available data to take advantage of the knowledge a rigorous mathematical analysis may provide to practitioners. The question of generalisation refers to the ability to extract generic knowledge about complicated/complex models by studying simpler models. As an example, we show how the study of two-step models, and their increasing complexity over time, has led to a better understanding of a number of specific microbial interactions as commensalistic, syntrophic or mutualistic relationships from which specific properties may eventually emerge.
Two-step mass-balance models are often used in biotechnology. Under aerobic as well as anaerobic conditions, a number of fundamental biological processes have been shown to be adequately described by such models. Also known as commensalistic systems (i.e., a system where one species grows on the product of another one), they present the advantage of being complex enough to capture important process properties while being simple enough to be mathematically studied. In particular, the number of steady-states and their stability as a function of model inputs and parameters may be investigated. Two-step models are commonly used to describe commensalistic microbial systems that take the form of a cascade of two biological reactions where one substrate S 1 is consumed by one microorganism X 1 to produce a product S 2 , which serves as the main limiting substrate for a second microorganism X 2 , as shown in Table 2D. A more general two-steps model is indicated by the equations in Table 2E.
The different analyses of this class of models available in the literature essentially differ (i) in the way the growth rate functions are characterised and (ii) whether a specific input for S 2 is considered or not (i.e., the presence of a term S in 2 in the dynamic equation of S 2 ). Table 3 presents the main studies of these models (cf. notation in Table 2E) and summarises their contributions for practitioners.
In all cases, the mathematical analysis of the proposed models have allowed authors to identify hypotheses about growth rate functions that were likely to explain their experiments and to identify conditions under which particular communities are stable or not (Weedermann et al., 2015;Wade et al., 2016;Sari and Harmand, 2014).
The studies reported in Table 3 are examples of models using resource-based growth functions. However, these may have limited application or fail to capture the true dynamics of the system under investigation and, as such, density-dependent growth functions have been proposed as an alternative to these classical models (See Section 4.3).

Variations in chemostat theory
Several extensions of the basic chemostat model (Table 2A) dealing with multi-species, spatial inhomogeneity, and attached biomass, etc., have been addressed in the literature in order to have more realistic representation of real bioprocesses. The focus here is on recent results that deal with the effects of those considerations on the conversion yielding of the chemostat and its stability.
Influence of spatial configurations. For the modelling of tubular reactors or situations for which the diffusion characteristics are uniform in some direction, the mathematical analysis and numerical simulations that are available for partial differential equations, of reactor-diffusion type, are quite efficient. For less homogeneous media, such as soils or lakes with complex dynamics, finding satisfactory and tractable representations is still challenging. A simple alternative to PDEs is to consider a network of interconnected compartments.
As an example, for a given bacterial strain, three different motifs are considered, as shown in Fig. 2. These motifs have the same total volume and flow rate, and can theoretically show that there exists a threshold on the input concentration such that the following property is satisfied at steady-state (Haidar et al., 2011): above the threshold the most efficient configurations for converting the input substrate are serial, while below the best ones are parallel.
Bio-augmentation. The Competitive Exclusion Principle (CEP) (cf. Section 3.1) has been extended when the growth functions are not necessarily increasing, such as with Haldane's law (Butler and Wolkowicz, 1985). For those functions, the CEP states that (generically), at most, one species survives at steady-state, but the winner may depend on the initial condition. This feature can be exploited as a way to stabilise an unstable chemostat with bio-augmentation: consider a chemostat with a single species that presents a nonmonotonic growth rate, and a dilution rate allowing the existence of a survival attractive equilibrium (see for instance a Haldane growth function, (S) = ¯ S/(K S + S + S 2 /K i ), where the coefficient K i is related to the growth inhibition for large concentration of substrate, plotted in red in Fig. 3). In such a case, the wash-out equilibrium can be attractive leading to the (undesired) property of bi-stability of the dynamics. A common way to stabilise globally the dynamics about the (desired) positive equilibrium is to act on the dilution rate. An alternative is to add another species, with a monotonic growth rate that performs worse at steady-state but is more robust in the transient phase: Fig. 3 show two possible candidates depicted in green and blue.
Recently, it has been shown that another way to stabilise globally the chemostat model with non-monotonic growth function is to consider particular spatial configurations with a buffer (Fig. 4) while it is not possible with serial or parallel configurations (Rapaport et al., 2014). The buffer chemostat acts as a protection zone during the initialisation of the bioreactor.  (Benyahia et al., 2012) 1 (S1, S2) 2 (S1, S2) 0 0 1 Mutualistic relationships with explicit expressions of growth rates (Kreikenbohm and Bohl, 1986) 1 (S1, S2) 2 (S1, S2) 0 0 1 Study of model in Kreikenbohm and Bohl (1986) (Burchard, 1994) 1 (S1, S2) 2 (S1, S2) 0 0 1 General study of model in Kreikenbohm andBohl (1986) (El Hajji et al., 2009) 1 (S1, S2) 2 (S1, S2) 0 0 1 Analysis of an extended class of models considering an input for S2  1 (S1, S2) 2 (S1, S2) Study of the influence of maintenance on process stability (Xu et al., 2011) 1 (S1, S2) 2 (S1, S2) Generalisation of results from Xu et al. (2011) (Sari and Over-yielding with fluctuations. The chemostat device has been originally designed to be operated at steady-state, but the effect of a non-constant periodic input of resource can be studied. If we compare biomass growth in one period with the one observed with the averaged inputs during the same period then, due to nonlinearities, the results are different. It can be shown in the mathematical model of the chemostat that when the growth rate is a concave density dependent function over-yielding is not possible while, when the growth rate is a ratiodependent function such as the Contois law (one the simplest laws that is ratio dependant: (S) = max S/(K x · X + S)), over-yielding does occur (Caraballo et al., 2015).

Density/ratio dependence in growth functions
Understanding the rate at which microorganisms grow is clearly one of the main challenges in microbiology over the last century. However, as opposed to physical phenomena, the laws of biology are unknown although there exist a number of attempts to give a physical basis to rate expressions (cf. Section 4.4) most models rely on completely heuristic functions, such as Monod (monotonic), Haldane (non-monotonic) or Contois functions. In their review, Bastin and Dochain have identified more than 60 different expressions used for modelling growth rates in biotechnology (Bastin and Dochain, 1990). However, it should be noted that the use of one function instead of another may completely change the qualitative behaviour of a mathematical model. For instance, using a Haldane function instead of a Monod function in the well-known chemostat model introduces bistability in the system as long as D < (S in ). Validating, or more often invalidating, biological models on the basis of qualitative considerations as proposed in the previous section appear to be of much higher value (because of this qualitative character) than considering criteria based on the fact a model is or not able to match data.
Studying whether the density of biomass influences its own growth has been debated for decades and actually remains questionable. In particular, it was noted that Monod's prediction for growth was confirmed for pure cultures growing on glucose (Grady et al., 1972), but the results consistently diverged from this prediction when working with mixed cultures (e.g., in wastewater treatment or fermentation processes) (Grady et al., 1972;Daigger and Grady, 1977). However, efforts to identify more complicated models (predator-prey dependent growth rates vs. prey-dependent ones) have meant that resource-dependent models have emerged and the use of density-dependent ones has declined (Jost, 2000). Regarding processes in play, it was hypothesised that the higher the prey density, the lower the growth rate, as predicted by densitydependence (Lobry and Harmand, 2006). Following an idea by Arditi to validate qualitatively such a hypothesis, Harmand and Godon have reviewed the literature and shown that membrane bioreactors (MBRs) in series could be used to investigate this question (Harmand and Godon, 2007). In a recent book presenting the state-of-the-art of their work on the subject, developed over the last 40 years, Arditi and Ginzburg argue that density-dependence, and ratio-dependence in particular (e.g., Contois' function), would be more general than is actually recognised currently in the ecological field (Arditi and Ginzburg, 2012). It should be stressed that such a switch in modelling growth rate functions is more than just a simple story: modelling the competition of microorganisms on a single substrate with a density-dependent growth function allows for species coexistence, which is observed in practice. Doing so would fundamentally change the intrinsic properties of bioprocesses and yield new control and optimisation strategies. Previous theoretical and experimental work with macro-organisms (crustaceans Saïah, 1992 andwolves Jost et al., 2005) suggest that spatial heterogeneity is the fundamental mechanism of ratio dependence. This heterogeneity can be imposed externally (Arditi and Saïah, 1992;Poggiale et al., 1998), e.g., with a physical refuge in which the resource is produced but into which the consumer cannot enter. The heterogeneity can also arise because of internal causes in the system, as for example with movement rules of the consumer that accelerates in the direction of the resource gradient (Arditi et al., 2001;Tyutyunov et al., 2008). Further work with microorganisms is likely to bring additional explanations, such as resource exhaustion in the immediate vicinity of individual bacteria (Lobry and Harmand, 2006).

Beyond thermodynamics: linking energy balances with biological rates
The modelling of microbial growth relies on empirical laws. Microbes, as the simplest living things on earth, can also be seen as structures representing a transition between physical and biological systems. The study of their growth in physical terms, therefore, constitutes an ideal thinking ground for crossing the disciplinary boundaries between biology and physics. Microbial thermodynamics has especially been a matter of study since the 1960s (McCarty, 1965). Energy balances of microbial growth have been investigated in detail by several authors (McCarty, 1965;Roels, 1980;Rittmann and McCarty, 2001;Kleerebezem and van Loosdrecht, 2010). Microbial metabolism is described as a combination of a catabolic energy yielding reaction and an anabolic reaction leading to the synthesis of new cell material. The overall metabolic reaction was defined as a linear combination of the catabolic and the anabolic reactions and is an irreversible process leading to the dissipation of energy. Significant progress came from studying growth yields and energy dissipation in greater detail. It indeed appeared that the energy dissipated per unit of newly formed biomass was strongly constrained and was dependent primarily on the carbon source used for growth. For heterotrophs, the amount of energy dissipated per unit of newly formed biomass varied within a relatively narrow range that depended mostly on thermodynamic properties of the electron donor molecule, mainly its oxidation degree and the length of the carbon chain (Heijnen and Vandijken, 1992). By taking into account these constraints on the amount of dissipated energy per unit of biomass formed, it was possible to link energy and matter balances during microbial growth and to predict growth stoichiometry (for a review see Kleerebezem and van Loosdrecht, 2010). However, in order to go further and to model microbial community dynamics using thermodynamic principles, the link between these energy balances and growth rates was still missing. Desmond-Le Quemener and Bouchez recently proposed a thermodynamic theory of microbial growth that allows for the link between thermodynamic balances and growth rates to be made explicit (Desmond-Le Quemener and Bouchez, 2014). In their work, they followed the initial intuition of Alfred Lotka that suggested that the similarity between individuals was an invitation to imagine a "statistical mechanics of living beings" (Lotka, 1922). They showed how systems composed of microbes in contact with molecules could be likened to ensembles described by the laws of statistical physics. Based on thermodynamic balances established by previous authors (see review by Kleerebezem and van Loosdrecht, 2010), an "activation energy" of a microbe was defined and, thus, the probability for an elementary division act to be triggered could be determined. A growth equation could then be proposed, which links a flux (the growth of microbes) to a force (the energy density): where stands for the microbial growth rate, [S] is the substrate concentration, and E M , E dis , E cat are, anabolic, dissipated and catabolic exergies associated to an elementary microbial division act, respectively. These values can be calculated using previously developed methods (Kleerebezem and van Loosdrecht, 2010). In addition to these variables, the equation comprises two parameters: max , which represents the frequency at which an activated microbe is able to divide and V harv , which represents the volume that a microbe can explore in order to harvest substrates during an elementary division act.
The proposed equation allows to adequately model experimental data of microbial growth. More importantly, it also allows original predictions in relation to the microbial isotopic fractionation phenomenon to be made and its dependence on energy variations induced by the use of different isotopic isomers (isotopomers). For a given atom, the composition of the nucleus in terms of proton and neutron indeed influences the stability of chemical bonds with other atoms. It, therefore, slightly changes the "Zero Point Energy" between molecules composed of different isotopes (Zeebe and Wolf-Gladrow, 2001). According to the flux force relationship derived by Desmond-Le Quemener and Bouchez (2014), these differences in energy should induce differences in rates, and the isotopic fractionation phenomenon could thus be viewed a kinetic consequence of the differences in energy contents of isotopomers. Theoretical predictions have been questioned using experimental data and first evaluations were found to support the theory. A more in depth evaluation is however needed in order to precisely evaluate the actual validity of this theoretical formulation and to eventually circumvent the points requiring additional work and clarification. More generally, this trial to relate available energy gradients to biological rates can be seen as an incentive for microbial ecologists to tackle the challenge of thinking about the fundamental principles underlying the phenomena they study.
Such scientific approaches are of course risky and often require crossing disciplinary boundaries between biology, mathematics and physics. However, the development of such initiatives is today urgently required to elaborate more predictive models and to progressively set up the principles of a sound ecological engineering of microbial communities.

Beyond the species
Why do we need to kill the species concept in the bacterial world?
The role of the biologist is not only to provide data for modelling life but also to define the "living entities" that are its cornerstone. Based on the definition of ecology: "scientific analysis and study of interactions among organisms and their environment", microbial ecology claims similar approaches to those for macro ecology. However, 'organism' (the "living entities") does not have the same meaning in micro-and macro-ecology. In macro-ecology, organisms are pluricellular and are mainly associated with defined species. Thus, microbial ecologists have applied the same reasoning to microorganisms. However, the transfer of the macro-organism model for micro-organisms, whilst trivial to implement, is not realistic because the criteria for species definition do not exist within the micro-organism domain.
Bacteria have a clonal reproduction, whereas sexual reproduction defines the boundary for species in macro-organisms. Moreover, the rate of mutation and broad DNA transfer prevent stable conservation of genetic information. As an example, the size of the pan-genome in the 'defined Escherichia coli species' linearly correlates to the number of genomes sequences (2000 in January 2015) and reaches 90,000 unique gene families. Whereas the corresponding core genome has only 3188 gene families and has not changed much since 2012 (Land et al., 2015).
Despite this, a lot of energy, almost all the energy of microbiologists, was used for: identification of species, classification of species, and function of species, and for microbial ecologists: species dynamics and functions expressed by the different species. With new sequencing techniques, the amount of available data has increased tremendously, but the questions remain roughly the same; identification of species, species dynamics and functions expressed by the different species. The definition of the molecular bacteria species switches from phenotype to a percentage similarity between their 16S rRNA gene, and more recently to co-occurrence of genes (Kim and Price, 2011). Thus, nowadays, the vision of each microbial ecosystem contains more than thousands of different 'species' and for most of them unknown associated functions. Mathematical models in microbial ecology must be developed to cope with this falsely defined artificial diversity.
To highlight the opinion that "bacterial species" are probably not a good entity, consider the classical chemostat model with an abundance of bacterial species. The exclusion principle generally applies and all but one species become extinct. However, if we replace a species (X b ) by the pair it forms with some specific virus (X v ) and consider the same modelling assumptions, then for this three trophic chemostat model (S → X b → X v ), it has been shown that coexistence is possible (Wolkowicz, 1989). Despite this, the model is not suitable with respect to the present discussion and many other aspects, where it can be considered the focus on bacterial diversity is not pertinent.
The opinion in this paper is that the prevalence of this incorrect concept sterilises our understanding of microbial ecology and the modelling approaches. Thus, the first step in a new way of thinking is to kill the concept of species and thereafter to envisage a new paradigm. The new challenge is to identify the entities that really exist and can be "individualised". These entities must make sense in terms of ecology or, in other words, to understand and describe their coupled interactions.
The candidate entities should be the individual cell (individual) or the entire system biocoenosis, but maybe also the expressed functions through genes, proteins or metabolites. Such entities are not so tangible or easy to comprehend as species, and lead to a less intuitive understanding of microbial ecosystems. Acceptance of these new entities will lead to somewhat of a renunciation of the current definition and understanding of diversity in the microbial world.

Conclusions
In this manuscript, we have presented a number of mathematical modelling methods that are relevant to the field of microbial ecology. In terms of applications, we have focused on scientific areas where it is hoped that we will be able, in the future, to overcome particular challenges related to understanding, observation and control, specifically with regard to ecosystem orientation of the metabolic pathways of species within it. We have presented several examples and show how the use of discrete and relatively simple continuous dynamic systems make it possible to deduce new knowledge or identify assumptions about the function of ecosystems. The advantages and limitations of these approaches were presented, accordingly.
Many have predicted that the twenty-first century would be the century of biology (Venter and Cohen, 2004). But the rules of life still remain unknown in many respects and several schools of thought have opposing views to the consensus within the scientific community. Supporters of studying organisms physiologically, as was the case during the 19th and 20th centuries, i.e., by reductionist approaches, face having to handle the complexity revealed by modern microbiological methods. Proponents of the study of natural systems additionally have the task of collating and absorbing the large amount of data that is generated, without the means to handle a lack of appropriate concepts or theories.
Mathematical modelling is a diverse field with form and function dependent on the area to which it is applied and the objectives pursued by the user. We are at a key period of its development where, on the one hand, researchers from very different scientific fields are able to meet and exchange ideas about the components of their study, and on the other, technical and digital advances have provided enormous possibilities for pursuing, often independently, much greater insights into microbial systems. Some dream of new models incorporating large degrees of complexity that would allow the behaviour of their objects of study to be perfectly predicted, without realising that such goals may obfuscate the ability to glean any knowledge of the object itself. Alternatively, those studying simplistic models and their variants ad infinitum, due to the inherent difficulties in analysing more complex models, face a further challenge of the practicality of these models for real-world application. Finally, there are others who believe that the acquisition of ever more data will provide a better understanding of ecosystems, however, the rate of accumulation (e.g., sequence data) puts a greater burden on storage systems, such that these become a limiting factor. More than ever, the parsimony principle should guide the search for the most relevant concepts to advance research in mathematical modelling of microbial ecosystems. Program, July-December 2014, CIB/EPFL, Lausanne, Switzerland. The authors gratefully acknowledge the support of CIB/EPFL in facilitating the workshop that has resulted in this manuscript. M.J. Wade acknowledges the support provided by the Biotechnology and Biological Sciences Research Council UK (BB/K003240/2 Engineering synthetic microbial communities for biomethane production).