Competitive exclusion in a DAE model for microbial electrolysis cells

1.   Introduction

Microbial electrolysis cells (MECs) are an emerging technology that employs microorganisms to recover energy and resources from organic waste ^[1]. Bacteria on an electroactive anode biofilm oxidize biodegradable substrate and transfer electrons, thereby generating electrical current and releasing protons (H$^{+}$) ^[2]. The protons then recombine to form hydrogen at the cathode. A small voltage (0.2–0.8 V) is needed to overcome the thermodynamic barrier, which is much lower than traditional water electrolysis (1.8–3.5 V) and can be supplied by a small solar panel, low-grade heat, or microbial fuel cells (MFCs), all of which can be available onsite ^[1,2]. The gap of energy input between microbial and pure electrochemical electrolysis is provided by chemical energy stored in the organics. While the electroactive bacteria facilitate hydrogen production, methanogenic archaea consume the same substrate to produce methane, a product which is less energy efficient ^[3]. As a result, methanogenesis leads to decreased efficiency of the system. Since methane is suboptimal as an energy source, operators should seek conditions where electroactive bacteria outcompete methanogens. This is the motivation for our investigation of the competitive exclusion principle in MECs. MEC technology has several advantages over other resource recovery and hydrogen production methods. Microbial electrolysis reduces energy use compared to water splitting because some of the energy is derived from embedded energy in the waste biomass ^[4]. MECs are also more efficient than other methods using renewable wastewater, such as fermentative hydrogen production ^[5]. In fact, ^[6] demonstrated up to 96% recovery of the maximum theoretical yield of hydrogen in MECs operated using fermentation effluent.

In our previous work ^[7], we analyzed and validated a regular, semi-explicit, index 1 differential-algebraic equation (DAE) model for a single substrate MEC ^[8]. The DAE system is an extended version of an ordinary differential equation (ODE) model for chemostats, also known as continuous stirred tank reactors (CSTRs). Besides an ODE system that describes the rate of change of concentrations of the microorganism populations, the biodegradable substrate, and an extracellular mediator involved in electron transfer, the system also includes an algebraic constraint that relates electric current through the external circuit to the concentrations of the electroactive bacteria and the mediator molecule. This constraint accounts for voltage losses that occur in practice. This construct had been used previously ^[8] to avoid solving a partial differential equation for the transfer of electric current. A constraint of this form is also commonly used in chemical fuel cell models ^[9]. Our group has demonstrated computationally that transcritical bifurcations in the dilution rate determine whether electroactive bacteria or methanogens or both will survive as stable equilibria ^[7]. The outcome of competition for substrate is a key question because the types of microbes that exist at the system's stable equilibria determine the electric current and hydrogen production rate. Hydrogen production at the stable equilibrium is possible only if electroactive bacteria are present to generate the needed current. Our efforts here provide answers by characterizing stability of equilibria for two versions of the MEC model, without reference to specific parameter values.

Our goal is to build upon extensive mathematical literature on chemostats to characterize stability of equilibria in the MEC model. One of the main conclusions in the chemostat literature is that if one or more microbes can grow at a lower substrate concentration than the others, then there is a globally asymptotically stable equilibrium where only those microbes have nonzero concentration. This phenomenon is often referred to as competitive exclusion and it holds under a variety of model assumptions. Unfortunately, the MEC analysis is complicated by the fact that growth of the mixed culture bacteria is a nonlinear function of two interdependent variables, the concentrations of both substrate and mediator molecules. In spite of this, we demonstrate that competitive exclusion by methanogens is globally asymptotically stable. We also provide additional conditions which are necessary for local asymptotic stability of electroactive-only equilibria in a simplified model. These conditions suggest that the criteria for competitive exclusion by electroactive bacteria are not as straightforward.

1.1.   Model

Versions of the semi-explicit index 1 DAE model for the MEC have been described previously in ^[7,8,10]. Our model differs from ^[8] by not including fermenting microorganisms that convert a complex substrate into a single compound such as acetate. Additionally, this model does not include a separate methanogen-only outer biofilm layer. We consider the following system, extended to include finitely many microbes of each type:

with initial conditions $0 < S(0)\leq S_{0}$, $0 < \sum_{j}X_{m, j}(0)$, $0 < \sum_{j}X_{e, j}(0)$, $0 < M(0) < M_{\text{total}}$, and $0 < I(0)$. The differential equations represent concentrations, so we are only interested in solutions with nonnegative concentrations. We assume that initial substrate concentration, $S(0)$, is less than or equal to the influent concentration, $S_{0}$. We also assume that initial current, $I(0)$, is positive, due to the nonzero electroactive bacteria concentration and a startup operational period. All of the model parameters are positive.

A single substrate with concentration $S$ flows into the tank at constant rate $DS_{0}$, where $D = F/V$ is the flow rate per volume and $S_{0}$ is the fixed influent substrate concentration. Let ${\bf{X}}_{m}\in\mathbb{R}^{n_{m}}$, and ${\bf{X}}_{e}\in\mathbb{R}^{n_{e}}$ represent the concentrations of microorganisms. The anodic biofilm contains $n_{m}$ methanogen species with concentrations $X_{m, j}$ for $j = 1, 2, \dots, n_{m}$, and $n_{e}$ electroactive bacteria species with concentrations $X_{e, j}$ for $j = 1, 2, \dots, n_{e}$. Substrate consumption is proportional to monotonically increasing microbial growth rates, $\mu_{m, j}(S)$ or $\mu_{e, j}(S, M)$, with constants of proportionality $1/y_{m, j}$ or $1/y_{e, j}$, respectively. Each microbe also has a constant decay rate, $K_{d, m, j}$ or $K_{d, e, j}$. $R$ is the ideal gas constant, $T$ is the temperature, $m$ is the number of moles of electrons transferred per mole of mediator, and $F$ is Faraday's constant. $\gamma$ is a shorthand for the mediator molar mass divided by the product of the MEC volume and the number of electrons transferred per mole of mediator, $mV$.

While methanogens only consume the substrate, electroactive bacteria also consume the oxidized form of a mediator molecule, $M_{\text{ox}}$, that is involved in electron transfer. The mediator exists in oxidized and reduced forms, $M_{\text{ox}}$ and $M_{\text{red}}$, respectively, such as NAD$^+$ and NADH ^[10]. The mediator has a constant maximum concentration, $M_{\text{total}} = M_{\text{ox}}+M_{\text{red}}$. Following ^[10], electroactive bacteria are assumed to transfer electrons via oxidation reduction reactions of the form

These reactions are represented in the diagram in Figure 1. The oxidized mediator is replenished at a rate proportional to the electric current in the device, $I$, as the reduced mediator molecules transfer electrons to the anode. $Y_{M}$ is the mediator yield of the reactions (1.6) and (1.7). To simplify notation in the differential equations, we will refer to the oxidized mediator, $M_{\text{ox}}$, simply as $M$.

The electric current is related to hydrogen production ^[7]. This current can be determined by accounting for the voltage losses in the system. Microbial electrolysis is endothermic, so some small external voltage, $E_{\text{applied}}$, is required. The applied voltage may be opposed by some counter-electromotive force, $E_{\text{CEF}}$. There are also several sources of voltage losses that occur in practice. Ohmic losses, $\eta_{\text{ohm}}$, arise from various types of resistance in the circuit. Activation losses, $\eta_{\text{act}}$, arise from the activation energy of the oxidation-reduction reactions occurring in the cell. Concentration losses, $\eta_{conc}$, arise from certain processes that limit the concentration of reactants at the anode and the cathode ^[9,11]. All of these can be expressed in the following electrochemical balance equation,

where subscripts A and C represent the anode and cathode, respectively. Previous models have ignored concentration losses at the cathode, $\eta_{\text{conc, C}}$, because hydrogen molecules should diffuse away from the cathode rapidly. They have also neglected activation losses at the anode, $\eta_{\text{act, A}}$, under the assumption that the MEC operates with higher voltage losses at the cathode. However, these voltage losses could be included by the general constraint in Section 3. Ohmic losses can be calculated from Ohm's law, $\eta_{\text{ohm}}\left({\bf{X}}_{e}, I\right) = R_{\text{int}}({\bf{X}}_{e})I$ where $R_{\text{int}}(\bf{X_{e})}$ is the internal resistance. $R_{\text{int}}(\bf{X_{e})}$ is a decreasing function of the total electroactive bacteria population because less electroactive bacteria is effectively greater resistance in the circuit. Also, $R_{\text{min}}\leq R_{\text{int}}({\bf{X}}_{e})\leq R_{\text{max}}$ with maximum resistance when there are zero electroactive bacteria and minimum resistance as $\sum_j X_{e, j} \to \infty$. Following ^[12], concentration losses at the anode are modeled by the Nernst equation,

Equation (1.9) assumes that the reference reduced mediator concentration is equal to the total extracellular mediator concentration, $M_{\text{total}}$ ^[8,10]. Activation losses at the cathode are calculated using an approximation to the Butler-Volmer equation for the relationship between electric current and potential at an electrode ^[8]. Standard simplifying assumptions are that the reaction occurs in one step and that the symmetry coefficient (or the fraction of activation loss that affects the rate of electrochemical transformation) is $\beta = 0.5$. With these assumptions we can write

where $A_{\text{sur, }A}$ is the anode surface area and $I_{0}$ is the reference exchange current density. See ^[13] for a discussion of approximations to the Butler-Volmer equation.

The following section discusses previous work that has characterized the equilibria of ODE systems resembling equations (1.1) - (1.2) when ${\bf{X}}_{e}\equiv{\bf{0}}$, $M = 0$, and $I = 0$. This paper extends some of those results to analyze local asymptotic stability in a DAE system with additional equations (1.3) - (1.5), representing current production by electroactive bacteria bacteria using an extracellular mediator.

1.2.   Mathematical Background

There is a large body of literature proving that competitive exclusion occurs for ODE models resembling equations (1.1) - (1.2) with ${\bf{X}}_{e}\equiv{\bf{0}}$, $M = 0$, and $I = 0$. Essentially, stable equilibria may exhibit either competitive exclusion (one species remains), coexistence (multiple species remain), or total extinction (no species remain). Competitive exclusion is generic because coexistence requires multiple species to share an identical parameter value and extinction requires all species to be inadequate competitors. ^[14] proved the competitive exclusion principle for an ODE model of chemostats with microbial growth determined by Monod kinetics. ^[15] then provided a more elegant proof using a Lyapunov function to guarantee global stability. These papers showed that, if survival was possible at all, only the microorganism(s) that could grow at the lowest substrate concentration would survive at the stable equilibria. The work also showed that coexistence was only possible if multiple species could grow at the same smallest substrate concentration. ^[16] provided experiments verifying the theory of competitive exclusion. Subsequently, dozens of authors have proven competitive exclusion occurs in chemostats with various special assumptions. See the monograph ^[17] or the more recent paper ^[18] for more details. Besides microbial growth determined by Monod kinetics, we will also clarify our results by focusing on a simplified case of general monotonically increasing growth rates with equal washout rates and a general constraint ^[19].

To explain microbial electrolysis, the DAE model includes a differential equation for an extracellular mediator involved in electron transfer to the anode as well as an algebraic constraint that determines the electric current. The constraint turns this model into a regular, semi-explicit, index 1 DAE system which does not fit into the previous ODE frameworks. In particular, the constraint (1.5) used in ^[7,8] requires a local representation; global results may not be possible unless the electric current constraint can be solved globally for the electrical current, $I$. In this paper, we extend results from the chemostat literature ^[15,19] to analyze local asymptotic stability in the DAE system given by (1.1) - (1.5). The following section provides an overview of the structure of the rest of the paper.

1.3.   Overview

Section 2 reviews asymptotic stability in semi-explicit DAEs, which is essential to the analysis in the following sections. In particular, we review the relationship between local asymptotic stablility and the spectrum of the matrix pencil, as well as LaSalle's invariance principle for global stability. In Section 3, we use the spectrum of the matrix pencil of the DAE to characterize local asymptotic stability of equilibria in a simplified model with 1 species of each type, general monotonically increasing kinetics, equal decay rates, and a general constraint. This reveals that competitive exclusion by electroactive bacteria is not locally asymptotically stable unless the spectrum of the matrix pencil satisfies certain conditions. Section 4 proves that competitive exclusion by methanogens is globally asymptotically stable for the full MEC system (1.1) - (1.5) with multiplicative Monod kinetics, different decay rates, and a constraint based on the Nernst and Butler-Volmer equations. However, the corresponding result for electroactive bacteria is not likely to be true, as illustrated in the simple case in Section 3. Numerical simulations supporting Theorem 2 and Corollary 2.1 are provided in Section 5. The conclusion in Section 6 summarizes results and points out that the conditions in Section 3 can be used to evaluate numerically whether competitive exclusion by electroactive bacteria will be locally asymptotically stable or not. The conclusion also indicates that MEC operators will want to avoid operating conditions where methanogens can survive at the lowest substrate value because those conditions make competitive exclusion by methanogens globally asymptotically stable.

2.   Asymptotic stability in semi-explicit DAEs

Before presenting results on asymptotic stability of MEC equilibria corresponding to extinction and competitive exclusion, we will briefly review methods for determining asymptotic stability in DAEs ^[20,21]. The DAE framework is necessary because the constraint (1.5) does not admit a global solution. The MEC system in equations (1.1) - (1.5) can be represented as a semi-explicit DAE,

where $f:\mathbb{R}^{r}\times\mathbb{R}^{p}\to\mathbb{R}^{r}$ and $g:\mathbb{R}^{r}\times\mathbb{R}^{p}\to\mathbb{R}^{p}$. More generally, (2.1) can be viewed as a quasilinear DAE,

where $A\in C^{2}(W_{0}, \mathbb{R}^{r+p\times r+p})$ and $F\in C^{2}(W_{0}, \mathbb{R}^{r+p})$ on some open set $W_0$. Both of these perspectives will be useful. For quasilinear DAEs (2.2), local asymptotic stability of equilibria can be determined from the spectrum of the matrix pencil, $\{\sigma A(z)-F'(z):\sigma\in\mathbb{C}\}$ ^{[21,22,23,24]}. An equilibrium point $z^{*}$ of a regular DAE is asymptotically stable if $\text{Re}(\sigma) < 0$ for all elements in $\{\sigma\in\mathbb{C}:\ \text{det}{{\left({\sigma}A(z^{*})-F'(z^{*})\right)} = 0\}}$ ^[21]. In the following section, we will use the spectrum of the matrix pencil to analyze local asymptotic stability in a simplified model with 1 species of each type, general monotonically increasing kinetics, equal decay rates, and a general constraint. The analysis shows that competitive exclusion by methanogens is locally asymptotitcally stable, but the corresponding result is not necessarily true for competitive exclusion by electroactive bacteria.

In the following, we assume that there is some open connected set $\Omega\subset\mathbb{R}^{r+p}$ on which $f$ and $g$ are twice continuously differentiable, and $g_{z}$ is nonsingular on $\Omega$.

Condition 2.1. Suppose that, for some open, connected set $\Omega \in \mathbb{R}^{r+p}$, the following assumptions hold:

1. $f, g \in C^2(\Omega)$;

2. $g_z(y, z)$ is nonsingular on $\Omega$.

Points where $g_{z}$ is nonsingular are called regular. Under the assumption of nonsingularity of $g_{z}$ on all of $W_{0}$, (2.1) is a regular DAE, with index 1. By Condition 2.1, the implicit function theorem allows one to describe $g(y, z) = 0$ as $z = \psi(y)$ on some open neighborhood of $y^{*}$ in $\mathbb{R}^{r}$, where $\psi$ is twice continuously differentiable. At least locally near $(y^{*}, z^{*})$, dynamics of the differential $y$-variables can be described by a reduced ODE,

Solutions of (2.3) which satisfy the constraint, $z = \psi(y)$, are solutions of the semi-explicit DAE (2.1) ^[20,21].

We will also use LaSalle's Invariance principle ^[25] to analyze global asymptotic stability. A well known property of regular semi-explicit, index 1 DAEs (2.1) is that they define a smooth vector field on a smooth manifold ^[21].

Theorem 1. (LaSalle's invariance principle) Consider the smooth dynamical system on an $n-$manifold given by $\dot{x} = X(x)$ and let $\Omega$ be a compact set in the manifold that is (positively) invariant under the flow of $X$. Let $V: \Omega \to \mathbb{R}$, $V \geq 0$, be a $C^1$ function such that

in $\Omega$. Then every solution with initial point in $\Omega$ tends asymptotically to largest invariant set in $\Omega$ where $\dot{V}(x) = 0$. If this set is an isolated equilibrium, it is asymptotically stable.

The function $V$ in Theorem 1 is called a Lyapunov function because LaSalle's theorem generalizes one due to Lyapunov where $\dot{V}$ must be strictly less than zero. Section 4 applies a Lyapunov function modified from ^[15] to the regular semi-explicit, index 1 DAE with finitely many species, multiplicative Monod kinetics, different decay rates, and a $C^\infty$ constraint that is not explicitly solvable on $\Omega$.

3.   Local asymptotic stability in a simplified model

In this section, we consider local asymptotic stability of equilibrium points corresponding to extinction and competitive exclusion in a simplified MEC model. In contrast to chemostats, competitive exclusion by electroactive bacteria is not necessarily locally asymptotically stable, even when electroactive bacteria can grow at the lowest substrate concentration. This is due to several issues, including the nonlinear dependence of the growth of electroactive bacteria on both mediator and substrate concentrations, as well as the form of the algebraic constraint which determines the electric current. We present this analysis to promote clarity in in Section 4.

Here we simplify the model by assuming that there is one compartment for each type of microbe (i.e., $n_{m} = n_{e} = 1$), that the biofilm decay rates are equal to the dilution rate (i.e., $K_{d, m, 1} = K_{d, e, 1} = D$). This is a mathematical simplification that allows for the simple representation of the dimensionless system (3.1) - (3.5). This greatly simplifies the presentation of conditions for local asymptotic stability. However, it may not truly represent the biology of the biofilm. We allow for general monotonically increasing kinetics with $\mu_{m, 1}(S)$ and $\mu_{e, 1}(S, M)$. We also allow for a general constraint, $0 = g(X_{e, 1}, M, I)$, with certain reasonable derivative conditions summarized below. These assumptions allow concise conditions for local asymptotic stability. The presentation of results is simplified if we rescale the model variables. We set

This yields a system of the form

where

and $g(x_{e}, m, I)$ represents the rescaled constraint. Derivatives are taken with respect to rescaled time $\tau$. Note that we can form the new variable $u = s+x_{m}+x_{e}$ to obtain the system

Since $u$ must satisfy $u(t) = 1+Ce^{-t}$, solutions of $u(t)$ will approach $1$ as $t\to\infty$. Therefore, any asymptotically stable equilibria will satisfy $u = s+x_{m}+x_{e} = 1$. We focus on equilibria corresponding to extinction, where $s = 1$, and competitive exclusion, where either $s+x_{m} = 1$ or $s+x_{e} = 1$.

In this section, we will allow for general monotonically increasing growth functions and a general constraint. We consider a class of kinetics where the growth rates and substrate consumption rates of the microbes will increase with substrate concentration and also mediator concentration in the case of electroactive bacteria. We require that attainable equilibria exist. We also make several physically reasonable assumptions based on (1.5). These assumptions state that ohmic voltage losses in $g$ will decrease with $x_{e}$ (because a decrease in electroactive bacteria concentration is effectively an increase in resistance in the circuit), that concentration losses will increase with the oxidized mediator $m$ as the reduced mediator becomes limited at the anode, and that activation losses will increase with $I$. Additionally, the absence of electroactive bacteria means that no electric current is present. Finally we require unique solutions to $g$ at two points. The requirements are summarized by the following condition.

Condition 3.1. The kinetics and constraint in this section satisfy:

1. $f_{m}(0) = f_{e}(0, m) = f_e(s, 0) = 0$;

2. $f_{m}(s)$ and $f_{e}(s, m)$ are continuously differentiable and monotonically increasing;

3. $\exists$ $\lambda_{m}$ and $\lambda_e(m) < 1$ such that $f_{m}(\lambda_{m}) = 1$ and $f_e(\lambda_e(m), m) = 1$ for some $m \in (0, m_0]$;

4. $\frac{\partial g}{\partial x_{e}} > 0$ and $\frac{\partial g}{\partial m}, \frac{\partial g}{\partial I} < 0$;

5. $x_e = 0 \Rightarrow I = 0$;

6. $\exists$ a unique positive $m_0 > 0$ such that $0 = g(0, m_0, 0)$;

7. $\exists$ a unique positive $m^* > 0$ such that $0 = g\left(1-\lambda_{e}(m^{*}), m^{*}, \frac{1-\lambda_{e}(m^{*})}{\Gamma}\right)$.

The rescaled model (3.1) - (3.5) has several important equilibrium points corresponding to extinction of all microbes or competitive exclusion by one type of microbe. These depend on the substrate (and mediator) values where each microbe attains zero net growth. $\lambda_{m}$ and $\lambda_{e}(m)$ are the substrate concentrations where $f_{m}\left(s\right) = 1$ and $f_{e}\left(s, m\right) = 1$, respectively. In the case of the electroactive bacteria, there is a curve $s = \lambda_{e}(m)$ that satisfies $f_{e}\left(s, m\right) = 1$. Biologically meaningful solutions will be located in the interval $(0, 1)$, which corresponds to attainable substrate concentrations. If $\forall m$, $\lambda_{m}, \lambda_{e}(m)\leq0$, then no microbe can ever have positive net growth and if $\lambda_{m}, \lambda_{e}(m)\geq1$, then all microbes require more substrate than is available given the influent substrate concentration. The model equilibria exhibiting extinction and competitive exclusion are shown in Table 1. The mediator concentrations in these equilibria are the solutions to the constraint, $0 = \ g(x_{e}, m, I)$, at the corresponding points. In particular, $m_{0}$ is the unique positive solution to $0 = g\left(0, m_{0}, 0\right)$ and $m^{*}$ is the unique positive solution to $0 = g\left(1-\lambda_{e}(m^{*}), m^{*}, \frac{1-\lambda_{e}(m^{*})}{\Gamma}\right)$.

The following mutually exclusive cases make one of the equilibrium points locally asymptotically stable. In the final case, local asymptotic stability of the electroactive-only equilibrium also depends on a discriminant that appears in two elements of the spectrum of the matrix pencil at $p_{e}$; we denote this discriminant by

where $x_{e}^{*} = 1-\lambda_{e}(m^{*})$.

1. (Total Extinction): If $\lambda_{m}, \lambda_{e}(m_{0})\notin(0, 1)$, then $p_{0}$ is locally asymptotically stable.

2. (Competitive Exclusion by Methanogens): If $\lambda_{m} < \lambda_{e}\left(m_{0}\right)$, then $p_{m}$ is locally asymptotically stable.

3. (Competitive Exclusion by Electroactive): If $\lambda_{e}(m^{*}) < \lambda_{m}$, and either (1) $\delta < 0$ or (2) $\text{Re}\left(\sqrt{\delta}\right) < -\left(x_{e}^{*}\frac{\partial f_{e}}{\partial m}\frac{\partial g}{\partial I}+x_{e}^{*}\frac{\partial f_{e}}{\partial s}\frac{\partial g}{\partial I}+\Gamma\frac{\partial g}{\partial m}\right)\big|_{p_{e}}$, then $p_{e}$ is locally asymptotically stable.

The extinction equilibrium, $p_{0}$, will be unstable in general because we expect at least one microbe to be an adequate competitor with a $\lambda_\star$ value in $\left(0, 1\right)$. The spectrum of the matrix pencil at $p_{0}$ is

so $p_{0}$ is unstable as long as either $\lambda_{m}$ or $\lambda_{e}(m_{0})$ are in the interval $\left(0, 1\right)$. In this case, a microbe introduced into the system may grow in the presence of plentiful substrate. If the methanogen can grow at the lowest substrate value, then the corresponding methanogen-only equilibrium, $p_{m}$, will be locally asymptotically stable. The spectrum of the matrix pencil at $p_{m}$ is

If Case 2 holds, then only $x_{m}$ will be able to attain positive net growth near the corresponding equilibrium. Finally, the spectrum of the matrix pencil at $p_{e}$ is $\left\{ \sigma_{1}, \sigma_{2}, \sigma_{3}, \sigma_{4}\right\}$ where

If Case 3 holds, then $p_{e}$ is locally asymptotically stable.

Cases 1 - 3 provide conditions for local asymptotic stability of each of the equilibria points $p_{0}, \ p_{m},$ and $p_{e}$. These apply for any continuously differentiable and monotonically increasing growth functions $f_{m}(s)$ and $f_{e}(s, m)$ and any $g$ that satisfies Condition 3.1. The conditions in Cases 1 - 3 can be checked numerically to determine if a parametrized model permits a locally asymptotically stable equilibrium where only the most competitive electroactive species persists. Unlike the methanogen-only equilibria, local asymptotic stability of the electroactive-only equilibria is not guaranteed when the electroactive bacteria can survive at the lowest substrate value. This means it is unlikely that chemostat results extend to electroactive-only equilibria in an MEC. However, we can assert that for Monod kinetics with different decay rates and a constraint based on the Nernst and Butler-Volmer equations, methanogen-only equilibria are globally asymptotically stable when methanogens can grow at the lowest substrate concentration. The proof of this assertion in Section 4 relies on LaSalle's invariance principle.

4.   Global asymptotic stability with Monod kinetics

In this section, we consider the full MEC system (1.1) - (1.5) with multiplicative Monod kinetics, different decay rates, and a constraint based on the Nernst and Butler-Volmer equations. LaSalle's invariance principle allows us to show that competitive exclusion by methanogens is globally asymptotically stable. The proof uses a Lyapunov function adapted from ^[15]. Suppose that growth rates for the methanogens are

and growth rates for the electroactive bacteria are

$\mu_{\text{max}, m, j}$ and $\mu_{\text{max}, e, j}$ are the maximum growth rates; $K_{S, m, j}$ and $K_{S, e, j}$, are half rate constants for consumption of substrate; $K_{M, j}$ is the half rate constant for consumption of mediator. As before, the equilibria depend on parameters or functions that are the $S$ solutions to $\mu_{m, j}\left(S\right) = K_{d, m, j}$ and $\mu_{e, j}\left(S, M\right) = K_{d, e, j}$. Denote these solutions by

Assuming that the microbe concentrations are not zero, then $\lambda_{m, j}$ and $\lambda_{e, j}(M)$ are the substrate concentrations at which each microbe has zero net growth. The difference between the two types of microorganisms is that each electroactive bacteria has zero net growth on a curve $S = \lambda_{e, j}(M)$. Examples of these curves are shown in Figure 3. The fact that electroactive bacteria have dual substrate-mediator limitation complicates the type of analysis that has appeared in the chemostat literature.

The equilibrium points corresponding to extinction and competitive exclusion in the full unscaled model are given in Table 2. Microbe concentrations at these equilibria are solutions to $\dot{S} = 0$ where only one microbe survives. If $X_{m, 1}$ or $X_{e, 1}$ are the most competitive microbes, i.e., $\lambda_{m, 1}$ or $\lambda_{e, 1}(m_0)$ are the smallest $\lambda$ values, then equilibrium concentrations are

The mediator concentrations $M_0$ and $M^*$ are the solutions to the constraint (1.5) at the corresponding points:

Finally, in the case of competitive exclusion by electroactive bacteria, the electric current is the solution to $\dot{M} = 0$ when ${\bf{X}}_{e}(t) = {\bf{X}}_{e}^{*}$ and $M(t) = M^{*}$:

This analysis focuses on several sets of interest. Let

In practice, $\Omega$ is bounded because the dynamical system is dissipative, as shown in appendix 6. The maximum concentration of each microbe must be bounded because it is not biologically possible to have infinite concentration. Although the upper bound for each microbe concentration is not clear, concentrations of each species will be bounded as $t \to \infty$. Let $\Omega_1 \subset \Omega$ be the bounded set containing these dynamics. Let $G \subset \overline{\Omega}_1$ be the closed set where the $C^\infty$ constraint (1.5) is satisfied. To obtain consistent initial conditions for the DAE, we will assume or the remainder of this section that initial conditions lie in the compact set $\Omega_G := \overline{\Omega}_1 \cap G$.

Condition 4.1. Suppose that $(S(0), {\bf{X}}_{m}(0), {\bf{X}}_{e}(0), M(0), I(0))\in\Omega_G$.

Our first theorem relies on the following lemma regarding positivity and boundedness of the DAE solutions.

Lemma 1. $\Omega_G$ is positively invariant for (1.1) - (1.5) with Monod kinetics given by (4.1) - (4.2).

We defer the proof of Lemma 1 to appendix 6. Lemma 1 will be used in the proofs of the theorem later in this section. The next lemma identifies conditions under which a microorganism cannot survive at the stable equilibrium.

Lemma 2. If a microbe species cannot obtain zero net growth for $\left(S, M\right)$ values in $(0, S_{0}]\times(0, M_{0}]$, then the concentration of that species will go to zero as $t\to\infty$.

We leave the proof in appendix 6. The intuition behind Lemma 2 is that the substrate concentrations where each microbe has zero net growth are $S$-coordinates of equilibria points and they must be attainable in the interval $(0, S_{0}]$. For electroactive bacteria, the equilibrium substrate concentrations must be in the interval $\left(0, S_{0}\right]$ for obtainable mediator concentrations, $M\in\left(0, M_{0}\right]$. If $\lambda_{m, j}$ and $\lambda_{e, j}(M)$ are not in this interval, definitions (4.3) and (4.4) tell us that either (a) the maximum growth rate is less than or equal to the decay rate, or (b) the microbe requires more substrate than is flowing into the device. In other words, for the microorganisms to survive, they must be able to attain positive net growth and must not require more substrate than is available. We will assume without loss of generality that the following condition holds for each microbe species; otherwise, the corresponding concentration will approach zero concentration as $t\rightarrow\infty$.

Condition 4.2. Suppose that $\lambda_{m, j}\in\left(0, S_0\right]$, and ${\lambda_{e, j}(M)}\in(0, S_0]$ for some $M \in (0, M_0]$.

We now present the main result of this section, a theorem describing the competitive exclusion principle in the MEC. Limiting behavior of the DAE system (1.1) - (1.5) with Monod kinetics (4.1) - (4.2) is determined by the smallest element in $\Lambda:=\left\{ \lambda_{m, j}\right\} {}_{j = 1}^{n_{m}}\cup\left\{ \lambda_{e, j}(M_{0})\right\} _{j = 1}^{n_{e}}$, the set of smallest substrate concentrations where each microbe has zero net growth. Intuitively, when $S$ approaches $\min\left(\Lambda\right)$ from above, all but the most competitive microbes will have negative net growth.

Theorem 2. (Competitive Exclusion by Methanogens) Suppose that Conditions 3 and 4 hold. Suppose also that $\lambda_{m, 1}$ is strictly smaller than all other elements of $\Lambda$ (i.e., methanogen $X_{m, 1}$ can grow at the lowest substrate concentration). Then all solutions of (1.1) - (1.5) with Monod kinetics (4.1) - (4.2) will approach the point $P_{m}$ as $t\to\infty$.

Proof. (^[15] provided a Lyapunov function for equations (1.1) - (1.2) with ${\bf{X}}_{e} = {\bf{0}}$. That function is extended to provide these results.) $\Omega_G\subset\mathbb{R}^{n_{m}+n_{e}+3}$ is compact and, by Lemma 1, it is positively invariant for (1.1) - (1.5). Suppose that $\lambda_{m, 1}$ is strictly smaller than all other elements of $\Lambda$. Let

for some unspecified constants $c_{m, j}$, $c_{m_{2}, j}$, and $c_{e, j}$. Then

To show that $\dot{V}\leq0$ in $\Omega_{G}$, we will use the following rearrangements of equations (1.2) - (1.3):

Choose positive constants

Then

Since

we can combine corresponding sums to obtain

By assumption, each of the sums is less than or equal to zero. Substituting $X_{m, 1}^{*} = \frac{D\left(S_{0}-\lambda_{m, 1}\right)y_{m, 1}}{\mu_{m, 1}\left(\lambda_{m, 1}\right)}$ yields

Thus, Theorem 1 tells us that solutions to equations (1.1) - (1.5) that start in $\Omega_{G}$ approach the largest invariant set in

From equation (1.4), we know that $\dot{M} = \gamma I\geq0$ in $E$. Thus, solutions approach $M = M_{0}$ and $I = 0$ or

Theorem 2 indicates that the methanogen that can grow at the lowest substrate concentration will outcompete as $t\to\infty$. It is possible, albeit unlikely, that multiple methanogen species obtain zero net growth at the same smallest substrate concentration. The following corollary predicts that all of these methanogens will coexist while competitively excluding the other microbes as $t\to\infty$.

Corollary 2.1. (Coexistence of Methanogens) Suppose that Conditions 3 and 4 hold and that $\lambda = \lambda_{m, 1} = \dots\lambda_{m, k}$ is strictly smaller than other distinct elements of $\Lambda$. Then as $t\to\infty$, all solutions of (1.1) - (1.5) with Monod kinetics (4.1) - (4.2) will approach the invariant set

Proof. If $\lambda = \lambda_{m, 1} = \dots\lambda_{m, k}$ are strictly smaller than all distinct elements of $\Lambda$, then terms corresponding to $X_{m, 1}, \dots, X_{m, k}$ vanish from $\dot{V}$ in equation (4.6) in the proof of Theorem 2. We still have $\dot{V}\leq0$, but now solutions approach the largest invariant set in $E = \{(x, y)\in\Omega_{G}:\dot{V}(x, y) = 0\} = \{(\lambda, \tilde{{\bf{X}}}_{m}, {\bf{0}}, M, I)\}$. As before, $M \to M_0$ and $I \to 0$ in $E$, so solutions the largest invariant set approach $W = \{(\lambda, \tilde{{\bf{X}}}_{m}, {\bf{0}}, M_{0}, 0)\}$.

Theorem 2 and Corollary 2.1 predict that the methanogen(s) that can grow at the smallest substrate concentration will outcompete the other species. It is possible that multiple species will share the same smallest $\lambda$ value, in which case coexistence of methanogens is possible. Theorem 2 also tells us that electroactive bacteria are guaranteed to lose the competition if $\forall j$, $\lambda_{e, j}(M_{0}) > \min\left(\Lambda\right)$. Unfortunately, the Lyapunov function used in the proof of Theorem 2 does not suffice to prove analogous results about global asymptotic stability for electroactive-only equilibria. As shown in Section 3, it is likely that additional conditions on the spectrum are required to reach conclusions about global stability of competive exclusion by electroactive bacteria. Unfortunately, due to the constraint (1.5), the spectrum of (1.1) - (1.5) cannot readily be evaluated, even when there is only one species of each type.

The chemostat literature has shown that competitive exclusion by the microbe that survives at the lowest substrate concentration is globally asymptotically stable in a variety of cases. Several authors have considered limitation by two complementary or substitutable substrates ^{[26,27,28,29]}, using either monotone or minimum Monod kinetics. To our knowledge, no one has considered multiplicative Monod kinetics of the form given by (4.2), particularly when one of the limiting substrates is a mediator molecule whose concentration depends on an algebraic constraint. Figures 3 and 4 show that the behavior of the electroactive bacteria is more complicated than a microbe growing on one or two substrates. One also needs detailed information about mediator concentration because electroactive bacteria grow much more slowly when the mediator concentration is low. Given the results in Section 3, it is not surprising that global asymptotic stability of competitive exclusion by electroactive bacteria is more complicated. In the next section, we provide numerical simulations supporting Theorem 2 and Corollary 2.1. These simulations also show that competitive exclusion by electroactive bacteria may occur if electroactive bacteria can grow at the lowest substrate concentration.

5.   Numerical Simulations

In this section, we consider solutions of (1.1) - (1.5) with Monod kinetics (4.1) - (4.2) and with one species of each type. We demonstrate that when Theorem 2 and Corollary 2.1 are satisfied, then solutions behave as expected. That is, if one or more species of methanogens can survive at the lowest substrate concentration, then the model exhibits competitive exclusion by those methanogens. We also demonstrate that if electroactive bacteria survive at the lowest substrate concentration, then the model may exhibit competitive exclusion by electroactive bacteria.

For the simulations below, we generally use parameters from the first table of ^[7], with the following exceptions. The influent substrate concentration is $S_{0} = 100$ and the maximum substrate consumption rates for each species are $14$. The maximum growth rates vary between simulations to change which microbe can grow at the lowest substrate concentration; the parameters are given in Table 3. The new parameters were chosen to demonstrate competitive exclusion by methanogen 1, coexistence two of methanogens, and competitive exclusion by electroactive bacteria. Numerical solutions were generated using the variable-order, variable-coefficient backward differentiation formula in fixed-leading coefficient form ^[30] from the IDAS package of SUNDIALS suite of nonlinear and DAE solvers ^[31]. Initial conditions were set as $S(0) = 100$, $X_{m, 1}(0) = 1$, $X_{m_{2}, 1}(0) = 1$, $M(0) = 25$, and $I(0) = 6$ while $X_{e, 1}(0)$ is the solution to the algebraic constraint (1.5). Figure 4a demonstrates that if a methanogen can grow at the lowest substrate concentration, then all solutions converge to a methanogen-only equilibrium, $P_{m}$. Figure 4b shows that if multiple methanogens can survive at the same substrate concentration, then solutions converge to a set with only those methanogens. These simulations support Theorem 2 and Corollary 2.1. Figure 4c shows that if an electroactive bacteria can grow at the lowest substrate concentration, then solutions may converge to an electroactive-only equilibrium.

6.   Conclusion

In Section 3, we characterized local asymptotic stability of equilibria in a model with one species of each type, equal decay rates, general monotone kinetics, and a general constraint. Subsequently, in Section 4, we showed that competitive exclusion by methanogens is globally asymptotically stable in a model with finitely many species, multiplicative Monod kinetics, different decay rates, and a constraint based on the Nernst and Butler-Volmer equations. Our results also show that certain operating conditions should be avoided. In both models, if a methanogen species can grow at the lowest substrate value, then competitive exclusion by methanogens is either locally or globally asymptotically stable. These results on competitive exclusion provide a recipe for MEC operation that offers the best chance for long term electrical current and hydrogen production:

1. Determine which microbe can grow at the lowest substrate concentration to ensure that methanogens will not outcompete eventually. Theorem 2 from Section 4 indicates that methanogens will competitively exclude the other microbe species if they can survive at the lowest substrate concentration.

2. Compute the spectrum of the matrix pencil at the electroactive-only equilibrium to ensure that it is locally asymptotically stable. Recall from Case 3 in Section 3 that if one of the electroactive bacteria has zero net growth at the lowest substrate concentration, $\lambda_{e, 1}(m^{*})$, and that the discriminant (3.6) satisfies either $\delta < 0$ or $\text{Re}\left(\sqrt{\delta}\right) < -\left(x_{e}^{*}\frac{\partial f_{e}}{\partial m}\frac{\partial g}{\partial I}+x_{e}^{*}\frac{\partial f_{e}}{\partial s}\frac{\partial g}{\partial I}+\Gamma\frac{\partial g}{\partial m}\right)\big|_{p_{e}}$, then competitive exclusion by electroactive bacteria is locally asymptotically stable. Then solutions near an electroactive-only equilibrium will approach that equilibrium.

In summary, if electroactive bacteria can survive at the lowest substrate concentration and the technical condition on the discriminant of the matrix pencil is satisfied, then electroactive bacteria are most likely to outcompete methanogens and the MEC is most likely to provide long term electrical current and hydrogen production. On the other hand, if the less energy efficient methane production is desirable, one should guarantee that a methanogen can grow at the lowest substrate concentration; then the model will exhibit competitive exclusion by methanogens. These results provide insight about whether the microbial electrolysis system will produce methane or electric current and resulting hydrogen.

Acknowledgments

We would like to thank Patrick De Leenheer, Jim Meiss, and Juan Restrepo for their valuable feedback and insight.

Conflict of interest

The authors declare that they have no competing interests.

Appendix

1. Proof of Lemma 1

Proof. (See ^[32] for elements of the proof pertaining to $S$ and $X_{\star, j}$.) $S(t)$ is positive because $S(\tau) = 0 \Rightarrow\dot{S}(\tau) = DS_{0} > 0$. For $\star \in \{m, e\},$ each $X_{\star, j}(t)$ is positive because boundaries where $X_{\star, j} = 0$ are invariant and cannot be reached in finite time if $X_{\star, j}(0)$ is positive. $S(t)$ and $X_{\star, j}(t)$ coordinates are bounded because of the following. Define

Let $\bar{D}$ be the minimum of $D$ and each $K_{d, *, j}$. Then

Thus, $S$ and each $X_{\star, j}$ is positive and bounded. In fact, $S(t) \leq S_0$ because $S(t^*) = S_{0} \Rightarrow \dot{S}(t^*)\leq0$.

$M$ and $I$ are positive and bounded because of the algebraic constraint (1.5). $M$ cannot be negative because $M = 0\Rightarrow\dot{M} = \gamma I\geq0$. Exponentiate both sides of (1.5) and solve for $M$ to obtain

Since the exponential is positive, $M < M_{\text{total}}$. In fact, there is an even tighter bound. Define

From equation (6.1),

Also, equation (1.4) tells us that $I = 0\Longleftrightarrow\dot{M}\leq0$. Thus, $I(0)\geq0\Rightarrow M(0)\leq M_{0}\Rightarrow M(t)\leq M_{0}$.

$I(t)$ must be positive because $0 < M(t) < M_{0}$. It remains to show that $I(t)$ is bounded. Because $M(t) < M_{0}$ we know that

or

Since $\text{arcsinh}(\frac{I}{2I_{0}}) > 0$ when $I > 0$ and $R_{\text{min}}\leq R_{\text{int}}({\bf{X}}_{e})\; \leq R_{\text{max}}$, it must be true that

Thus, $M$ and $I$ are positive and bounded. Let $\Omega_1$ be the set where these variables are bounded. To be consistent with the $C^\infty$ constraint (1.5), solutions of the DAE must lie in the closed set $G \subset \overline{\Omega}_1$. Thus, the compact set $\Omega_G := \overline{\Omega}_1 \cap G$ is positively invariant for (1.1) - (1.5).

2. Proof of Lemma 2

Proof. (See ^[14] for the proof when ${\bf{X}}_{e}(t) = {\bf{0}}$ and there is no mediator or current.) If $\lambda_{m, j}\leq0$, then $\mu_{\text{max}, m, j}\leq K_{d, m, j}$. We can rearrange the integral representation of $X_{m, j}$ and use the fact that $S(t)\leq S_{0}$ to get

$X_{m, j}(t)$ is positive and bounded by a decaying exponential, so $\lim_{t\rightarrow\infty}X_{m, j}(t) = 0$.

If $\lambda_{m, j} > S_{0}$, then $\mu_{\text{max}, m, j} > K_{d, m, j}$. We can rearrange the integral representation of $X_{m, j}$ and use the same fact as before to get

$X_{i, j}(t)$ is positive and bounded by a decaying exponential, so $\lim_{t\rightarrow\infty}X_{m, j}(t) = 0$.

Now consider $X_{e, j}$. If $\lambda_{e, j}(M)\leq0$ for all $M\in(0, M_{0})$, then $\mu_{\text{max}, e, j}\left(\frac{M}{K_{M, j}+M}\right)\leq K_{d, e, j}$ and $\mu_{e, j}(S, M)\leq K_{d, e, j}$ for all $S\in\left(0, S_{0}\right]$ and $M\in(0, M_{0})$. Thus,

$X_{e, j}(t)$ is positive and bounded by a decaying exponential, so $\lim_{t\rightarrow\infty}X_{e, j}(t) = 0$. If instead $\lambda_{e, j}(M) > S_{0}$ for all $M\in(0, M_{0})$, then $\mu_{\text{max}, e, j}\left(\frac{M}{K_{M, j}+M}\right) > K_{d, e, j}$ and $S_{0} < \lambda_{e, j}(M_{0})$. Rearrange the integral representation to find

$X_{e, j}(t)$ is positive and bounded by a decaying exponential, so $\lim_{t\rightarrow\infty}X_{e, j}(t) = 0$.