CompLaB v1.0: a scalable pore-scale model for flow, biogeochemistry, microbial metabolism, and biofilm dynamics

Microbial activity and chemical reactions in porous media depend on the local conditions at the pore scale and can involve complex feedback with fluid flow and mass transport. We present a modeling framework that quantitatively accounts for the interactions between the bio(geo)chemical and physical processes and that can integrate genome-scale microbial metabolic information into a dynamically changing, spatially explicit representation of environmental conditions. The model couples a lattice Boltzmann implementation of Navier–Stokes (flow) and advection–diffusion-reaction (mass conservation) equations. Reaction formulations can include both kinetic rate expressions and flux balance analysis, thereby integrating reactive transport modeling and systems biology. We also show that the use of surrogate models such as neural network representations of in silico cell models can speed up computations significantly, facilitating applications to complex environmental systems. Parallelization enables simulations that resolve heterogeneity at multiple scales, and a cellular automaton module provides additional capabilities to simulate biofilm dynamics. The code thus constitutes a platform suitable for a range of environmental, engineering and – potentially – medical applications, in particular ones that involve the simulation of microbial dynamics.

1 Introduction

Biogeochemical turnover in Earth's near-surface environments is governed by the activity of microbes adapted to their surroundings to catalyze reactions and gain energy. In turn, these activities shape the environmental composition, which feeds back on metabolic activities and creates ecological niches. Such feedbacks can be captured by reactive transport models that compute the evolution of geochemical conditions as a function of time and space and simulate microbial activities in porous media (Meile and Scheibe, 2019). Commonly used macroscopic reactive transport models simplify small-scale features of natural porous media. For example, heterogeneous pore geometry and transport phenomena are represented by only a few macroscopic parameters such as porosity, permeability, and dispersivity (Steefel et al., 2015). However, such simplifications can lead to a disparity between model estimations and actual observations because these models do not resolve the physical and geochemical conditions at the scale that is relevant for microbial activity (e.g., Molins, 2015; Oostrom et al., 2016). Furthermore, microbial reaction rates are often formulated using Monod expressions, which describe a dependency of metabolic rates on nutrient availability but substantially simplify the complex metabolic adaptation of microbes in changing environments. This recognition has prompted the development of constraint-based models including, for example, COMETS (Harcombe et al., 2014), BacArena (Bauer et al., 2017), and IndiMeSH (Borer et al., 2019), which have enabled detailed descriptions of complex microbial metabolisms and metabolic interactions (Dukovski et al., 2021). However, most constraint-based models are not designed to capture combined diffusive and advective transport of metabolites in heterogeneous subsurface environments and are not optimized to handle such settings in a computationally efficient way. Notably, computational efficiency and the integration of adequate formulations of microbial function have been identified as critical aspects in pore-scale models of microbial activity (Golparvar et al., 2021).

To account for the feedback between environmental conditions, chemical processes, microbial metabolism, and structural changes in the porous medium caused by these activities, we introduce a novel pore-scale reactive transport modeling framework with spatially explicit descriptions of hydrological and biogeochemical processes. Our work complements existing efforts, encompassing both individual- and population-based spatially explicit microbial models reviewed by König et al. (2020), some of which take into consideration the structure of the porous medium. Our modular framework is developed to account for various chemical reactions and/or genome-scale metabolic models with advective and diffusive transports in porous media at the pore scale. The lattice Boltzmann (LB) method is used to compute fluid flow and solute transport in complex porous media, capable of simulating both advection- and diffusion-dominated settings. Microbial metabolism and chemical reactions are incorporated as source or sink terms in the LB method solving mass conservation equations. These sources or sinks can be described classically using approximations such as Monod kinetics (Tang et al., 2013) or can be derived from cell-scale growth and metabolic fluxes simulated with flux balance analysis (Orth et al., 2010). Biomass dynamics can be described by keeping track of cell densities (similar to chemical concentration fields) of different organisms or populations, with cell movement either based on an advection–diffusion formulation or using a cellular automaton approach. In addition, we incorporate a surrogate modeling approach to make larger-scale simulations possible. Thus, the framework provides options either to maximize computational efficiency via the use of surrogate models or to directly utilize well-established metabolic modeling environments without losing the inherent parallel scalability of the LB method. The model is validated by comparing model simulations to published simulation results. We demonstrate the flexibility of the new microbial reactive transport framework, its scalability, and the benefits of using surrogate models to circumvent computational bottlenecks posed by flux balance analysis. Our work therefore facilitates cross-disciplinary efforts that integrate bioinformatic approaches underlying cell models with descriptions suitable to resolving the dynamic nature of natural environments. This allows for the representation of microbial interactions, which is a major challenge to our current quantitative understanding of microbially mediated elemental cycling (Sudhakar et al., 2021).

2 Use of open-source codes

To establish a modeling framework that builds on the existing and future knowledge and know-how from multiple disciplines, our approach uses the open-source software Palabos (Parallel Lattice Boltzmann Solver) and integrates it with the open-source linear programming solver GLPK (GNU Linear Programming Kit) and the COBRApy (CPY) Python package for genome-scale metabolic modeling. Palabos is a modeling platform that has established itself as a powerful approach in the field of computational fluid dynamics based on the lattice Boltzmann (LB) method. The Palabos software is designed to be highly extensible to couple complex physics and other advanced algorithms without losing its inherent capability of massive parallelization (Latt et al., 2021). Palabos has been parallelized using the Message Passing Interface, where computational domains are subdivided while minimizing the inter-process communication. It has been used for building modeling platforms to simulate deformable cell suspensions in relation to blood flows (Kotsalos et al., 2019) and complex subsurface biogeochemical processes at the pore scale (Jung and Meile, 2019, 2021). It is highly scalable and hence was chosen as a high-performance modeling framework to be integrated with our representations of chemical and microbial dynamics. GLPK is an open-source library designed for solving linear programming (LP), mixed integer programming and other related problems (GLPK, 2022). It contains the simplex method, a well-known efficient numerical approach to solve LP problems, and the interior-point method, which solves large-scale LP problems faster than the simplex method. GLPK provides an application programming interface (API) written in C language to interact with a client program. COBRApy is an object-oriented Python implementation of constraints-based reconstruction and analysis (COBRA) methods (Ebrahim et al., 2013), which is suitable to be integrated with other libraries without requiring commercial software. Through a simple Python API, the fast evolving and expanding biological modeling capacity of COBRApy, which includes features such as flux balance analysis (FBA), flux variability analysis (FVA), metabolic models (M models), and metabolism and expression models (ME models), can be employed.

Figure 1 Flowchart of CompLaB for reactive transport simulations. Changes in pore geometry are assessed due to biomass changes. The state variables involved in the reactions can represent dissolved chemicals and planktonic microbes and solid phases or sessile microorganisms.

3 Model description

CompLaB simulates a fully saturated 2D fluid flow and solute transport at the pore scale based on the LB method implemented in Jung and Meile (2019, 2021). These earlier efforts established some of the underlying model developments, such as the simulation of the flow field, mass transport, and biochemical processes including kinetic rate expressions and cellular automaton implementation of biofilm growth. This study expands on the previously established models to offer a much broader applicability by building the modular structure that makes the use of flux balance and surrogate models possible. The LB method is particularly useful for simulating subsurface processes because boundaries between solid and fluid can be handled by a simple bounce-back algorithm (Ziegler, 1993) in addition to its massive parallelization efficiency (Latt et al., 2021). For these reasons, the LB method has been applied to simulate a broad range of pore-scale reactive transport processes (e.g., Huber et al., 2014; Kang et al., 2007; Tang et al., 2013). The simulations are guided by an input file, CompLaB.xml, that sets the scope of the simulation and capabilities utilized through command blocks that define the model domain, chemical state variables, microorganisms, and model input and output (Fig. 1). Below, the features of the model are presented, and we refer to the manual in the code repository for examples of the implementation.

3.1 The lattice Boltzmann flow and mass transport solvers

The LB method retrieves the numerical solutions of the Navier–Stokes (NS) for fluid flow and advection–diffusion-reaction equations (ADREs) for solute transport by solving the mesoscopic Boltzmann equation across a defined set of particles (Krüger et al., 2017). CompLaB obtains a steady-state flow field by running a flow solver with a D2Q9 lattice BGK (Bhatnagar–Gross–Krook; Bhatnagar et al., 1954) model defined as

$$\begin{array}{}\text{(1)}& f_i(\mathit{r}+{\mathit{c}}_i\mathrm{\Delta }t,t+\mathrm{\Delta }t)=f_i(\mathit{r},t)-{\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\tau }}_{\mathrm{f}}}}\left[f_i-f_i^{\text{eq}}\right],\end{array}$$

where f_i(r,t) is the ith discrete set of particles streamed from a position r to a new position r+c_i after a time step with lattice velocities c_i(c₀ = [0,0], c₁ = [1,0], c₂ = [0,1], c₃ =  [−1,0], c₄ = [0,−1], c₅ = [1,1], c₆ = [−1,1], c₇ = [−1,−1], c₈ = [1,−1]). τ_f is the relaxation time related to the fluid viscosity ν_f (${\mathit{\nu }}_{\mathrm{f}}=c_{\mathrm{s}}^{\mathrm{2}}({\mathit{\tau }}_{\mathrm{f}}-\frac{\mathrm{\Delta }t}{\mathrm{2}})$; c_s is a lattice-dependent constant; here, $c_{\mathrm{s}}^{\mathrm{2}}$ = $\mathrm{1}/\mathrm{3}$). The equilibrium distribution function for fluid flow ($f_i^{\text{eq}}$) is given by

$$\begin{array}{}\text{(2)}& f_i^{\text{eq}}={\mathit{\omega }}_i\mathit{\rho }+{\mathit{\omega }}_i{\mathit{\rho }}_{\mathrm{0}}\left({\displaystyle \frac{\mathit{u}\cdot {\mathit{c}}_i}{c_{\mathrm{s}}^{\mathrm{2}}}}+{\displaystyle \frac{{\left(\mathit{u}\cdot {\mathit{c}}_i\right)}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{4}}}}-{\displaystyle \frac{\mathit{u}\cdot \mathit{u}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{2}}}}\right),\end{array}$$

where ω_i is the lattice weights (ω₀ = $\mathrm{4}/\mathrm{9}$, ω₁₋₄ = $\mathrm{1}/\mathrm{9}$, ω₅₋₈ = $\mathrm{1}/\mathrm{36}$), ρ is the macroscopic density (ρ=Σf_i), ρ₀ is the rest state constant, and u is the macroscopic velocity calculated from the momentum (ρu = Σc_if_i). The steady-state flow field is then imposed on a transport solver defined as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{g}_i^j(\mathit{r}+{\mathit{c}}_i\mathrm{\Delta }t,t+\mathrm{\Delta }t)=& g_i^j(\mathit{r},t)-{\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\tau }}_{\mathrm{g}}^j}}\left[g_i^j-g_i^{j,\text{eq}}\right]\\ & +{\mathrm{\Omega }}_i^{\text{RXN}}(\mathit{r},t),\end{array}\end{array}$$

where g_i(r,t) represents the discrete particle set i of a transported entity j at position r and time t. ${\mathit{\tau }}_{\mathrm{g}}^j$ is the relaxation time for solute transport related to the diffusivity (D^j = $c_{\mathrm{s}}^{\mathrm{2}}({\mathit{\tau }}_{\mathrm{g}}^j-\frac{\mathrm{\Delta }t}{\mathrm{2}})$). The equilibrium distribution function for transport ($g_i^{j,\mathrm{eq}}$) using a D2Q5 lattice for numerical efficiency, which, satisfying the isotropy requirement for a LB transport solver, is given by

$$\begin{array}{}\text{(4)}& g_i^{j,\text{eq}}={\mathit{\omega }}_i{C}^j\left(\mathrm{1}+{\displaystyle \frac{\mathit{u}\cdot {\mathit{c}}_i}{c_{\mathrm{s}}^{\mathrm{2}}}}\right),\end{array}$$

with the lattice weights ω₀ = $\mathrm{1}/\mathrm{3}$ and ω₁₋₄ = $\mathrm{1}/\mathrm{6}$, lattice velocities c₀₋₄, and the solute concentration C^j ($C^j=\sum g_i^j$). In solving an advection–diffusion problem, CompLaB adjusts the value of ${\mathit{\tau }}_{\mathrm{g}}^j$, which controls the length of a time step, to obtain a user-provided Péclet number (Pe^j = $UL/D^j$) for a given average flow velocity U and a user-provided characteristic length L. With a steady-state flow field obtained from the solution of Eq. (1) and a reaction step ${\mathrm{\Omega }}_i^{\text{RXN}}=\mathrm{\Delta }t{\mathit{\omega }}_{i}R$, the LB transport solver recovers an ADRE with the following form:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial j}{\partial t}}=\mathrm{\nabla }\cdot \left(D^j\mathrm{\nabla }j\right)-\mathit{u}\cdot \mathrm{\nabla }j+R.\end{array}$$

Here, the transported entity j includes solute concentrations (C) and planktonic biomass densities, and R is a reaction term computed by the reaction solver of CompLaB as described below.

3.2 Reactions

The reaction step (${\mathrm{\Omega }}_i^{\text{RXN}}$) computation is separated from a transport computation via the sequential non-iterative approach (Alemani et al., 2005). A unique feature of CompLaB is that its reaction solver can compute biochemical reaction rates R through (1) kinetic rate expressions, (2) flux balance analysis, and (3) a surrogate model such as a pre-trained artificial neural network or combinations thereof, by summing their contributions to the net reaction rates of individual state variables.

3.2.1 Kinetic rate expressions

CompLaB provides a C$++$ template that users can adapt to formulate kinetic rate expressions using metabolite concentrations and biomass densities (defineReactions.hh). This is designed to accommodate user-specific needs and to enable simulating various microbial dynamics including Monod kinetics, microbial attachment and detachment, and arbitrary rate expressions defined by the user. Reactions can be restricted to particular locations using material numbers (mask) differentiating fluid, biomass, and grain surfaces.

Local biomass densities and concentrations calculated after the collision step of transport solvers are transferred to the function as vectors B and C, where the vector elements follow the order defined in the user interface (CompLaB.xml). The biomass density and metabolite concentrations are updated according to

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{B}}_{t+\mathrm{\Delta }t}={\mathit{B}}_{\mathrm{t}}+{\mathit{\gamma }}_{\mathrm{t}}{\mathit{B}}_{\mathrm{t}}\mathrm{\Delta }t,\text{(7)}& {\displaystyle }{\mathit{C}}_{t+\mathrm{\Delta }t}={\mathit{C}}_{\mathrm{t}}+{\mathit{R}}_{\mathrm{t}}\mathrm{\Delta }t,\end{array}$$

where γ_t values are the cell-specific biomass growth rates, R_t values are the microbially mediated reaction rates, expressed as the product of metabolite uptake/release rates per cell (F_t) and the cell density B_t (i.e., R_t=F_tB_t), calculated every time step for every pore and biomass grid cell.

3.2.2 Flux balance analysis

For genome-enabled metabolic modeling, CompLaB loads metabolic networks and calculates microbial growth rates as well as metabolite uptake/release rates through an FBA method (Orth et al., 2010). FBA investigates the metabolic capabilities by imposing several constraints on the metabolic flux distributions. Assuming that metabolic systems are at steady state, the system dynamics for a metabolic network are described by the mass balance equation Sv = 0. Here, S is a m × n matrix with m compounds and n reactions, where the entries in each column are the stoichiometric coefficients of the metabolites composing a reaction and v is a n × 1 flux vector representing metabolic reactions and uptake/release of chemicals by the cell. Most metabolic models have more reactions than compounds (n > m), meaning that there are more unknowns than equations. To solve such underdetermined systems, FBA confines the solutions to a feasible set by imposing constraints on metabolic fluxes lb (lower bounds) ≤ v ≤ ub (upper bounds) and applies an objective function f(v) = c^Tv, where c is the vector of weights for the objective function to identify an optimal solution. Commonly used objective functions include maximization of biomass yield, maximization of ATP production, and minimization of nutrient uptake (Nikdel et al., 2018). CompLaB utilizes the stoichiometric matrix S from standard metabolic databases such as BiGG and KEGG, which are widely used in FBA simulation environments (e.g., COBRA toolbox, Heirendt et al., 2019; COBRApy, Ebrahim et al., 2013; KBase, Arkin et al., 2018). Therefore, CompLaB can integrate many existing in silico cell models.

CompLaB computes the solution of the metabolic models at each point in space and time for each organism or microbial community (if the model represents multiple microorganisms) and updates biomass density and metabolite concentrations according to Eqs. (6) and (7). The metabolic uptake fluxes are set through imposing constraints by defining the lower bound (lb) of a chemical (uptake fluxes are negative) through one of the following approaches. The first is the parameter-based method employed by Harcombe et al. (2014), setting the metabolic fluxes in analogy with Michaelis–Menten kinetics using a maximum uptake rate (V_max; e.g., $\mathrm{mmol}{\mathrm{g}}_{\text{DW}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$, where g_DW is gram dry weight):

$$\begin{array}{}\text{(8)}& \text{lb}=-V_{max}\left({\displaystyle \frac{C}{C+K_{\mathrm{s}}}}\right),\end{array}$$

where C is a local metabolite concentration (e.g., mM) and K_s is a half-saturation constant (e.g., mM). The second is the semi-linear approach employed by Borer et al. (2019). This method replaces V_max with C(BΔt)⁻¹, where B is a local biomass density (e.g., g_DW L⁻¹) and Δt is the length of a time step (e.g., h):

$$\begin{array}{}\text{(9)}& \text{lb}=-{\displaystyle \frac{C}{B\mathrm{\Delta }t}}\left({\displaystyle \frac{C}{C+K_{\mathrm{s}}}}\right).\end{array}$$

If K_s is set to 0, then the uptake flux estimation becomes a linear function to local concentrations. Note that the units in the fluid flow and mass conservations model simulation must match those of the FBA bounds, which in our case were $\mathrm{mmol}{\mathrm{g}}_{\text{DW}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$. With lower bounds defined, the solution of an FBA problem outputs biomass growth rate (γ; e.g., h⁻¹) and uptake/release rates of metabolites (F; e.g., $\mathrm{mmol}{\mathrm{g}}_{\text{DW}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$).

3.2.3 Surrogate model

CompLaB also provides a C$++$ template (surrogateModel.hh) where users can incorporate a pre-trained surrogate model for calculating biogeochemical reactions, including artificial neural networks (ANNs). This functionality can be used to replace FBA, which requires solving many computationally expensive linear optimization problems (Sect. 3.2.2). In the example shown in Sect. 5, CompLaB provides local metabolite concentrations and biomass densities as inputs and the surrogate model outputs microbial growth rate (γ) and uptake/production rates of metabolites (F). While our demonstration is based on ANN models, any pre-trained statistical surrogate model (e.g., De Lucia and Kühn, 2021) that describes the sources and sinks – or their parameterization – can be used to enhance computational efficiency and accommodate various user-specific needs.

Figure 2 Schematic representation of the cellular automaton method for biomass spreading (a) before and (b) after biomass redistribution. The solid, pore, and biomass are color-coded gray, white, and blue, respectively. Brown solid arrows indicate the randomly chosen path which excess biomass travels. The numbers in the biofilm grid cells are the Manhattan distances to the nearest pore grid cell. Because all the neighbors of the dark-blue grid cell are saturated with biomass, the excess biomass is first moved to a neighboring cell randomly chosen amongst those with the shorter Manhattan distance. This allocation process is repeated until the excess biomass is moved to the grid cell with the Manhattan distance of 1 (color-coded white) where the excess biomass first encounters unsaturated neighbor cells. The biomass holding capacity ($B_{max}-B$) of the randomly selected neighboring grid cell (outlined green) is first evaluated. If the excess biomass exceeds the holding capacity, the excess biomass is redistributed to the first chosen neighbor up to the maximum holding capacity and the remaining excess biomass is placed to the next available neighbor grid cell (outlined yellow), finalizing the cellular automaton algorithm. The dashed black horizontal lines indicate the minimum biomass level required for a grid cell to be designated as biofilm.

3.3 Biomass redistribution

To explicitly model the spatial biomass expansion, CompLaB utilizes a cellular automaton (CA) with a predefined maximum biomass density (B_max) based on the CA algorithm developed by Jung and Meile (2021). After updating local biomass densities, the CA algorithm checks at every time step if there is any grid cell exceeding B_max and redistributes the excess biomass ($B-B_{max}$) to a randomly selected neighboring grid cell. If the selected grid cell cannot hold the excess biomass, the first chosen grid cell is filled up to the maximum holding capacity (B_max, a value defined by the user), and then the remaining excess biomass is allocated to a randomly chosen second neighbor cell. If all the neighboring grid cells are saturated with biomass ($B\ge B_{max}$) and hence the excess biomass cannot be placed, the Manhattan distances of biomass grid cells to the closest pore cells are evaluated. The remaining excess biomass is then placed in a neighboring grid cell that is closer to pores. and this biomass allocation process is repeated until all the excess biomass is redistributed. Figure 2 shows an example of the cellular automaton process for biomass redistribution. Note that this biomass redistribution method is a simple approximation for biomass density conservation with room for improvement (e.g., Tang and Valocchi, 2013).

When the sessile biomass reaches a threshold density ($B\ge \mathit{\psi }{B}_{max}$, where ψ is a user-defined threshold biomass fraction; $\mathrm{0}\le \mathit{\psi }\le \mathrm{1}$), the pore grid cell is designated as biomass; if the biomass density falls below the threshold due to microbial decay or detachment (B < ψB_max), then a biomass grid cell is converted back to a pore. Pore grid cells allow for both advective and diffusive transport, but in the biofilm, sessile biomass hinders (i.e., permeable biofilm) or prevents (i.e., impermeable biofilm) flow, and the feedback between biomass growth/decay and advective flow conditions is accounted for by rerunning the flow solver to steady state after updating biomass distribution and corresponding pore geometry (Jung and Meile, 2021). The reduced advective transport efficiency in permeable biomass grid cells is implemented by modifying local fluid viscosity in the biofilm (ν_bf) with ${\mathit{\nu }}_{\text{bf}}={\mathit{\nu }}_{\mathrm{f}}/X$, where X is a user-defined viscosity ratio ($\mathrm{0}\le X\le \mathrm{1}$), while for impermeable biomass, a bounce-back condition is imposed (Pintelon et al., 2012). After imposing the new steady-state flow field, a streaming step of the transport solver is executed (Fig. 1).

4 Model verification

To verify the CompLaB implementation, the engineered metabolic interaction between two co-dependent mutant strains (Escherichia coli K12 and Salmonella enterica LT2) originally established by Harcombe (2010) and implemented in COMETS (Harcombe et al., 2014) and IndiMeSH (Borer et al., 2019) was chosen as a test case. E. coli K12 is deficient in producing methionine and relies on the release of methionine by the mutant S. enterica LT2. In turn, S. enterica LT2 requires acetate released by E. coli K12 because of its inability to metabolize lactose. As a result, these genetically engineered strains are obliged to engage in mutual interaction where neither species can grow in isolation. The ratio of the two strains converges to a stable relative composition after 48 h in all the in vitro and in silico experiments at both initial ratios of 1:99 and 99:1.

Both COMETS (Harcombe et al., 2014) and IndiMeSH (Borer et al., 2019) integrate flux balance cell models of these two microorganisms into a two-dimensional environment in which metabolites are exchanged via diffusion. The initial and boundary conditions of these simulations were mirrored in CompLaB, with 100 grid cells containing 3 × 10⁻⁷ g biomass each (total = 3 × 10⁻⁵ g biomass) distributed randomly across a two-dimensional domain of 25 × 25 grid squares. The grid length was set to 500 µm, and the initial distributions of the two species were allowed to overlap. Three replicate simulations were carried out for each initial microbial ratio of E. coli and S. enterica (1:99 and 99:1). For the exchange metabolites acetate and methionine, fixed boundary concentrations of 0 mM were imposed on the left and right side of the domain, respectively, with no flux conditions at the top and bottom boundaries. The concentrations of lactose (2.92 mM) and oxygen (0.25 mM) were imposed at all domain boundaries. Solute and biomass diffusion coefficients were fixed at 5 × 10⁻¹⁰ and 3 × 10⁻¹³ m² s⁻¹, respectively. Consistent with the IndiMeSH model implementation, the simulation was carried out using Eq. (9), a finite-difference method for biomass diffusion, and reduced metabolic models of E. coli K12 and S. enterica LT2 in which the number of metabolites and reactions of the original metabolic models was systematically reduced by 1 order of magnitude (Borer et al., 2019).

Figure 3 illustrates the average ratio of E. coli and S. enterica of all six simulations (triplicates for each initial composition ratios of 1:99 and 99:1) after 48 h of simulation time. CompLaB simulation results agree with both the observations and the model results of COMETS and IndiMeSH, demonstrating the metabolic inter-dependence of two strains, and the convergence to a stable composition ratio (Appendix A, Fig. A1).

Figure 3 The relative composition of six in vitro experiments and COMETS, IndiMeSH, and CompLaB simulation results after 48 h. Error bars represent 1 standard deviation.

5 Surrogate model integration

A major issue in fully coupling genome-scale metabolic networks to reactive transport models is the large computational demand due to the repeated calculation of the LP problems at every biomass grid cell and every time step. Previous studies have alleviated this issue by interpolating the solutions of LP problems from a lookup table generated in advance to a reactive transport simulation (Scheibe et al., 2009), dynamically creating a solution pool of the LP problems during the reactive transport process (Fang et al., 2011), or systematically reducing the size of the matrix encoding the metabolic network (Ataman et al., 2017; Ataman and Hatzimanikatis, 2017). Here, we introduce a statistical surrogate modeling approach using a pre-trained artificial neural network (ANN) following the approach presented in Song et al. (2023). A trained ANN model directly relates input parameters (i.e., uptake rates of substrates) to outputs (i.e., biomass and metabolite production rates) through a set of nonlinear algebraic equations. As computing such input–output relationships from a pre-trained ANN model is several orders of magnitude faster than running FBA using a fully fledged metabolic model, the use of such surrogate models to achieve a significant speed-up has attracted much attention recently (e.g., De Lucia and Kühn, 2021; Prasianakis et al., 2020).

Here, we use the metabolic network of Geobacter metallireducens GS-15 (iAF987; Lovley et al., 1993; http://bigg.ucsd.edu/models/iAF987, last access: September 2022), a strict anaerobe capable of coupling the oxidation of organic compounds to the reduction of metals such as iron and manganese, using ammonium as its nitrogen source to train an ANN model. The dataset used for training a surrogate ANN model was obtained by collecting FBA solutions of the base model obtained using the IBM ILOG cplex optimizer with the objective function chosen to maximize biomass production. The solution set was prepared by randomly sampling 5000 combinations of two growth-limiting metabolites – acetate (C_Ac) and ammonium ($C_{{\mathrm{NH}}_{\mathrm{4}}}$) – within the concentration ranges of $\mathrm{0}\le C_{\mathrm{Ac}}\le$ 0.5 mM and $\mathrm{0}\le C_{{\mathrm{NH}}_{\mathrm{4}}}\le$ 0.05 mM via Monte Carlo simulations. Substrate concentrations were converted to uptake fluxes via the parameter-based approach (Eq. 7) using the parameters from Fang et al. (2011) (V_max =  10 $\mathrm{mmol}\mathrm{acetate}{\mathrm{g}}_{\mathrm{DW}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ and 0.5 $\mathrm{mmol}\mathrm{ammonium}{\mathrm{g}}_{\mathrm{DW}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$; K_s = 0.01 mM acetate and 0.1 mM ammonium). These fluxes were used as lower bounds, and Fe³⁺ was allowed to be consumed without limitation. The collected FBA solution dataset was split and used to train (70 %), validate (15 %), and test (15 %), and we developed an ANN model using MATLAB's neural network toolbox. As key hyperparameters, the number of layers and the number of nodes in the ANN model were, respectively, determined to be 4 and 10 through grid search. Ensuring the accuracy of a surrogate model is of critical importance in reactive transport models because even small errors accumulating over successive time steps can lead to a substantial error. In this simple example, the trained ANN estimates the biomass growth rate of the full FBA model almost perfectly (R > 0.999) against training, validation, and test datasets. This shows that the fully fledged metabolic model can be replaced by the surrogate ANN model without substantial loss of accuracy, boosting the simulation speed as shown next.

6 Model performance

CompLaB inherits the massive scalability of Palabos which decomposes a simulation domain into multiple subdomains and assigns them to individual computational nodes. In the following, the scalability of various components of CompLaB is assessed for a simplified microbial dynamics problem.

Figure 4 (a) Simulation domain of the model performance simulation, with the steady-state flow velocity distribution and (b) an enlarged subsection showing the distribution of material numbers around a solid.

Table 1 Parameters used for all the model performance simulations.

6.1 Test case

The simulation domain (Fig. 4) was prepared by taking a subsection of 500 × 300 square elements from the porous medium of Souzy et al. (2020). The length of each element was 16.81 µm, and material numbers were assigned to solid (0), pore (1), and solid–pore interfaces (solid side of interface: 2; pore side of interface: 3). In total, 10 % of the interface grid cells (pore side) were randomly assigned as sessile biomass grid cells (4) initially. Flow was induced from left to right by imposing a fixed pressure gradient between left- and right-side boundaries, and no flow conditions were set at the top and bottom boundaries as well as on the grain surfaces. The steady-state flow field was then provided to the CompLaB transport solver for mass transport and reaction simulations (Péclet number is 1 for a characteristic length scale of 2 mm). Two growth-limiting metabolites – acetate (CH₃COO⁻) and ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) – were considered for the mass transport simulations. Acetate was injected at the inlet (left) boundary with the fixed concentration of 0.45 mM to the simulation domain initially filled with the same concentration. Ammonium concentration in the inflowing fluid and initially in the domain was 0 mM, but it was produced at solid surfaces assuming a zeroth-order mineralization rate (Table 1). For both metabolites, no gradient conditions were imposed at the outlet, top, bottom, and grain boundaries. No external source and initial planktonic biomass were assumed, so that all planktonic biomass is detached sessile biomass.

The biogeochemical problem is described by the following set of ADREs:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial C_{\mathrm{Ac}}}{\partial t}}=\mathrm{\nabla }\cdot \left(D_{\mathrm{Ac}}\mathrm{\nabla }{C}_{\mathrm{Ac}}\right)-\mathit{u}\cdot \mathrm{\nabla }{C}_{\mathrm{Ac}}-F_{\mathrm{Ac}}{B}_{\mathrm{S}},\text{(11)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial C_{{\mathrm{NH}}_{\mathrm{4}}}}{\partial t}}=& \mathrm{\nabla }\cdot \left(D_{{\mathrm{NH}}_{\mathrm{4}}}\mathrm{\nabla }{C}_{{\mathrm{NH}}_{\mathrm{4}}}\right)-\mathit{u}\cdot \mathrm{\nabla }{C}_{{\mathrm{NH}}_{\mathrm{4}}}+k_{{\mathrm{NH}}_{\mathrm{4}}}{\mathit{\delta }}_{\mathrm{S}}\\ & -F_{{\mathrm{NH}}_{\mathrm{4}}}{B}_{\mathrm{S}},\end{array}\text{(12)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial B_{\mathrm{P}}}{\partial t}}=& \mathrm{\nabla }\cdot \left(D_{{\mathrm{B}}_{\mathrm{P}}}\mathrm{\nabla }{B}_{\mathrm{P}}\right)-\mathit{u}\cdot \mathrm{\nabla }{B}_{\mathrm{P}}-{\mathit{\mu }}_{\mathrm{B}}{B}_{\mathrm{P}}-k_{\text{att}}{\mathit{\delta }}_{\mathrm{B}}{B}_{\mathrm{P}}\\ & +k_{\text{det}}{B}_{\mathrm{S}},\end{array}\text{(13)}& {\displaystyle \frac{\partial B_{\mathrm{S}}}{\partial t}}=\mathit{\gamma }{B}_{\mathrm{S}}-{\mathit{\mu }}_{\mathrm{B}}{B}_{\mathrm{S}}+k_{\text{att}}{\mathit{\delta }}_{\mathrm{B}}{B}_{\mathrm{P}}-k_{\text{det}}{B}_{\mathrm{S}}.\end{array}$$

C_Ac and $C_{{\mathrm{NH}}_{\mathrm{4}}}$ are the concentrations of acetate and ammonium, respectively; B_P is the planktonic biomass density; B_S is the sessile biomass density; u denotes the flow field; δ indicates the presence (δ = 1) or absence (δ = 0) of a grain surface (δ_S) and a biomass grid cell (δ_B); and D_Ac, $D_{{\mathrm{NH}}_{\mathrm{4}}}$, and D_B are the diffusion coefficients of acetate, ammonium, and planktonic biomass, respectively, which differ between pore, biomass, and solids. For simplicity, the diffusivities of all the metabolites and planktonic biomass were set to 10⁻⁹ m² s⁻¹ in the pores, 8 × 10⁻¹⁰ m² s⁻¹ in biomass grid cells, and 0 in the solid. Biomass attachment and detachment (k_att and k_det), biomass decay (μ_B), and organic matter mineralization ($k_{{\mathrm{NH}}_{\mathrm{4}}}$) were simulated using the reaction kinetics solver, with the corresponding rate constants summarized in Table 1. The simulation assumed that G. metallireducens grows only on solid surfaces and planktonic biomass attaches only to existing surface-attached aggregates (Grinberg et al., 2019).

The computational efficiency and parallel performance of CompLaB was tested by executing four independent simulations utilizing each reaction solver. The cell-specific metabolic fluxes (F, Eq. 9) and biomass growth rates (γ, Eq. 12) were calculated either through FBA (CPY or GLPK), ANN, and/or reaction kinetics (KNS). The KNS solver was combined with other solvers (CPY, GLPK, and ANN) or used as a stand-alone reaction solver with the cellular automaton algorithm invoked (CA) and B_max set to 100 g_DW L⁻¹. The model performance simulations with the FBA solvers (CPY and GLPK) were prepared with the same conditions used for training the ANN model (Sect. 5). The pre-trained ANN model from Sect. 5 was used for the separate simulation ANN. For the CA simulation, we create a situation similar to the above examples but in which substantial biofilm growth over a short simulation time is artificially induced. To that end, F and γ were computed as

$$\begin{array}{}\text{(14)}& {\displaystyle }F=k_{\text{kns}}\left({\displaystyle \frac{C_{\mathrm{Ac}}}{C_{\mathrm{Ac}}+K_{\mathrm{Ac}}}}\right)\left({\displaystyle \frac{C_{{\mathrm{NH}}_{\mathrm{4}}}}{C_{{\mathrm{NH}}_{\mathrm{4}}}+K_{{\mathrm{NH}}_{\mathrm{4}}}}}\right),\text{(15)}& {\displaystyle }\mathit{\gamma }=Y{F}_{\mathrm{Ac}},\end{array}$$

where F denotes the metabolic fluxes (for simplicity assuming F_Ac = $F_{{\mathrm{NH}}_{\mathrm{4}}}$), γ is the biomass growth rate, k_kns = 2.5 × 10⁻⁶ mM s⁻¹ (Marozava et al., 2014) is the maximum uptake rate, and K_Ac = 0.1 mM and $K_{{\mathrm{NH}}_{\mathrm{4}}}$ = 0.01 mM are the half-saturation constants for acetate and ammonium, respectively. The growth yield Y is set to 40 000 ${\mathrm{g}}_{\mathrm{DW}}{\mathrm{mmol}}_{\mathrm{Ac}}^{-\mathrm{1}}$, an arbitrarily large number used only to invoke the CA algorithm within 10 000 time steps (440 s). The flow field was updated every 10 time steps when the CA algorithm was invoked (as, e.g., Thullner and Baveye, 2008; Jung and Meile, 2021).

The performance tests were carried out on the computing nodes of the Georgia Advanced Computing Resources Center. Each node has an AMD EPYC 7702P 64-core processor with a 2.0 GHz clock cycle (AMD Rome), with 128 GB of RAM. The nodes are interconnected via an EDR InfiniBand network with 100 GB s⁻¹ effective throughput and run a 64-bit Linux operating system (CentOS 7.9 distribution). The elapsed (wall-clock) time for 10 000 time steps was recorded.

Figure 5 The wall-clock time recorded in seconds for (a) the total, (b) flow and transport, (c) all reactions, and (d) each reaction part of the CompLaB algorithm. Four different simulations were carried out deploying each reaction solver: three simulations used COBRApy (CPY), GLPK, artificial neural network (ANN) solvers for microbial metabolism calculation, and a kinetic solver (KNS) for the attachment/detachment of biomass (AT-DT). One simulation used only a kinetic solver for both microbial metabolism and AT-DT with the cellular automaton algorithm invoked (CA). The simulation time for CA only can be inferred by subtracting the time for AT-DT with CA (dashed green line with orange triangle symbols) from the time for KNS including CA (solid green line with gray triangle symbols) in (d). Each symbol represents the average of two simulations. The simulations exhibiting no or limited speed-up with 128 and 192 cores except the ANN simulation were excluded when drawing the regression lines. The negative numbers are the slopes of the solid regression lines. The average slope of all the dashed regression lines in panel (d) is −1.14.

6.2 Performance

A comparison of simulation times for flow, transport, and reaction shows that most compute time is used for simulating the reactions, in particular when integrating in silico cell models into a reactive transport framework (Fig. 5). This highlights the benefit of using efficient surrogate models. The surrogate ANN model substantially improves computational efficiency compared to CPY and GLPK (about 2 orders of magnitude in total elapsed time; Fig. 5a) because calculating the pre-trained ANN is much faster than solving the linear programming problem every time step in FBA (Fig. 5c). For reaction calculations (Fig. 5d), the ANN simulation (solid sky-blue line with gray square symbols) even exhibits comparable but slightly shorter simulation times than the traditional reaction kinetics calculation (KNS; dashed lines with orange symbols) because the ANN implementation only operates on biomass grid cells while KNS operates on both biomass and pore grid cells.

In addition to the computational efficiency, negligible errors introduced by the surrogate ANN model assure the use of a surrogate ANN model (Fig. B1). Although the error in biomass calculation accumulates over simulation time steps, it is kept to very low values throughout the simulation (the relative error is on the order of 10⁻⁹; Appendix B) and has practically no influence on metabolite concentration calculations. This observation illustrates that CompLaB can calculate microbial metabolic reactions in porous media with a heterogeneous distribution of pore, biomass, organic matter, and minerals, based on the genome-scale metabolic model with the superior computational efficiency of a surrogate model without losing accuracy. Applying Monte Carlo simulations to generate FBA data to train a neural network model required solving the linear programming problem 5000 times with a set of randomly chosen uptake rates of acetate and ammonium, which does not add a significant computational burden. In our application, we determined the number of layers and nodes in the neural network to be 4 and 10, respectively, because no further improvement in model performance was observed beyond those values.

The computational efficiency of ANN also works in favor of scalability. The reaction processes of CompLaB are inherently an embarrassing parallel task because calculating biogeochemical sub-processes is completely independent of the neighbors (except CA) and all model performance simulations show the reasonable scalability up to 64 cores (Fig. 5). However, the scaling behaviors of all simulations except ANN illustrate suboptimal scalability with no or limited speed-up when using more compute resources. The loss of efficiency originates mostly from the calculation of reaction processes (Fig. 5c) because the domain decomposition applied to the heterogeneous simulation domain (Fig. 4) resulted in an uneven distribution of biofilm grid cells per subdomain and hence a variable size of the problem to be solved on each core. In fact, in our simple example problem (a total of 500 × 300 computational grid cells, constant random seeds), a domain decomposition when using 64, 128, and 192 cores produced 6, 38, and 76 subdomains, respectively, that have no initial biomass grid cells. Such subdomains do not contribute to computing microbial metabolisms (FBA) and biomass redistribution (CA) that consumes most of the computational cost (e.g., GLPK calculation consumes ∼ 99 % of the total computational cost), preventing ideal parallel performance (Fig. 5d). While this is also true for ANN, computational efficiency of ANN reduces the time wasted by such subdomains. As a result, the suboptimality is not readily apparent in ANN (Fig. 5c).

The CA algorithm implemented in CompLaB is a nonlocal process requiring information from neighboring grid cells, but the result still exhibits a good scalability when using up to 64 cores and suboptimal scalability when more cores are used. The CA simulation required less time than FBA simulations because the metabolic reactions were calculated using KNS. But CA spent extra time on updating the flow field (Fig. 5b) and redistributing excess biomass (Fig. 5c and d). This illustrates that the actual time required for the CA algorithm depends on the nature of the biomass expansion. For example, more time will be required for a system with rapid biofilm growth excess because excess biomass has to travel a longer distance through a thick biofilm. Furthermore, flow fields will need to be updated more often to reflect the influence of rapid biofilm growth on flow.

7 Conclusions

The numerical modeling platform CompLaB for simulating 2D pore-scale reactive transport processes is capable of utilizing the quantitative implementation of the microbial metabolism through the coupling of genome-scale metabolic models. The integration of in silico cell models with reactive transport simulations makes this framework broadly applicable and enables the integration of knowledge gained from the “omics”-based characterization of microbial systems. For example, the successful reproduction of experimentally observed convergences to a stable composition of a two-species consortium (S. enterica and E. coli) demonstrates the capability of CompLaB to investigate metabolic interactions between multiple microbial species.

Our novel numerical framework based on the lattice Boltzmann method allows simulating advection as well as diffusion of metabolites in complex porous media. A wide range of simulation domains can be used to represent soil structure and fractured rock images directly obtained from various imaging techniques (e.g., μ-CT, FIB-SEM) or numerically generated porous media with material numbers assigned to pore, solid, and source/sink grid cells for biogeochemical reactions which include but are not limited to biofilms. The inherent parallel efficiency of CompLaB facilitates simulating dynamic flux balance analysis capturing the microbial feedback on flow and transport in porous media, as done previously using Monod-type representations of microbial activity (Jung and Meile, 2021). Furthermore, the versatile simulation environment of CompLaB allows utilizing surrogate models, such as an artificial neural network. The resulting speed-up enables the investigation of complex biogeochemical processes in natural environments.

Appendix A: Convergence of the verification model to a stable ratio after 100 h

The six simulation cases used in Sect. 4 for model verification were run for 100 h of simulation time to further evaluate whether the observed convergence to an average composition ratio is stable (Fig. A1). The composition ratio observed after 48 h (0.75) is largely maintained through the extended simulation period (increases only to 0.78 after 100 h).

Figure A1 The evolution of the fraction of E. coli relative to the total cell density (E. coli + S. enterica) over 100 h. Dashed and dotted lines denote an initial abundance of 99 and 1 % E. coli, respectively.

Appendix B: Surrogate model accuracy

The surrogate modeling approach inevitably introduces errors in model estimations. The errors should be maintained sufficiently low throughout the surrogate simulation otherwise errors can propagate in successive iterations and result in unphysical results. To quantify the errors in surrogate model estimations, the solutions of our artificial neural network (ANN) were compared to the reference simulation COBRApy (CPY) by calculating the arithmetic mean of the root mean squared errors:

$$\begin{array}{}\text{(B1)}& {\text{error}}_t={\displaystyle \frac{\mathrm{1}}{m}}\sum _j^m{\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{n_j}{\sum }_i^{n_j}({\text{CPY}}_{i,j}-{\text{ANN}}_{i,j}{)}^{\mathrm{2}}}}{max({\text{CPY}}_j{)}_t}},\end{array}$$

where j is the variable type (B_S, B_P, C_Ac, and $C_{{\mathrm{NH}}_{\mathrm{4}}}$), m is the number of variables j used in calculating the error, n_j is the number of grid cells for each variable j, and t is the time step where the error is evaluated. The differences between the FBA simulations using GLPK and CPY solvers are negligible, so only the CPY solution was chosen as the reference (Fig. B1).

Figure B1 The discrepancy (relative error) between the surrogate ANN and the fully fledged physical model simulation using COBRApy calculated via Eq. (A1). The surrogate model overestimates biomass (B_S+B_P), and the errors accumulate over time. But the errors are kept low and negligible throughout the simulation, as evidenced by no influence on metabolite concentrations ($C_{\mathrm{Ac}}+C_{{\mathrm{NH}}_{\mathrm{4}}}$).