Quantifying Geobacter sulfurreducens growth: A mathematical model based on acetate concentration as an oxidizing substrate

1.   Introduction

Electrogenic bacteria, a subset of prokaryotic microorganisms, possess the unique ability to donate electrons to a conductive surface. Within this category, Geobacter stands out for its proficiency in electron transfer to molecules like Fe(Ⅲ) and Mn(Ⅳ), nitrogen oxides, and high molecular weight substances that result from the chemical and biological degradation of plant and animal remains ^[1]. Geobacter also exhibits the remarkable ability to directly transfer electrons to electrodes. In this process, the bacterium adheres to the conductive surface (electrode), channeling electrons to the electrode through its cytochromes or pili, which function as nanowires, thereby generating small electric currents ^[2]. The electrogenic property of Geobacter finds practical application in microbial fuel cells (MFCs). These small reactors convert chemical energy into electrical energy through the process of microbial communities catalyzing chemical reactions ^[3]. MFCs involve a complex interplay of physical, chemical, and biological components, influencing the metabolic functions of electrogenic bacteria and the interactions and properties of individual substrates within the MFC ^[4].

Geobacter sulfurreducens, initially identified in 1993 by Dr. Derek Lovley in freshwater sediment from the Potomac River in Washington USA, belongs to the family Geobacteraceae. This family is situated within the domain of bacteria, the phylum Proteobacteria, the class Deltaproteobacteria, and the order Desulfuromonadales. Characterized as Gram-negative, anaerobic, and chemoorganotrophic organisms ^[1], the Geobacter species are noteworthy components of microbial communities.

The two most extensively studied Geobacter species are G. metallireducens and G. sulfurreducens. G. metallireducens, the first Geobacter species to be isolated in pure culture, played a pivotal role in elucidating the genetic system of Geobacter and garnered attention for its exploration of anaerobic degradation pathways. In contrast, G. sulfurreducens, renowned as a sulfur reducer, holds the distinction of being the first genetically manipulated Geobacter species ^{[5,6,7,8,9,10,11,12,13,14,15]}. Since then, G. sulfurreducens has been a focal point in genomic studies, contributing significantly to our understanding of the metabolism and distinctive traits within the genus Geobacter.

All Geobacter species exhibit Gram-negative rod morphology and possess the ability to oxidize acetate while reducing Fe(Ⅲ). These bacteria commonly demonstrate additional traits, including the reduction of Mn(Ⅳ), U(Ⅵ), elemental sulfur, and humic substances. Moreover, many Geobacter strains display versatility by utilizing other low molecular weight organic acids and ethanol as electron donors ^[16]. Specifically, G. sulfurreducens primarily utilizes acetate as its main electron donor, although studies suggest its potential utilization of hydrogen (i.e., H₂). Notably, among the electron acceptors, fumarate plays a significant role alongside Fe(Ⅲ). Biochemical investigations have identified proteins with fumarate reductase activities in G. sulfurreducens, indicating a dual function that is similar to that observed for succinate dehydrogenase ^[17].

Geobacter species, exemplified by G. sulfurreducens, find diverse applications, particularly in environmental remediation and bioelectronics. Notably, these bacteria are recommended for their high functionality in sludges and water that have been contaminated with aromatic compounds, hydrocarbons, and metals, catalyzing their conversion to a more manageable form, i.e., both soluble and precipitate. Geobacter plays a pivotal role as a key component in bioremediation processes. In the realm of bioelectronics, Geobacter pili, also referred to as "nanowires, " exhibit electrical conductivity that is comparable to metallic wires, rendering them valuable in electronic applications. Additionally, G. sulfurreducens has emerged as a highly recommended species for MFCs. Positioned in proximity to an electrode, G. sulfurreducens engages in electrogenic metabolism, leveraging the electrode to anchor and release electrons that were generated during its metabolism, particularly in the presence of certain compounds like acetate ^{[18,19,20,21,22]}. The utilization of electrogenic bacteria in MFCs represents a compelling and innovative technology within the bioenergy field.

In addition, Geobacter demonstrates significant potential for application in the bioremediation of environments contaminated with organic pollutants and heavy metals, particularly under anaerobic conditions. The mechanism underlying bioremediation via Geobacter is characterized by its unique ability to reduce and degrade environmental pollutants through the direct transfer of electrons. This process leads to the formation of biofilms and the production of less toxic final residues ^[23,24]. Commonly targeted pollutants encompass organic contaminants such as hydrocarbons and solvents, as well as heavy metals like U(Ⅳ) and Fe(Ⅲ) ^[25,26].

The implementation of functional MFCs demands the intricate management of multiple components. These elements, working in tandem, enable the proper functioning of the cell, influencing aspects such as cell design (size, shape, material), electrode type (size, conductive material), choice of inoculum or biological organism (cultivation conditions for the electrogenic bacterium), biofilm formation and thickness, and the electron transfer mechanism (mediator-assisted or mediator-free). Each of these components plays a crucial role in bioelectricity generation. Although the essential components of a functional MFC have been identified, the interaction of microorganisms with the electrode and biofilm remains incompletely understood. Therefore, the development of appropriate mathematical models to simulate the microbial growth of Geobacter species is an actual task that could improve our understanding of complicated mechanisms underlying bacterial growth.

Thus, scant attention has been directed towards the development of a mathematical model for analysis of the growth of G. sulfurreducens, the formation of a biofilm, and the electricity generation in a MFC. Functioning in a similar manner to a biological battery, an MFC involves a complex interplay of various components and features, including the microorganism, electrode, cell type, substrates, membrane, and more. Consequently, the creation of a mathematical model to analyze the microbial growth of Geobacter is deemed valuable, offering the potential to optimize the intricate relationships among these components and features that dictate the overall operation of the cell. The reciprocal influence between mathematical and biological sciences has been pivotal in shaping the development of each discipline.

Mathematical models hold particular significance as they facilitate quantitative predictions, inference about underlying mechanisms, the testing of biological hypotheses, the design of experiments with living organisms, and quantitative depiction of relationships among diverse components and features within a biological entity ^[27]. Mathematics continues to assume an increasingly pivotal role in shaping the development of life, social, and biochemical sciences, with these sciences reciprocally influencing the evolution of mathematical methodologies. Mathematical modeling of biological systems has become exceptionally expansive, reflected in an extensive body of literature that is dedicated to the theoretical description of various biological phenomena. The surge in interest is attributed to advancements in diverse fields, facilitating the collection and management of substantial datasets, and thereby fostering the application of mathematical techniques to explore intricate biological systems ^[28]. Importantly, it is crucial to recognize that mathematical modeling cannot replace experiments with living organisms; rather, it serves as a complementary tool. Experiments remain indispensable for validating mathematical models and analyzing their predictions, thus highlighting the collaborative nature of mathematical modeling in tandem with other scientific disciplines.

The paper is structured as follows. Section 2 provides details on the experimental setup and presents the experimental results. In Section 3, the mathematical model is introduced, elucidating the biological mechanism of Geobacter proliferation and offering a comparison with experimental observations. Section 4 delves into the qualitative properties, encompassing the analysis of fixed points and saddle-node bifurcation, while Section 5 presents the main conclusions of the study.

2.   Experimental setup and results

2.1.   Conditions for activating the strain

The freeze-dried strain of G. sulfurreducens (ATCC 51573, Cientifica Senna SA de CV) was utilized in this experiment. To activate the strain, a sterile saline solution with a NaCl:H₂O ratio of 1:9 was prepared, and a 1/10 dilution was made. This sample was stored at 4C until needed. It is essential to note that G. sulfurreducens is classified as a Biosafety Level 1 (BSL-1) organism and has not been reported to cause disease in healthy adults ^[29].

2.2.   Preparation of the culture medium

We prepared 100 mL of a nonaqueous base fluid (NABF) medium ^[30] using the components listed in Table 1. The components in rows 1 to 5 were weighed and dissolved in deionized water. The pH was then adjusted to 7 with NaOH (1N); lasty, the remaining reagents listed in Table 1 (rows 6 to 10) were added, and the volume was adjusted with deionized water. All media were autoclaved at 121℃ and 15 lb of pressure for 15 minutes. They were kept refrigerated until use. For the experiment, the acetate concentration was varied, Medium 1 with 20 mM acetate and Medium 2 with 61 mM acetate. Note that Medium 2 comprised the NABF medium without modification and Sigma-Aldrich vitamin kit. We used the following compounds to fabricate the NB salt solution: KH₂PO$_4$, NH$_4$Cl, KCl, NaCl, and applied the following mineral elixir composition: NTA, FeSO$_4 \times 7$H₂O, CoCl$_2 \times 6$H₂O, ZnSO$_4 \times 7$H₂O, CuCl$_2 \times$ 2H₂O, AlK(SO$_4$)$\times 12$H₂O, H$_3$BO$_3$, NaH₂MoO$_4 \times 2$H₂O, NiSO$_4 \times 6$H₂O, NaH₂Wo$_4 \times 2$H₂O.

Table 1.  Composition of culture Medium 1 and Medium 2 for the proliferation of G. sulfurreducens. The X indicates that the components were contained in the medium. The difference between the two media is the concentration of sodium acetate, i.e., 20 mM for Medium 1 and 61 mM for Medium 2.

No	Components	Medium 1	Medium 2
$1$	Sodium acetate	20 mM	61 mM
$2$	Sodium carbonate heptahydrate	X	X
$3$	Calcium chloride monohydrate	X	X
$4$	Magnesium sulfate heptahydrate	X	X
$5$	Sodium bicarbonate	X	X
	pH 7 with NaOH 1N
$6$	Fumaric acid	X	X
$7$	NB Salt	X	X
$8$	Mineral elixir	X	X
$9$	Vitamin mix	X	X
$10$	Deionized water	X	X

---

 | Show Table

DownLoad: CSV

2.3.   Inoculation

The inoculum for this study consisted of G. sulfurreducens cells, prepared at a concentration of 10⁴ cells/mL. The cell concentration was accurately determined by using the CytoSMART cell counting chamber (Corning brand), combined with the Cell counter algorithm. 20 $\mu$L of each sample was taken and diluted with 20 $\mu$L of trypan blue before being loaded into the counting chamber where the dilution parameters were adjusted as needed. Subsequently, the captured images were subjected to analysis by using the cell counting algorithm provided by the cell counter software.

For the inoculation step, 2 mL glass tubes that had been previously filled with the culture media (Medium 1 and Medium 2) were utilized. The inoculation process was conducted within a semi-oxygen-free chamber, ensuring a constant COH₂ concentration of 5000 ppm. The entire test was performed in triplicate, and uninoculated media were included as controls to validate the experimental results.

2.4.   Bacterial growth

Measurements were conducted by using the Multiskan Go 1510-01413C automated device (Thermo Scientific, USA). A 96-well plate was loaded into the device containing Medium 1 and Medium 2, as well as the controls (non-inoculated media). The analysis was carried out over a period of three days at a temperature of 36℃ without shaking. Measurements were taken every hour, specifically at 580 nm.

2.5.   Experimental results

G. sulfurreducens was cultured in the NABF medium ^[30], which included components that serve as sources of carbon, energy, electrons, mineral salts, and vitamins. Among these components, acetate plays a crucial role. Acetate, an organic molecule, can be oxidized by bacteria to generate metabolic energy (adenosine triphosphate, ATP) and carbon structures for biomass production ^[31]. In our study, G. sulfurreducens was cultured in the NABF medium with different acetate concentrations: 20 mM and 61 mM.

The proliferation of G. sulfurreducens cells is illustrated in Figure 1. The cell densities reached the maxima of $OD_{580} = 0.14$ for the 20-mM acetate concentration and $OD_{580} = 0.15$ for the 61-mM acetate concentration.

Figure 1.  Experimental data depicting the G. sulfurreducens growth for acetate concentrations $A_{c} = 20$ mM (black dots) and $A_{c} = 61$ mM (dark yellow dots).

DownLoad: Full-Size Img PowerPoint

One can see that the Geobacter cell growth increases faster at $A_c = 61$ mM and saturates at higher values than that at $A_c = 20$ mM. It is important to note that some studies have reported bacterial growth at lower acetate concentrations of 5 mM and 5.5 mM, although these experiments often involved genetically modified bacteria ^[32] or increased acetate concentrations analyzed through protein quantification ^[33]. These various methodologies make direct comparison of results challenging. Another study utilized acetate concentrations of 10 mM and 15 mM, which are lower than the 20 mM base of the NABF medium. However, microbial growth curve analysis was not the primary focus of these experiments ^[34,35].

We should mention that our experimental protocols adhered strictly to established standards ^{[30,31,32,33]}. This adherence ensured that variations in experimental conditions were kept to a minimum, thereby enhancing the reliability and reproducibility of our results.

3.   Theoretical model

Every scientific study requires a theoretical model that simulates the relevant aspect of the scenario under examination. However, the reliability of such models is not always guaranteed, and they may sometimes produce inaccurate predictions. As a result, models may undergo modifications, be discarded, or continue to be used despite imperfections, solely because they provide a better understanding than having no model at all. The science of modeling involves aligning computational processes with series of measurements, consequently converting graphical and abstract representations into constructive sets of numbers. Our objective was to translate real-world processes into numerical form and compare them to the calculated values derived from potential theories.

3.1.   Model equations

Geobacter species in microbial culture media possess the ability to reduce various compounds, including metals and organic substances like acetate, and to utilize these reduced compounds as electron donors. The growth of Geobacter species involves intricate metabolic processes and adaptations that allow them to thrive in anaerobic environments, playing pivotal roles in soil biogeochemistry and bioremediation.

The growth of Geobacter species significantly surpasses that of the components of the microbial growth medium, rendering these components essentially fixed or changing very slowly. Consequently, the minimum number of variables includes the Geobacter population as the fast variable and the quantity of the microbial growth medium as the slow variable. Geobacter species derive their energy from the oxidation of organic compounds, such as acetate, or reduced metals, serving as electron donors. In a sequence of redox reactions, electrons are transferred to electron acceptors, which can be inorganic compounds like iron oxides or heavy metals, or even other microbes (acetates) through direct electron transfer processes.

Here, we introduce a straightforward mathematical model to simulate the experimental data obtained for the Geobacter species population $G$ evolution as depicted in Figure 1. Given the relatively slow rate of environmental change as compared to the rapid growth of the Geobacter species, we can treat the slowly varying environmental factors as effectively constant over short time intervals. Consequently, the slow variable can be adiabatically eliminated. In addition, in our model, we specifically emphasize the temporal dynamics, as our access was limited to experimental data that focuses solely on the temporal evolution of G. sulfurreducens.

The microbial growth in the absence of acetate oxidation can be expressed as a logistic function as follows ^[36,37]:

where $r_{G}$ represents the growth rate and $k_{G}$ denotes the carrying capacity. The carrying capacity is contingent on the proportions of various compounds, such as metals and organic substances like acetate. Geobacter species are renowned for their capacity for direct electron transfer, wherein bacteria can transfer electrons directly to other microbial species or reduced compounds without the need for chemical mediators. This mechanism fosters intimate and cooperative interactions within microbial communities.

If we consider the effects of the oxidation of organic compounds, such as acetate, in the microbial nutrient medium, the saturation term $f(G)$ should be added to the right side of equation (3.1). A decline in acetate oxidation prevents Geobacter species from deriving energy through the oxidation of organic compounds, leading to Geobacter saturation. Consequently, $f(G)$ should approach an upper bound $\gamma$ as $G$ approaches infinity. This limit $\gamma$ may vary depending on the microbial nutrient medium. At relatively low Geobacter species densities, acetate oxidation increases, providing energy to Geobacter species. This phenomenon can be modeled as an $S$-shaped function that resembles a saturation reaction. Therefore, a plausible form for $f(G)$ that allows saturation at a level $\gamma$ can be given as follows ^[36,37]:

where the parameter $\delta$ dictates the scale of Geobacter proliferation at which saturation occurs, reflecting a decrease in $\delta$ that is attributed to the acetate oxidation according to $A_c$. Subtracting the saturation function given by equation (3.2) from equation (3.1) results in the following differential equation that governs the rate of change of $G$:

The logistic and saturation functions respectively outlined in equations (3.1) and (3.2) elucidate fundamental properties of dynamical systems that are applicable to various biological models, including those concerning bacterial growth. Specifically, Ludwig et al. ^[36] employed these functions to characterize the dynamics of insect outbreak systems, focusing on the spruce budworm and balsam fir. This involved categorizing state variables into fast and slow components. Notably, the form of the saturation term in equation (3.2) corresponds to a type-Ⅲ S-shaped functional response, as described by Holling in 1959 ^[38]. These papers equipped with vital mathematical tools that can be leveraged to gain insights at three distinct levels: (ⅰ) qualitative and non-quantitative, (ⅱ) approximate parameter values, typically within the grasp of informed biologists when prompted by the appropriate inquiries, and (ⅲ) highly detailed quantitative data, which is infrequently available.

3.2.   Parameter optimization and comparison with the experimental results

Geobacter growth is influenced by the initial acetate concentration $A_c$ in the microbial nutrient medium. In Section 2 we considered two distinct initial concentrations of acetate, namely $A_{c} = 20$ mM and $A_{c} = 61$ mM. This difference in concentration can be considered in the model by varying in the parameters $r_G$, $k_G$, $\gamma$, and $\delta$ in equation (3.3) for each acetate concentration $A_c$.

To determine the values of the parameters in equation (3.3) that optimally describe Geobacter growth, we attempted to optimize the parameter values by minimizing the square root of the mean squared error between experimental and modeledGeobacter growth. In essence, our objective was to identify the error function, as follows ^[39]:

where $G_{exp-i}$ represents the individual experimental Geobacter growth data obtained in Section 2, while $G_{mod-i}$ denotes the Geobacter growth derived from the model equation (3.3). It is essential to emphasize that the time points used in the mathematical model were consistent with the times employed in the experiment. The numerical simulations of equation (3.3) were performed by using the classic Runge-Kutta method of 4th order MATLAB R2022's fmincon function facilitated the determination of the model parameters in equation (3.4).

The optimal values for the parameters $r_{G}$, $k_{G}$, $\gamma$, and $\delta$ in model equation (3.3) are detailed in Table 2. These values represent the best fit between the theoretical Geobacter cell density $G_{mod}$ and the experimental data denoted by $G_{exp}$. Figure 2 illustrates the dynamics of Geobacter cell density, with the red dots representing the experimental data $G_{exp}$ and the blue lines denoting the Geobacter cell density $G_{mod}$ derived from the solution of equation (3.3). Two distinct scenarios are depicted, corresponding to different acetate concentration values $A_c$ in the microbial nutrient medium: Figure 2(a) for $A_c = 20$ mM and Figure 2(b) for $A_c = 61$ mM. Furthermore, Figure 2(c) provides a comparison of the model data for the evolution of Geobacter growth under both acetate concentrations. The distinction is evident: while for $A_c = 20$ mM, the Geobacter cell density $G_{mod}$ exhibited rapid growth up to 45.5 hours, beyond point the Geobacter cell density began to saturate (black line), for $A_c = 61$ mM, $G_{mod}$ increased almost linearly (dark yellow line).

Table 2.  Optimized parameters for equation (3.3) to model experimental Geobacter cell density evolution ($G_{exp}$) for two acetate concentrations ($A_c$).

$A_{c}$	$r_{G}$	$k_{G}$	$\gamma$	$\delta$
20 Mm	0.7345	0.2209	0.0431	0.0389
61 Mm	0.0485	0.2225	0.0190	0.4831

---

 | Show Table

DownLoad: CSV

Figure 2.  Evolution of experimental $G_{exp}$ (red dots) and theoretical $G_{mod}$ (blue dots) (calculated by using equation (3.3)) for two acetate concentrations in the microbial nutrient medium: (a) $A_{c} = 20$ mM and (b) $A_{c} = 61$ mM. (c) Evolution of $G_{mod}$ as calculated by using equation (3.3) for $A_c = 20$ mM (black line) and $A_c = 61$ mM (dark yellow line).

DownLoad: Full-Size Img PowerPoint

The mechanism underlying the Geobacter species growth for the low acetate concentrations ($A_c = 20$ mM) can be understood as follows. Given that Geobacter species derive energy from the oxidation of organic compounds, the growth rate according to equation (3.3) reached significant value, specifically $r_{G} = 0.7345$. Consequently, this led to a decrease in the acetate oxidation according to equation (3.3) to the values $\delta = 0.0389$, resulting in Geobacter saturation ($\gamma = 0.0431$). Simultaneously, the carrying capacity $k_{G} = 0.2209$ in equation (3.3) remained relatively stable because the primary contribution to the increase in Geobacter growth was derived from oxidized acetates ($A_c$) rather than a combination of various compounds, including the metals and organic substances present in the microbial nutrient medium.

For the high acetate concentrations ($A_c = 61$ mM), Geobacter growth was driven by a distinct mechanism. In this scenario, the growth rate according to equation (3.3) assumed smaller values, specifically $r_{G} = 0.0485$. Consequently, Geobacter species derive their energy from the oxidation of acetate ($A_c$) at a slower rate, resulting in $\delta$ in equation (3.3) reaching higher values ($\delta = 0.4831$). This led to the Geobacter species not reaching saturation ($\gamma = 0.0190$), indicating sustained growth of Geobacter cell as denoted by $G_{mod}$.

4.   Model analysis

In this section, we present both the linear and nonlinear analyses of the model equation (3.3). This facilitates a deeper understanding of the behavior of Geobacter growth as well as more accurate predictions, under varying experimental conditions.

4.1.   Qualitative properties

To analyze the dynamics of equation (3.3), reducing the number of parameters is essential. The original equation (3.3) contains four parameters: $r_{G}$, $k_{G}$, $\gamma$, and $\delta$. Utilizing the fact that the parameter $\delta$ and the Geobacter cell density $G$ have the same dimension, equation (3.3) can be scaled through the transformation $g = G/\delta$. Dividing equation (3.3) by $\gamma$, we obtain:

Considering the following change $\tau = \dfrac{\gamma t}{\delta}$, $r = \dfrac{r_{G} \delta}{\gamma}$, and $k = \dfrac{k_{G}}{\delta}$, equation (4.1) in a dimensionless form becomes

where $r$ and $k$ denote the dimensionless growth rate and carrying capacity, respectively. The scaleless nature of equation (4.2) is advantageous because it involves only two parameters, facilitating the examination of its qualitative properties, which includes an analysis of fixed point stability and the determination of the bifurcation diagram.

4.2.   Fixed point analysis

Equation (4.2) possesses an unstable trivial fixed point (FP) at $g^{*} = 0$. This can be explained as follows. Geobacter rate of growth $g$ derives its energy from the oxidation of acetate $A_c$, which acts as an electron donor. When $A_c$ is large and $g$ is small, the Geobacter cells grow very rapidly when $g$ is close to zero. The remaining FPs of equation (4.2) satisfy the following equation:

In order to find the FP of equation (4.3), we rewrite equation (4.3) as

where the left-hand side represents a straight line with $k$ as the $x$-intercept and $r$ as the $y$-intercept. The right-hand side of this equation depicts a curve that originates from the origin and asymptotically approaches the $x$-axis for high values of $g$; notably, this side is independent of the parameters $r$ and $k$.

To determine the FP of equation (4.4), a geometric approximation method, as illustrated in Figure 3, was employed. This method involves plotting the graphs of $y_{1} = \dfrac{g}{1+g^{2}}$ and $y_{2} = r\left(1-\dfrac{g}{k}\right)$ on the same axes with respect to $g$ and observing their points of intersection between the results of $y_{3} = r\left(1-\dfrac{g}{k}\right)-\dfrac{g}{1+g^{2}}$ and the $g$-axis. The number and positions of these FPs depend on the parameters $k$ and $r$ from the straight line $y_{2} = r\left(1-\dfrac{g}{k}\right)$.

Figure 3.  FP $\left\{y_{1} = \dfrac{g}{1+g^{2}}, y_{2} = r\left(1-\dfrac{g}{k}\right), y_{3} = r\left(1-\dfrac{g}{k}\right)-\dfrac{g}{1+g^{2}}\right\}$ results for equation (4.4) as a function of $g$ for different values of $r$. (a, b) $r = 0.1619$. The intersection FP $g_{1}$ is stable. The flow directs to the right, where $y_{3} > 0$, as indicated by the arrows in (b). (c, d) $r = 0.2619$. Two intersection FPs, $g_{1}$ and $g_{S1}$, can be observed. $g_{1}$ remained stable, while $g_{S1}$ is a saddle-node point. The flow tends to the left, where $y_{3} < 0$, as shown in (d). (e, f) $r = 0.4372$. Three intersection FPs, $g_{1}$, $g_{2}$, and $g_{3}$, are apparent. $g_{2}$ is unstable, and the flow goes both to the left, where $y_{3} < 0$ and to the right, where $y_{3} > 0$, as illustrated in (f). In contrast, $g_{3}$ became stable, with the flow going to the right, where $y_{3} > 0$ and to the left, where $y_{3} < 0$. (g, h) $r = 0.5372$. Two intersection FPs, $g_{s2}$ and $g_{3}$, can be observed. $g_{S2}$ is a saddle-node point with the flow directed to the right, where $y_{3} > 0$, as shown in (h). $g_{3}$ is a stable fixed point. (i, j) $r = 0.6$. A stable intersection FP $g_{1}$ can be observed.

DownLoad: Full-Size Img PowerPoint

In Figure 3(a), where $k = 15$ is sufficiently large, we observe the existence of an intersection FP $g_{1}$ at $r = 0.1619$. This FP $g_{1}$ is stable, as indicated by the flow going to the right when $y_{3} > 0$ and to the left when $y_{3} < 0$, as shown by the arrows in Figure 3(b). Furthermore, when the value of $r$ was increased to 0.2619, as depicted in Figure 3(c), two intersection FPs $g_{1}$ and $g_{S1}$ become apparent. The stability analysis revealed that FP $g_{1}$ remained stable, while $g_{S1}$ became an unstable saddle-node FP. The flow directs towards the left where $y_{3} < 0$, as illustrated by the arrows in Figure 3(d). Additionally, the saddle node $g_{S1}$ resulted from the collision of two new FPs, i.e., an unstable FP $g_{2}$ and a stable FP $g_{3}$ when the parameter $r$ was decreased from values greater than 0.2619. As $r$ was further increased to $r = 0.4372$, three intersection FPs $g_{1}$, $g_{2}$, and $g_{3}$could be observed, as shown in Figure 3(e). The FP $g_{2}$ was determined to be unstable, as evidenced by the flow directing both to the left, where $y_{3} < 0$, and to the right, where $y_{3} > 0$, as illustrated in Figure 3(f). In contrast, the FP $g_{3}$ was determined to be stable, with the flow going to the right, where $y_{3} > 0$ and to the left, where $y_{3} < 0$, as depicted in Figure 3(f).

For a larger value of $r$ ($r = 0.5373$), two intersection FPs, $g_{S2}$ and $g_{3}$, became apparent, as shown in Figure 3(g). The stability analysis revealed that $g_{S2}$ was an unstable saddle node, with the flow directed to the right where $y_{3} > 0$, as indicated by the arrows in Figure 3(h). The emergence of this saddle node $g_{S2}$ resulted from the collision of a stable FP $g_{1}$ and an unstable FP $g_{2}$ as a result of the parameter $r$ increasing beyond 0.2619, while $g_{3}$ maintained a stable FP (Figure 3(h)). Finally, for $r = 0.6$, a stable intersection FP $g_{1}$ can be observed, as depicted in Figures 3(i) and 3(j).

4.3.   Bifurcation analysis

A saddle-node bifurcation (SNB) occurs where the curve $y_{1} = \dfrac{g}{1+g^{2}}$ tangentially intersects the line $y_{2} = r\left(1-\dfrac{g}{k}\right)$, i.e., the following conditions must be fulfilled: $y_{1} = y_{2}$ and $\dfrac{dy_{1}}{dt} = \dfrac{dy_{2}}{dt}$. The derivative of equation (4.4) gives:

Using equations (4.5) and (4.6), the parameters $k$ and $r$ can be obtained as functions of $g$, i.e.,

The condition $k > 0$ in equation (4.8) implies that $g$ must be confined to $g > 1$.

In Figure 4, we present the bifurcation diagrams of Geobacter growth with $r$ and $k$ serving as control parameters. These diagrams showcase two branches of stable FPs, denoted as $g_1$ and $g_3$, which can be seen to emerge and disappear at SNB points $g_{S1}$ and $g_{S3}$, respectively. Furthermore, there exists a branch of unstable solutions, represented by $g_2$, which connects these SNB points.

Figure 4.  Bifurcation diagrams for the Geobacter growth with respect to (a) $r$ for $k = 8$ and (b) $k$ for $r = 0.55$. The solid lines $g_1$ and $g_3$ are branches of stable solutions and blue dots indicate unstable FPs $g_2$. $g_{S1}$ and $g_{S3}$ are SNB points.

DownLoad: Full-Size Img PowerPoint

In Figure 5, we delineate the codimension-two bifurcation diagram across the $(r, k)$ parameter space, as outlined by equations (4.7) and (4.8). The diagram depicts the coexistence of stable ($g_1$ and $g_3$) and unstable ($g_2$) FPs, whose positions span the entire surface contingent upon the parameters $r$ and $k$. The two SNB curves, labeled as $g_{S1}$ and $g_{S2}$, delineate the stability thresholds of these FPs. Along these curves, all three FPs coexist. Conversely, beyond these boundaries, only one stable FP exists, with $g_1$ occurring to the left of this FP and $g_3$ to its right. The intersection point of the SNB curves denotes the cusp point (CP), located at $(k, r) = (5.16, 0.65)$.

Figure 5.  Codimension-two bifurcation diagram for the results of equations (4.7) and (4.8) in the ($r, k$)-parameter space. The blue lines represent the SNB curves, with CP denoting the cusp point. Dashed lines indicate the values of the fixed parameters used to construct the bifurcation diagrams in Figure 4.

DownLoad: Full-Size Img PowerPoint

The CP occurring at the specific values of $k$ and $r$ was found to possess a unique property: changes in the parameters ($k, r$) lead to an abrupt transition of the results of solving system equation (4.1). This transition occurs between the FPs $g_{1}$ and $g_{3}$, bypassing the critical region delineated by the bifurcation curves $g_{S1}$ and $g_{S2}$. This abrupt shift can be likened to a discontinuous drop or a catastrophic event.

By utilizing the parameters $r_{G}$, $k_{G}$, $\gamma$, and $\delta$ from equation (3.3) for various acetate concentrations ($A_c$), as presented in Table 2, we can derive the parameters $r = \dfrac{r_{G} \delta}{\gamma}$ and $k = \dfrac{k_{G}}{\delta}$ for equation (4.2). For $A_{c} = 20$ mM, the corresponding dimensionless parameters were determined to be $k = 5.6787$ and $r = 0.6629$, resulting in a stable FP $g_3$. Conversely, for $A_{c} = 61$ mM, with $k = 0.4606$ and $r = 1.2332$, the solution lies below the CP and the two stable equilibria ($g_1$ and $g_3$) become indistinguishable, i.e., they collapse into an FP. Consequently, Geobacter growth is not saturated.

This observation can be described by considering the growth rate $r_{G}$ in equation (3.3) for $A_{c} = 61$ mM, as indicated by the dark yellow curve in Figure 2(c). At $r_{G} = 0.0485$, the Geobacter species derive energy from acetate oxidation at an extremely slow rate. Consequently, $\delta$ in equation (3.3) reaches high values ($\delta = 0.4831$), implying that the Geobacter species do not reach saturation ($\gamma = 0.0190$). In other words, the Geobacter cell density $G_{mod}$ continues to grow after the slow growth phase.

5.   Discussion

Mathematical models for bacterial growth span a spectrum of complexity, ranging from simple one-term exponentials to highly nonlinear functions ^{[40,41,42,43,44,45,46,47,48,49,50]}. A robust model should employ the minimum necessary parameters to produce predicted values that are consistent with experimental data, with these parameters directly linked to significant biological properties of the system under examination.

5.1.   Overview of bacterial growth models

Some researchers have employed generic growth models that comprise merely four parameters, and they were found to be capable of effectively elucidating numerous biological systems ^[40]. These models are formulated by imposing suitable limits or conditions on selected parameters within the constraints of the generic curve. For instance, under conditions in which population growth is consistently proportional to its size without limiting factors, the generic curve simplifies to a simple exponential curve. Introducing a single limiting factor results in the logistic curve. Specialized cases of the generic growth model have been employed in the analysis of diverse experimental data, such as bacterial agglutination (hyper-Gompertz), complement-mediated hemolysis (hyperlogistic), and transplantable tumor growth (logistic) ^[41].

In exponential models, each subset of cells diminishes in proportion to its size at any given moment, without any discernible limiting factor or initial delay. However, if there is a delay in this diminishment, survival curves become more intricate, making the simple exponential model unsuitable for accurately fitting the data. In such scenarios, logistic models prove invaluable for analyzing these complex survival curves, accommodating both homogeneous and heterogeneous cell populations. Particularly, when examining a homogeneous cell population with growth rates that are related to the cell concentration and the concentration of a single growth-limiting factor, the logistic model emerges as a fitting choice for analysis.

Logistic models have been extensively used to describe bacterial growth ^[42]. For example, Fujikawa et al. ^[43] developed a logistic model that effectively described the growth of Escherichia coli and Salmonella, with subsequent improvements ^[44] enhancing predictive accuracy of the growth rate and lag time. Tashiro and Yoshimura ^[45] proposed a neo-logistic model that incorporates inducible enzyme synthesis, providing clearer parameter interpretation. Pinto and Shimakawa ^[46] introduced a compressed logistic model for bacterial growth, integrating a time-dependent rate, which they contend offers a more robust physiological framework. They assert that this model adeptly replicates experimental findings across diverse bacterial species. Additionally, they investigated the potential fractal-like characteristics of growth rates and proposed the model as a Gaussian representation of bacterial temporal evolution. Recently, Lo Grasso and colleagues ^[47] introduced a generalized logistic model, demonstrating strong performance in terms of fitting experimental data for Pseudomonas fluorescens. These models offer diverse enhancements and applications for predicting bacterial growth.

In our model, alongside bacterial logistic growth, we account for the oxidative impact of organic compounds, such as acetate, within the microbial nutrient medium, as reflected by the saturation term $f(G)$. Specifically, all strains of Geobacter possess a Gram-negative rod-shaped morphology and the ability to oxidize acetate while concurrently reducing Fe(Ⅲ). These bacteria commonly showcase additional capabilities, including the reduction of Mn(Ⅳ), U(Ⅵ), elemental sulfur, and humic substances. Notably, G. sulfurreducens primarily utilizes acetate (Ac) as its primary electron donor, deriving energy from the oxidation of organic compounds like acetate or reduced electron-donor metals. The saturation effect resulting from the oxidation of organic compounds within the microbial nutrient medium, as demonstrated for acetate, is encapsulated within our mathematical model through the term $f(G)$.

It is noteworthy that several researchers have conducted bifurcation and stability analyses of bacterial growth. For instance, Aviram and Rabinovitch ^[51] discovered the coexistence of bacteria and bacteriophage in the presence of bacterial debris, revealing various periodic behaviors and chaotic attractors. Subsequently, Giverso et al. ^[52] investigated pattern formation in the expansion of bacterial colonies, uncovering linear and branching instabilities. Furthermore, Ren and Yuan ^[53] delved into the dynamic behavior of a microbial continuous culture model, focusing on fold and Hopf bifurcations in the unforced system, as well as periodic solutions in the forced system. In a recent study, Ma et al. ^[54] studied the structure and stability of steady states within a bacterial colony model featuring cell density-suppressed motility. They established the existence and structure of positive solutions, along with criteria for determining the stability or instability of bifurcation branches. While these studies collectively contribute to our understanding of the complex dynamics of bacterial growth, the fixed-point stability and bifurcation analyses conducted in the present study further deepen our comprehension of bacterial growth dynamics.

5.1.1.   Comparison with a linear model

It is important to emphasize that our model offers a significant advantage over existing methodologies that are commonly utilized among biotechnological and bioprocess researchers studying bacterial kinetics. For instance, the linear model of bacterial growth introduced by Zwietering et al. ^[55], which is widely employed in biotechnological and bioprocess applications, is expressed as the following equation:

where $G$ represents Geobacter cell density and $r_{G}$ is the growth rate.

The optimal parameters describing the experimental evolution of Geobacter cell density $G_{exp}$ according to equation (5.1) were determined to be $r_{G} = 0.0227$ for $A_c = 20$ mM and $r_{G} = 0.0237$ for $A_c = 61$ mM. Illustrated in Figure 6 are the experimental $G_{exp}$ and the theoretically obtained $G_{mod}$, i.e., according to equation (5.1) for two acetate concentrations in the microbial nutrient medium: $A_c = 20$ mM (Figure 6(a)) and $A_c = 61$ mM (Figure 6(b)). It is evident that there is a significant lack of correlation between the experimental data (represented by the red dots) and the numerical solution (indicated by the blue lines).

Figure 6.  Experimental $G_{exp}$ (red dots) and theoretical $G_{mod}$ (blue lines) (calculated by using equation (5.1)) for two acetate concentrations in the microbial nutrient medium: (a) $A_c = 20$ mM and (b) $A_c = 61$ mM. (c) Evolution of $G_{mod}$ according to equation (5.1) for $A_c = 2$ (black line) and $A_c = 5$ (dark yellow line).

DownLoad: Full-Size Img PowerPoint

Figure 6(c) presents the trajectory of Geobacter cell density $G_{mod}$ according to equation (5.1) for $A_c = 20$ mM (depicted by the black line) and $A_c = 61$ mM (depicted by the dark yellow line). The differences between the curves were minimal, as the values of the growth rate $r_{G}$ were very similar. Notably, in Figure 6(c), the Geobacter cell density $G_{mod}$ displays continuous growth, which is biologically implausible, as it suggests an unlimited food source in the microbial nutrient medium. In reality, the availability of compounds such as metals and organic compounds (acetate) regulates this growth. Furthermore, equation (5.1) fails to elucidate how geobacteria derive energy from acetate oxidation, experience saturation, and react to the diverse compounds present in the microbial nutrient medium. In contrast, our mathematical model given by equation (3.3) can account for these processes through the parameters $r_{G}$, $k_{G}$, $\gamma$, and $\delta$ in the equation.

5.1.2.   Comparison with a logistic model without saturation

Since the traditional logistic model does not take into account the saturation effect, it has only two parameters ($r_G$ and $k_G$), as indicated by equation (3.1). The results of the parameter optimization of this model with the experimental data using the determination coefficient $R^2$ are presented in Figure 7.

Figure 7.  Experimental $G_{exp}$ (red dots) and theoretical $G_{mod}$ (blue lines) (calculated by equation (3.1)) for two acetate concentrations in the microbial nutrient medium: (a) $A_c = 20$ mM and (b) $A_c = 61$ mM.

DownLoad: Full-Size Img PowerPoint

Figure 7 shows the trajectories of Geobacter cell density $G_{mod}$ according to equation (3.1) and that obtained experimentally for $A_c = 20$ mM (Figure 7(a)) and $A_c = 61$ mM (Figure 7(b)). Although the theoretical curves indicate a tendency to reach saturation, the saturation was not consistent with that shown in the experimental curve, as determined by our four parametric model ($r_{G}$, $k_{G}$, $\gamma$, and $\delta$) given by equation (3.3) and displayed in Figure 2(c).

The results of the comparative analysis of linear and logistic models are presented in Table 3. The parameter values were adjusted by calculating the root mean square error, i.e.,

Table 3.  Optimal parameters, RMSE, and Pearson correlation coefficient R for the linear and logistic models under acetate concentrations $A_c = 20$ mM and $A_c = 61$ mM.

Parameters	Equation	RMSE	R
$A_c=20$ mM
Linear model:	(5.1)	0.198	0.8920
$r_G=0.0227$
Logistic model:	(3.1)	0.0049	0.9921
$r_G=0.0594, \ k_{G}=0.1436$
Our model:	(3.3)	0.0043	0.9950
$r_G=0.7345, \ k_{G}=0.2209$
$\gamma=0.0431, \ \delta=0.0389$
$A_c = 61$ mM
Linear model:	(5.1)	0.0152	0.9411
$r_G=0.0237$
Logistic model:	(3.1)	0.0055	0.9893
$r_G=0.0481, \ k_{G}=0.1650$
Our model:	(3.3)	0.0040	0.9942
$r_G=0.0485, \ k_{G}=0.2225$
$\gamma=0.0190, \ \delta=0.4831$

---

 | Show Table

DownLoad: CSV

and the Pearson correlation coefficient, i.e.,

between the experimental cell density $G_{exp}$ and the Geobacter cell density $G_{mod}$ derived from mathematical models of bacterial growth, i.e., a linear model, a logistic model without saturation (with two parameters), and our logistic model with a saturation term $f(G)$ (with four parameters), for two different values of acetate concentration ($A_c = 20$ mM and $A_c = 61$ mM).

6.   Conclusion

We have successfully developed a mathematical model to analyze the microbial growth of Geobacter cell density, providing a comprehensive framework to optimize the relationships among the various variables that dictate cellular function. The model for the growth of the Geobacter population consists of two segments. The first segment incorporates a logistic form that includes the growth rate ($r_{G}$) and carrying capacity ($k_{G}$). The carrying capacity is dependent upon the presence of diverse compounds, such as metals and organic compounds like acetate, and adheres to the condition of acetate oxidation absence. The second segment introduces the effect of organic compound oxidation (acetate) in the microbial nutrient medium as the saturation term $f(G)$. This segment assumes that the saturation level ($\gamma$) is dependent on the microbial nutrient medium, and the parameter $\delta$, which establishes the scale of Geobacter cell densities at which saturation commences. The decrease in $\delta$ is attributed to acetate oxidation $A_c$.

The mathematical model developed in this study effectively captures the microbial growth of Geobacter cell density, providing a comprehensive framework for optimizing the relationships among various variables that influence cellular function. By deriving the values of parameters $r_{G}$, $k_{G}$, $\gamma$, and $\delta$ from experimental data, we found that our model allows us to uncover the mechanisms hat regulate the evolution of Geobacter cell density. The parameters that best describe the experimental growth of Geobacter species were determined by minimizing the square root of the mean square error between the experimental data and the mathematical model equation (3.3). Our results highlight the significance of considering nonlinearity and saturation.For various parameter values, we were able to elucidate the mechanisms governing Geobacter cell density evolution, with crucial aspects corroborated by experimental data.

To analyze the system dynamics, we scaled equation (3.3) to a dimensionless form, represented by equation (4.2). This scaling allowed us to reduce the number of parameters, originally $r_{G}$, $k_{G}$, $\gamma$, and $\delta$, to two dimensionless parameters:, i.e., $r$ and $k$, signifying the dimensionless growth rate and carrying capacity, respectively. The significance of the scaleless equation (4.2) lies in its simplicity, as it is characterized by two parameters, which facilitates the exploration of its qualitative properties. This includes the analysis of fixed point stability and the determination of the codimension-two bifurcation diagram.

Finally, we conducted a comparative analysis between our model and a foundational linear model ^[55], which has been extensively utilized in biotechnological and bioprocess contexts to explore bacterial kinetics. Our findings revealed clear advantages of our model over the linear model in terms of predicting bacterial growth dynamics.

To conclude, we anticipate that our proposed model will prove to be a valuable tool for microbiologists engaged in simulating bacterial growth experiments. By capturing the intricate dynamics of Geobacter population growth, our model, with its nuanced consideration of nonlinear terms and saturation effects, offers a comprehensive framework for understanding the behavior of these microorganisms under varying conditions. Moreover, the advantages of our model, including its scalability to a dimensionless form and the resulting ease of parameter interpretation, contribute to its utility in diverse experimental scenarios. The insights gained from the model's ability to elucidate mechanisms of Geobacter cell density evolution, as supported by experimental data, enhance its applicability in real-world settings. We envisage that our model will not only serve as a theoretical foundation for studying microbial growth but that it will also find practical applications in the optimization of bacterial cultures for specific outcomes. As our understanding of microbial processes continues to evolve, the proposed model will stand as a versatile and insightful tool for researchers in the field of microbiology.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence tools in the creation of this article.

Conflict of interest

Alexander N. Pisarchik is an special issue editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.