A Computational Study of Amensalistic Control of Listeria monocytogenes by Lactococcus lactis under Nutrient Rich Conditions in a Chemostat Setting.

We study a previously introduced mathematical model of amensalistic control of the foodborne pathogen Listeria monocytogenes by the generally regarded as safe lactic acid bacteria Lactococcus lactis in a chemostat setting under nutrient rich growth conditions. The control agent produces lactic acids and thus affects pH in the environment such that it becomes detrimental to the pathogen while it is much more tolerant to these self-inflicted environmental changes itself. The mathematical model consists of five nonlinear ordinary differential equations for both bacterial species, the concentration of lactic acids, the pH and malate. The model is algebraically too involved to allow a comprehensive, rigorous qualitative analysis. Therefore, we conduct a computational study. Our results imply that depending on the growth characteristics of the medium in which the bacteria are cultured, the pathogen can survive in an intermediate flow regime but will be eradicated for slower flow rates and washed out for higher flow rates.

1. Introduction

Safe food is a major contributing factor to human health. The control of foodborne pathogenic bacteria plays an important role in food safety, for example by preservation of foods. Among the traditional methods of food preservation are fermentation, temperature treatment (heat/cold/freezing), or the addition of chemical adjuncts to foods, including spices [1]. A newer concept of food preservation or microbial food safety control could be described as “ecological control”, in which the microbial ecology of the food is augmented by bacteria that are generally regarded as safe (GRAS), with the goal to either out-compete the pathogen, or to change the local environmental conditions such that they become detrimental for the pathogen. Examples for this concept can be found in [1,2,3,4] and the references that they cite. Since this approach is based on the abilities of living microorganisms, it belongs to the so-called biopreservation techniques, in which in particular lactic acid bacteria play an important role [1]. Other biopreservation techniques include control based on microbially produced bacteriocins, such as natural antimicrobials, or so called “hidden fermentation”.

We study here a mathematical model of biopreservation, or biocontrol, of the pathogenic Listeria monocytogenes by the lactic acid bacteria Lactococus lactis, based on earlier work in [2]. The use of such an ecological control mechanism has been proposed for minimally processed refrigerated vegetable products [1,3].

L. monocytogenes is a pathogenic bacteria, which can cause the sometimes fatal disease listeriosis. While individuals with normal health may not develop symptoms, it can be deadly for fetuses, newborns, infants, the elderly, pregnant women, and immuno-compromised people [5]. An outbreak of listeriosis in Canada in 2008 that was traced back to a meat processing facility caused several deaths [6]. L. monocytogenes can be found in uncooked meat, milk, egg, seafoods, and fish as well as certain minimally refrigerated vegetables, but also heat-processed foods such as pasteurized milk and dairy products and ready-to-eat meat [5]. In January through May 2016, the Canadian Food Inspection Agency reported at least 20 food recalls and/or alerts because of this pathogen, the majority of which classified as Class I (or high risk) [7]. Due to its importance from a food safety perspecitve, mathematical modeling to characterize the growth L. monocytogenes has been a very active field of research for more than 20 years. Most studies focus on mono-culture settings, where a wide range of modeling tools is used, from data-driven descriptive fitting, to agent based models, to dynamic differential equations; some such examples are [8,9,10,11,12,13,14,15,16,17], but this list is by no means complete. Fewer models have been proposed to study the dynamics of interaction of L. monocytogenes in co-culture with other species [2,18,19,20], which is the setting that is relevant for the ecological control aspect that we focus on.

The control agent L. lactis is a lactic acid bacteria that is associated with certain dairy products [1]. While it plays an important role in food microbiology because of its fermentation abilities [5], we are here interested in the changes that it inflicts to the microbial ecology by production of lactic acids and subsequent changes of pH. Simply described, the underlying control mechanism is the following amensalistic principle: the control agent produces lactic acids which decrease the pH in the system. Both, decreasing pH and increasing lactic acid concentration are detrimental to the pathogen, while the control agent itself is much more tolerant to these self-inflicted environmental changes. In the mathematical formulation, this is described by a nonlinear system of five ordinary differential equations [2] for the dependent variables population sizes of pathogen and control agent, concentration of lactic acids, concentration of hydrogen ions (as a measure of pH), and malate. Unlike [2], where the model was studied quantitatively with relatively few computer simulations for a batch culture, we investigate it in a continuous chemostat setup aiming at a qualitative description of the longterm dynamics and how it depends on environmental conditions. The mathematical theory for the chemostat, in which several species compete for a limited number of resources, such as nutrients, is well understood, culminating in a mathematically rigorous formulation of the principle of competitive exclusion [21]. However, if additional effects of interaction between species are at play that can give one competitor an additional growth advantage over the other one, this result cannot be readily applied, because it makes heavy use of a certain functional relationship between microbial growth and substrate depletion that no longer holds. Such systems need to be studied individually, cf. [20,22,23,24]. This is also the situation for the amensalistic control system at hand. Moreover, while the principle of competitive exclusion implies that only one of two competitors for a shared resource in a competition based chemostat can survive, it is known from a variety of examples that pH values and dilution rate can determine the fate of co-cultures in chemostats, and that indeed both competitors can survive, cf. [25] for an experimental study. In our specific example, the pH is actually modified by the species themselves. In our model we implicitly assume that nutrient availability is not limiting bacterial growth. This assumption is made, because the set of reaction parameters that we use in our study was determined in [2,14] from experiments in vegetable broth, where nutrients are plenty. Including nutrient limitation effects and competition for nutrients in the model would be straightforward, but at the expense of introducing additional unknown model parameters and dependent variables, thus increasing the model complexity, and potentially shifting the focus away from amensalistic effects.

The algebraic expressions arising from our model are too involved to allow for a complete and rigorous analytical treatment. Therefore, we study the system with a mix of computational methods and analytical techniques, the latter restricted to relative simple special questions, such as stability of trivial equilibria, or to derive estimates on critical flow rates. In our model, the reaction terms contain in total 18 parameters. A complete set of values for these are given in [2]. With growth kinetics fixed, we will focus on the influence of the reactor operating conditions. These are the flow rate and the properties of the growth medium. In our case these are expressed in terms of the bulk concentrations of lactic acids and hydrogen ions. Thus, we explore the behavior of the system with respect to three parameters.

We end this brief introduction with two remarks. First, lactic acid bacteria play also a major role in fermentation. Fermentation is an important, by some accounts the most important, mechanism for food preservation. This is not addressed by our study, which solely focuses on the control of pathogens; Secondly, the biocontrol mechanism that we study is based on a simple ecological (amensalistic) principle: a (possibly invasive) species alters the environmental conditions such that they become less favorable and detrimental for a competitor, while the control agent itself is more tolerant toward these self-inflicted changes. This gives a natural growth advantage. Lactic acid bacteria are often also used as probiotics, which are defined to be live food ingredients which confer health benefits to the host if administered in sufficient quantities [26]. The same amensalistic control mechanism that we study here in the context of food preservations can also play a role as a probiotic control mechanism. Although the parameter values likely will be different, the model itself might apply to this situation as well.

2. Governing Equations

In [2] a mathematical model for the control of Listeria monocytogenes by Lactococus lactis was formulated for the case of batch cultures. We will study this system in the chemostat setting. It is formulated for the five dependent variables

: population size of L. monocytogenes,

: population size of L. lactis,

C: concentration of lactic acids,

P: concentration of hydrogen ions,

M: malate concentration.

The hydrogen ion concentration is equivalent to the pH value in the sense that if P is measured in moles.

In the chemostat setup, the dissolved substrates malate, lactic acids and hydrogen ions are added to the reactor at a constant reactor flow rate q, at bulk concentration levels and . The bacterial populations, and the dissolved growth limiting substrates are completely mixed and subject to convective transport into and out of the reactor. Thus the reactor is continuously replenished with fresh medium.

The bacterial populations grow if C and P are small; they decay if one of these concentrations becomes large. Growth and decay phases are separated by extended neutral phases.

Both bacterial species produce lactic acids until a saturation level is reached. Similarly, lactic acids increase the concentration of hydrogen ions (i.e., decrease the pH) until a saturation level is reached. Moreover, hydrogen ions are diminished by L. lactis, during which process also malate decays.

This is modeled by the nonlinear system of five ordinary differential equations

In (1) and (2) the constants are the maximum specific growth rates of L. monocytogenes and L. lactis. In (3), the constants and are the saturation levels for lactic acid production and and the production rates. Similarly in (4), parameter is the saturation level for production of P and the production rate. Constant is the rate at which malate is decayed by L. lactis, while the rate of decay of P in this process is given by .

The growth functions of L. monocytogenes and of L. lactis are defined by minimal inhibitory concentrations (MIC) for growth and metabolism. In [2], based on the earlier work [14], the following piecewise linear functions are suggested for non-negative arguments and

The coefficient functions are the piecewise linear functions where

Here the function H is defined by

These growth functions are sketched in Figure 1. Since L. lactis is more tolerant to high values of C and P than L. monytogenes, we may assume that . Furthermore we may assume that the saturation levels for production of C and P are in the decay range,

Figure 1

Piecewise linear, continuous net growth rate : The population grows, , for small values of C and P and decays, , if either C or P becomes large; the region in between marks the neutral, stationary phase.

The piecewise definition of the reaction kinetics is unusual in the context of Mathematical Biology in general and microbial growth modeling in particular, where normally sufficiently smooth response functions are the assumed. However, because these functional descriptions, along with quantitative parameters were identified from experiments [2,14], we chose to use them as reported in the experimental literature, rather than mollifying them for mathematical convenience, even if this comes at the expense of increased algebraic complexity.

The right hand side of our model (1)–(5) is continuous but not everywhere differentiable, due to the piecewise definition of and . However, the one-sided derivatives of and exist and are bounded from below by negative numbers that depend on , (note that are monotonously decreasing functions). In particular we have, because of , that

Thus a Lipschitz constant can be found such that (1)–(5) satisfies a Lipschitz condition, which guarantees unique solutions to the initial value problem. These solutions are continuous and at least once differentiable.

The solutions to (1)–(5) with non-negative initial data are bounded from above and remain non-negative.

First we remark, by comparison with the trivial solution of (1) and (2), that the bacterial population sizes and are non-negative.

From (5) it follows by comparison, because .

Let . Then it follows from (3), because of , and again by comparison that . Similarly, let to obtain ; moreover, for all sufficiently large t.

From (4) it follows by comparison with that .

Thus we have established that and . This allows us to improve our lower estimates on the population sizes and . More important, however, is that due to (6) we have and , which implies that and are bounded form above. To see this, we assume that this statement does not hold for one of them, say we assume that (the same line of argumentation applies to ). Then due to (3), . Due to (6), this would imply that is eventually declining, which is a contradiction to our assumption. Because of continuity this means that has a maximum, i.e., is bounded. ☐

The non-negativity of P cannot be established with standard arguments without additional restrictions on parameters. Indeed, one can construct non-negative initial data such that P becomes negative. This situation obviously marks a breakdown of the model. Recall, however, that means , wherefore for practical purposes this break down is not relevant. In none of the numerous simulations that we ran for our study did this breakdown situation occur or did solutions come close to this breakdown situation.

The mathematical theory of competition in the chemostat and related laboratory devices, such as gradostats, is well understood [21]. However, the underlying mathematical machinery makes heavy use of certain close relationships between bacterial growth rates and rates of substrate consumptions, which do not hold for the biocontrol model (1)–(5). Therefore, the mathematical results from traditional chemostat analysis cannot be carried over in a straightforward manner.

The model (1)–(5) is a five-dimensional nonlinear autonomous system of differential equations. For all practical purposes models of this type are typically difficult to study with the qualitative methods of dynamics systems theory if they do not have certain useful properties, such as quasi-monotonicity etc. Moreover, our model has 21 parameters, which could be somewhat reduced by 4 after non-dimensionalization. Furthermore, because of the piecewise definition of , , the Jacobian is only piecewise defined and only for and . A rigorous mathematical analysis would require us to distinguish between a multitude of cases, depending on parameters. The resulting algebraic expressions will be too complex to lend themselves to biologically insightful results. Instead, we will embark on studying the model behavior with a mix of analytical arguments and computational techniques.

In [2] a complete set of reaction parameters is given, that was derived from laboratory experiments and model simulations of the competitive growth of L. monocytogenes and L. lactis in vegetable broth. We will use these data in all our simulations. They are summarized in Table 1. Note that the less tolerant pathogen grows faster than the control agent, . This leaves us with four unknown parameters that were introduced by the chemostat description, flow rate q and the bulk concentrations . In our studies we pick the latter to be , and explore the system behavior with respect to the remaining parameters .

Table 1

Reaction parameters used in this study, from [2].

Parameter	Symbol	Unit	L. lactis	L. monocytogenes
Specific growth rate	μ_2,μ_1	h−1	0.1049	0.1471
Protonated acid production rate	δ,γ	Milimoles CFU−1·h−1	1.7 * 10−10	2.95 * 10−10
MIC acid (growth)	k_3,k_7	Milimolar	5.2	4.058
MIC acid (metabolism)	k_4,k_8	Milimolar	8.907	8.908
Maximum acid concentration	k_1,k_2	Milimolar	11.5	11.65
MIC proton ion (growth)	k_5,k_9	Milimolar	10−1.405	10−1.892
MIC proton ion (metabolism)	k_6,k_10	Milimolar	10−1.147	10−1.151
Maximum proton ion concentration	k_11	Milimolar	10−1.12	10−1.132
Malate decay rate	θ	Milimole CFU−1·h−1	1.69 * 10−10	0
Malate utilization rate	κ	Milimole−1	10−5.33	10−5.33
proton concentration change rate	ρ	Moles CFU−1·h−1	10−5.472	10−5.472

We assume for our simulations that the reactor is initially filled with the same liquid that is used as bulk liquid. This fixes the initial data

3. Analysis of the Single Species Pathogen Sub-Model

Before we study the complete dual-species model, in order to gain a better insight into self-limiting growth behavior, we study a simpler single species model, namely system (1)–(5) for . This is the model for the uncontrolled pathogen. It is easy to verify that then for all and the model reduces to

In order to investigate whether species persists or dies out, we investigate the stability of the equilibrium

The Jacobian in this point is obtained as the triangular matrix

The last two eigenvalues and are always negative. The sign of the first eigenvalue depends on the model parameters. We note that for large enough flow rates q the steady state value converges to the bulk concentration , while for small q it converges to the saturation concentration . For realistic parameters we have , i.e., we expect for small enough q that , i.e., the population dies out due to self inhibition at high P values. On the other hand, the concentrations and in the growth medium are usually small enough that . Thus for high enough flow rates the population dies out due to washout.

This can be formulated more precisely by substituting into , which leads to the following quadratic inequality for the persistence of species

From this we find that is positive if and negative if one of the inequalities is reversed.

In Figure 2 we show simulations of (8)–(10), using the parameters from Table 1. The bulk concentration and were chosen small enough to permit growth of the bacteria under ideal conditions. The initial conditions for C and P were chosen as the bulk concentration values, for we chose 107 CFU/mL. The simulations were conducted for various flow rates q, ranging from to . For the highest of these values , whence (12) predicts washout, for the smallest value , whence (12) predicts extinction due to self-inhibition.

Figure 2

Simulation of the single species model (8)–(10) for various flow rates q. Top panel: population size ; Bottom panel: lactic acid concentration and hydrogen ion concentration in the C-P-plane. The vertical and horizontal lines at and and and mark the transition from growth to neutral to decay regimes.

This is confirmed in the top panel of Figure 2. For the largest flow rate, the flow dominates growth and the population size decreases monotonically. On the other hand, for the smallest flow rate, the population initially grows, reaches a plateau and then eventually dies out. This is explained by the bottom panel of Figure 2. For the highest flow rate, remain in the growth region throughout, but growth is dominated by washout. For the smallest flow rate, however, we notice that initially both C and P increase as a consequence of bacterial activity. They first reach and pass through the neutral range and then, as C continues to be produced, reach the decay range. Eventually, when is small enough, the second term in (9) dominates over the first one and C declines. After some time also P starts to decline and the system converges to the equilibrium (11) with .

The behaviour for intermittent values of the flow rate can vary. In all tested cases, the population establishes itself, often after passing through a brief oscillatory phase which occurs when reach the decay range first and then return to growth range as a consequence of a declining bacterial population. The longterm behavior in these simulations is independent of the initial concentration for , although the transients might change.

The case requires some additional explanation. The criterion (12) predicts that the population will eventually establish itself. This is also confirmed in Figure 2. The population increases first, reaches a plateau, decreases, then increases again and starts oscillating, due to C and P oscillating between the growth and decay stage. For all practical purposes, however, we point out that at around the population drops down to 10 CFU/mL, i.e., to a level where the description of the bacterial population by a continuous variable breaks down. Thus, while the numerical simulation confirms the theoretical analysis, the model is not valid anymore, and does not allow a conclusion whether the theoretical result is rigorous from a practical point of view.

4. Numerical Experiments for the Complete Dual Species Biocontrol Model

4.1. Batch Cultures

We start with a simple simulation experiment of batch cultures, i.e., the special case of our model (1)–(5). This is the case that was originally studied in [2]. We will investigate in a simple simulation experiment how adding the control agent L. lactis to the batch culture initially will affect the eradication of the pathogen L. monocytogenes.

A typical simulation is shown in Figure 3 for initial values

Figure 3

Simulation of (1)–(5) with and CFU/mL, mM, mM and mM.

The bacterial growth curves show the typical three stages, a growth phase initially, followed by a stationary phase, and then the decay phase [5]. Since in these simulations , malate is not replenished. Therefore, M is continuously decreasing and eventually vanishes. On the other hand, the concentrations of lactic acids C and of hydrogen ions P increase continuously due to production, and in the absence of washout. Production of lactic acids slows down as the populations vanish and C reaches a plateau at saturation concentration level. Similarly, pH (i.e., P) reaches a stationary value eventually.

In a simple numerical experiment we investigate how the initial bacterial count of the control agent L. lactis affects the decay of the pathogen L. monocytogenes. To this end, we keep the initial data for the same as in the previous simulations, but vary to be for different values . We measure in our simulations the decay time for , which we define as the first time at which the population size falls below one percent of its initial value, i.e.,

In Figure 4 we plot the population size of the pathogen for different amounts of control agent, as well as as a function of the initial population size of the biocontrol agent. While we see the growth, stationary and decay phase in all simulations, it is of particular interest that the duration of the stationary phase seems to be not simply correlated with the initial amount of control agent. For example, for both and we observe slightly longer stationary phases than for . The onset of the stationary phase, however, happens slightly earlier for higher initial populations of control agents, and the population level of the stationary phase clearly decreases as the initial number of control agents increases. Thus in order to keep the pathogen cell count below a certain value above initial conditions, a sufficiently high number of control agents can be added. The latter two results confirm what one might expect intuitively.

Figure 4

Simulation of model (1)–(5) with , for initial data CFU/mL, CFU/mL, mM, mM, mM. The initial amount of control agent is varied by picking different values for k. The left plot shows the population size of L. monocytogenes for different initial population sizes of L. lactis (. In the right figure the decay time for L. monocytogenes is plotted for different initial population sizes of L. lactis ().

Similarly, increasing the number of control agents initially accelerates eradication of the pathogen, as expected. However, the stopping time does not decrease proportionally with the increase in the initial count of control agents. In fact doubling the number of control agents leads to a decrease in stopping time of less than 10%. Different choices of initial data for lead to qualitatively similar results.

4.2. Longterm Behavior of the Full Chemostat Model

As previously discussed, with the parameters in Table 1 and the assumptions on initial conditions (7), along with and mM, we are left with three free parameters that can determine the behavior of the solutions of model (1)–(5), namely . Since our focus will be on survival and persistence of bacterial populations in the chemostat, we can restrict the ranges of and to the range where both species can grow, in order to avoid a trivial outcome. Thus

Characteristic for the chemostat setup is that the bacterial populations wash out if the flow rate q exceeds the growth rates, independent of the initial bacterial concentration. This is expressed by the asymptotic stability of the trivial steady state, which for our system is obtained as with where again for small enough q we find and for large q we have . For complete washout to occur, we require in extension of the single species model above that

Indeed, the sign pattern of the Jacobian in this equilibrium point is then which implies asymptotic stability of . Note that this is in particular true for , due to the monotonicity of and . Therefore, we can restrict ourselves to

The results of three simulations of (1)–(5) with mM, mM and different values for the flow rate q are shown in Figure 5, Figure 6 and Figure 7.

Figure 5

Simulation of model (1)–(5), with and . Note that but .

Figure 6

Simulation of model (1)–(5), with and .

Figure 7

Simulation of model (1)–(5), with and .

In Figure 5, the flow rate q is smaller than but bigger than the growth rates at equilibrium concentrations, and . Thus, q is big enough to induce washout of both species. and are strictly decreasing functions and equilibrium is attained. The malate concentration decreases first, while sufficient is in the system for degradation but then goes back to due to replenishment. Both P and C first increase due to production by the bacteria and then drop to their equilibrium values when the bacterial populations vanish.

In Figure 6, the flow rate is decreased to , i.e., it is smaller than both maximum growth rates. Initially growth conditions are favorable for both species, which leads to growth of and , implying an increase of C and P, as well as a drop in M. While L. lactis is more tolerant than L. monocytogenes, it also has a smaller maximum growth rate, . Thus, eventually, the actual growth rate of the control agent, while still positive, drops below the dilution rate q due to growth in C and P and dies out. After the extinction of , the malate concentration goes back to its bulk value. The system eventually attains a steady state of type with

Thus, the less tolerant but faster growing pathogen survives. In Figure 7 the flow rate is further decreased to . Again, initially both populations begin to grow. Compared to the previous case in Figure 6, the population sizes reach higher level and accordingly also C and P grow faster. Hence, the growth conditions become more unfavorable and lead to the extinction of the less tolerant pathogen L. monocytogenes, while the more tolerant control agent L. lactis has the opportunity to establish itself, despite its lower growth rate but due to the also lower flow rate. With taking a plateau, the malate concentration declines. Eventually the system reaches an equilibrium of the form with

Thus, the control agent inhibits the establishment of the pathogen.

A fourth type of steady state that is admitted by model (1)–(5) is the co-existence equilibrium in which , and ,

In this case

Note that this condition of two equations for two unknowns and given q does not have a unique solution. Since L. lactis is more tolerant than L. monocytogenes, it does not exist if . i.e., it exists at most if the less tolerant pathogen has a higher maximum growth rate than the more tolerant control agent, . If exist that satisfy (15), then there are infinitely many such pairs, thus the calculation of can not be decoupled from the rest of the system when calculating the equilibrium points. One obtains which still depend on the yet undetermined . In principle, these can be determined from plugging and into the piecewise defined equations (15), however, the resulting expression are too unwieldy to be of practical value. The same holds for the Jacobian and its eigenvalues. Its sign pattern is

To the best of our knowledge no sign pattern criterion similar to those proposed in [27,28,29] exists that could be applied to this matrix to reach conclusions about stability or instability. Numerical calculations for some choices of parameters indicate a positive eigenvalue, and thus, instability.

The results in Figure 5, Figure 6 and Figure 7 show that the flow rate q can play a crucial role in (1)–(5), which is a consequence of a delicate balance between (i) replenishment of the reactor with fresh, growth permitting medium; (ii) bacterial population dynamics, and (iii) removal of bacteria from the reactor. Of these three, the first and last clearly are primarily controlled by q in the sense that both effects become stronger as q increases. These results on longterm behaviour are independent of initial bacterial counts, which however could affect transients, i.e., how quickly the equilibrium solutions are approached.

In order to illustrate the dependency of the long term behaviour on the flow rate, we solve (1)–(5) for q varying over the interval , where is chosen large enough to include washout steady states, . In Figure 8 the population sizes and at steady state are plotted in dependence of q. The case is the case that was studied in the previous Section 4.1. As discussed there, the populations grow and induce unfavorably high concentrations of C and P which first lead to self-inhibition and eventually to decay. This is primarily a consequence of the growth medium not being replenished but also of the fact that active cells are not removed, which increases C and thus P. When increasing q to a critical value , the control agent L. lactis survives while the pathogen L. monogytogenes dies out, i.e., an equilibrium of type is reached. Initially the steady state population size of the control agent increases until . This is primarily a consequence of increasing the flow rate at which fresh medium is replenished, i.e., controlled by (i). With the increasing population size also C and P increase; only if q becomes large enough to contribute considerably to the washout of the L. lactis the steady state population size decreases again. Thus, for the longterm behavior is dominated by (iii). As q passes through , the more tolerant but slower growing control agent is not viable anymore; the growth rate falls below q and the control agent is washed out. This allows the less tolerant but faster growing pathogen to establish itself as long as , i.e., an equilibrium of type is attained. Note that also and change with q. Hence, (iii) is stronger than (i). However, due to the dominance of (iii), increasing q further leads to decreasing steady state populations , until complete washout occurs and both populations die out and an equilibrium of type is reached. Equilibria of type are never attained in our simulations. For the transition case from to the eigenvalues of were computed and it was found that one of them is positive, i.e., is unstable.

Figure 8

Population sizes and at steady state for model (1)–(5) for varying flow rate q and .

So far all computations were carried out for more or less arbitrarily fixed small model parameters and . These bulk concentrations of the growth inhibitors were chosen small enough so that it can be anticipated that bacterial communities can establish themselves. The calculations above show that the steady state values explicitly depend on and . Therefore, in order to investigate whether the behavior seen in Figure 8 is generally observable, we repeat the simulations for various bulk concentration values. To this end the region given by and is discretized by a regular grid of dimension 39 × 37 and for each point on this grid the model was solved for 48 different values of q between . Thus in total, model (1)–(5) was solved 69,264 times for different parameter combinations of . In Figure 9 the steady state population sizes and that are reached in these simulations are visualized. The values for the pathogen L. monocytogenes are color coded using a yellow-green map, while the values for the control agent L. lactis are color coded in a blue-scale color map. The behavior that was seen in Figure 8 is only obtained for relatively small values of and : for small values of q the control agent survives while the pathogen cannot establish itself. After q passes through a critical value , which depends on and , the control agent is washed out and the pathogen can establish itself until q is finally large enough to wash it out, too. Note that the washout value for q also depends on , and . If and are chosen larger, the environmental conditions are too unwieldy for the pathogen to survive. For small flow rates the more tolerant control agent can establish itself but it is washed out if q becomes large. The less favorable the bulk conditions are, the smaller is the minimum flow rate that induces washout.

Figure 9

Exploration of the three dimensional parameter space q (x-axis), (y-axis), (z-axis). Shown are the sizes of the microbial populations and at steady state (4 different views of the same simulation data set).

4.3. Continuously Adding Control Agents to Eradicate the Pathogens

The question arises naturally whether the addition of control agents to the system will enable a quicker eradication of the pathogen. To study this question, the model (1)–(5) is slightly modified by adding a supply term to (2). We obtain where the new parameter is the amount of control agents added continuously to the system. The rigorous qualitative analysis of this model is as cumbersome and impractical as it was for the original model (1)–(5). Therefore, we will restrict ourselves a priori to a computational study and report here the results. The interesting question is how the dosage amount affects eradication times of .

Again we carry out the simulations for the case mM, mM as above, but vary q. First we note that the longterm behavior can depend on both q and . For small , the pathogen could not establish itself even in the absence of control agents being added, and this is of course what one expects and finds in the new simulations, cf. Figure 10a. For flow rates the long term behavior depends on and . For small values of coexistence of pathogen and control agent can be observed, see Figure 10b. This corresponds to the previous case where the pathogen only survived. If the value of is sufficiently increased, the pathogen dies completely out, cf. Figure 10c, i.e., it can be controlled.

Figure 10

Population dynamics for the control model (16)–(20) with continuously added control agents. For (a) small values of ; and values for (b) small and (c) larger amount of control agents added; In (d) eradication time is plotted as a function of the dosage value of the control agent for various .

In Figure 10d finally we plot the eradication time as a function of for four different flow rates . We notice that continuously adding a constant amount of control agents becomes more effective the higher the flow rate q, in the sense that for the lowest flow rate the eradication time reduces by approximately 20% while for the highest flow velocity it decreases by approximately 47% if the dosage is increased from the initial population to twice its value . Moreover, the simulations indicate that increasing the dosage more, to values greater than will not lead to considerably shorter eradication times. Eradication times will quantitatively be different for different initial values of but qualitatively these findings carry over. In applications one will try to balance the dosage between as high as necessary to see the desired effect and to keep it as low as possible, for reasons of process performance, economical cost and to avoid negative microbial side effects that might be introduced by the control agent. The latter could be aspects of microbial food safety in the context of food preservation or texture in the context of functional foods where such a control mechanism could be used, e.g., as probiotics in dairy products.

5. Conclusions

Compared to mathematical models of batch cultures, the chemostat setting has the advantage that it allows a focus on the method of interaction between species, minimizing the effect of the initial population sizes. The well established mathematical theory of the chemostat [21] applies in the first place to systems of microbial competition for common but limited resources. While its qualitative results are powerful and can be considered one of the greatest success stories in Mathematical Biology, they are based on certain functional relationships between bacterial growth terms and substrate consumption terms. However, mathematical models of microbial co-culture systems that are based on other ecological interactions than competition for foods may not have these properties. Therefore, the general theory of the chemostat usually cannot be readily applied and each such system must be studied individually, cf. [20,22,23,24,30]. This is also the case for the amensalistic biopreservation model for the control of the foodborne pathogen L. monocytogenes by the lactic acid bacteria L. lactis that is studied in this paper for nutrient rich environments where competition between both species for nutrients is not growth limiting. Due to the algebraic complexity of the five-dimensional system of ordinary differential equations with 21 parameters, it was necessary to conduct the model study computationally.

Our simulations imply that the persistence of the pathogen depends on the characteristics of the growth medium and on the dilution rate. In order for the pathogen to persist, the concentrations of the detrimental substances must be low enough to allow for bacterial growth and the flow rate must be small enough to not washout the cells before they can reproduce. Only in an intermediate range of the flow rate the pathogen can survive. This is due to the more tolerant but slower growing control agent being washed out of the system in this regime. At lower flow rates the pathogen is crowded out by the control agent which is generally regarded as safe for humans. This amensalistic control effectively enhances and accelerates the self-inhibition that we observe in the uncontrolled system. At higher flow velocities, both species are washed out. In practical applications one might not have unhindered control over flow rate and growth medium. In such a case eradication of the pathogen can be achieved by adding control agents to the system. Our analysis of the uncontrolled pathogenic system implies, however, that the overall longterm behavior as briefly described here is actually the result of the intricate interplay of self-inhibition, amensalism, and reactor conditions. In situations where nutrients become limited, these results might change because, depending on parameters, competition and amensalism may affect the dynamics between both species in opposite direction.

To the best of our knowledge this is the first study of this pH based amensalistic biopreservation mechanism. The model that we study here is somewhat idealized in the sense that it considers a completely mixed culture in the chemostat. While this is achievable in laboratory experiments, and while this is a reasonably good approximation of some bioreactors, many bacteria, including the pathogen L. monocytogenes that we consider in our study, in fact grow as bacterial biofilms, i.e., as spatially structured bacterial populations attached to a surface. While the biocontrol mechanism can be formulated in the biofilm context [31,32,33,34], such biofilm models presently cannot be upscaled to the reactor scale but only studied on the meso-scopic biofilm scale.

While out study was cast for L. monocytogenes and L. lactis and while we used functions of interaction and model parameters from the literaure specific to these two species, the qualitative results may apply to amensalistic systems more generally.