A simulation method for reactive solute transport in groundwater based on microbial metabolic constraints

Abstract

This invention relates to the field of hydrogeology and discloses a method for simulating reactive solute transport in groundwater based on microbial metabolic constraints. The method includes: collecting microbial genetic information to construct genome-scale metabolic models of at least two functional bacterial communities; introducing thermodynamic feasibility constraints to each metabolic model, calculating the Gibbs free energy changes of key metabolic reactions and determining flux feasibility to obtain a thermodynamically constrained feasible region for metabolic flux; establishing a substrate competition mechanism among different bacterial communities based on energy yield within the feasible region, allocating upper limits for substrate uptake flux, and calculating the metabolic flux distribution of each bacterial community through dynamic flux balance analysis; converting the metabolic flux distribution into solute reaction rates and coupling it into a reactive solute transport model to simulate microbially driven biogeochemical processes. This invention characterizes the energy competition and regulation of microorganisms at the metabolic mechanism level, improving simulation accuracy and model transferability.

Description

Technical Field

[0001] This invention relates to the field of hydrogeology, and in particular to a method for simulating groundwater reactive solute transport based on microbial metabolic constraints. Background Technology

[0002] The migration and transformation of pollutants in groundwater systems is one of the core issues in hydrogeological research. Pollutants in groundwater (such as nitrates, sulfates, and heavy metals) are not only controlled by physical processes such as convection and dispersion during their flow, but also affected by microbial-mediated biogeochemical reactions.

[0003] Microorganisms drive the transformation and degradation of pollutants by breaking down organic matter or inorganic electron donors through metabolic activities, transferring electrons to electron acceptors (such as nitrates, sulfates, and oxygen). This process involves complex metabolic networks, interactions among multiple microbial communities, and dynamic changes in environmental conditions (such as substrate concentration, redox potential, and pH). To quantitatively describe the migration and transformation patterns of pollutants and provide a scientific basis for groundwater remediation, researchers typically use reactive solute transport models (RTMs) for numerical simulations.

[0004] In existing technologies, reactive solute transport models typically employ reaction rate expressions based on empirical kinetics to describe microbial-mediated chemical reactions. For example, patent document CN117473884A discloses a "numerical simulation method and system for reactive solute transport in groundwater in a multiphase system." This method constructs a model coupling multiphase flow, solute transport, and chemical reaction processes through mass and energy conservation equations to simulate the migration and transformation of inferior components (such as arsenic, fluorine, and iodine) between gas, liquid, and solid phases. In this method, the kinetic reaction rate is described using kinetic equations based on rate constants, such as the mineral dissolution / precipitation kinetic reaction rate expressed as r_k = S_area·V0·(V / V0)^θV·Σ(1-Is)·θ_j,k·k_j, where k_j is the reaction rate constant. This type of method has been widely used in traditional geochemical reaction simulations and can characterize simple chemical equilibrium and reaction kinetic processes.

[0005] However, this method has limitations: reaction rate parameters (such as half-saturation constant and maximum reaction rate) are usually obtained through experimental fitting or by searching historical data, exhibiting strong scenario dependence. These parameters need to be recalibrated when environmental conditions (such as temperature, pH, and substrate type) change. More importantly, this type of method cannot reflect the metabolic regulation mechanisms of microorganisms under different environmental conditions. For example, how microorganisms allocate electron donors based on energy gain, how they selectively respond among different electron acceptors, and how they compete for limited substrates when multiple microbial communities coexist. These mechanistic deficiencies limit the model's predictive ability for microbially driven biogeochemical processes in complex groundwater environments.

[0006] In recent years, with the development of molecular biology techniques, some studies have attempted to incorporate microbial functional gene information into reactive solute transport models, forming so-called "gene center kinetic models." These methods attempt to reflect the impact of changes in microbial community structure on reaction rates by incorporating functional gene abundance or enzyme concentration into the reaction rate expression. However, these methods are essentially still based on empirical kinetic expressions, and their core still relies on a pre-defined reaction rate function form, failing to accurately describe the material transformation pathways and electron transport processes in microbial metabolic networks. Furthermore, resource allocation between different metabolic pathways still needs to be set through empirical parameters, making it difficult to explain the competitive relationships between different reaction pathways at the metabolic mechanism level. Another type of research uses flux balance analysis (FBA) to construct microbial metabolic network models to predict the allocation of metabolic fluxes in microorganisms under specific environmental conditions. However, these metabolic models are usually not coupled with groundwater solute transport processes, failing to characterize the dynamic feedback relationship between metabolic activity and changes in solute concentration in the environment, thus limiting their application in real groundwater systems. Summary of the Invention

[0007] The purpose of this invention is to propose a simulation method for reactive solute transport in groundwater based on microbial metabolic constraints. This method addresses the shortcomings of existing techniques, which cannot characterize the energy competition strategies of microorganisms with limited substrates in multi-microbial communities, cannot automatically simulate the switching of metabolic pathways under changing environmental conditions (such as electron acceptor depletion), and suffer from poor model parameter transferability due to empirical fitting. Therefore, there is an urgent need for a simulation method that can deeply integrate microbial metabolic mechanisms with groundwater solute transport processes.

[0008] Specifically, this invention provides a method for simulating groundwater reactive solute transport based on microbial metabolic constraints, the method comprising the following steps: S1. Collect microbial gene information in the target groundwater environment, and construct genome-scale metabolic models of at least two functional bacterial communities based on the microbial gene information to obtain a set of multi-community metabolic models. S2. For the key metabolic models in the multi-microbial community metabolic model set, introduce thermodynamic feasibility constraints, calculate the Gibbs free energy change of the key metabolic reaction, and determine the feasibility of the metabolic flux based on the Gibbs free energy change to obtain the feasible region of metabolic flux after thermodynamic constraints. S3. Within the feasible region of metabolic flux after the thermodynamic constraints, a substrate competition mechanism between different bacterial communities is established based on energy yield. The upper limit of substrate uptake flux is allocated according to the energy yield of each bacterial community using the same substrate. The metabolic flux distribution of each bacterial community under dynamic environmental conditions is calculated through dynamic flux balance analysis. S4. The metabolic flux distribution is converted into solute reaction rate and coupled into a reactive solute transport model. By iteratively updating the environmental solute concentration and microbial metabolic flux, the microbial-driven biogeochemical processes in the groundwater system are simulated.

[0009] A storage medium storing instructions and data for implementing a groundwater reactive solute transport simulation method based on microbial metabolic constraints.

[0010] A groundwater reactive solute transport simulation device based on microbial metabolic constraints includes: a processor and a storage medium; the processor loads and executes instructions and data in the storage medium to implement a groundwater reactive solute transport simulation method based on microbial metabolic constraints.

[0011] The beneficial effects of this invention are as follows: By collecting microbial gene information to construct a multi-community genome-scale metabolic model, this invention introduces thermodynamic feasibility constraints (based on Gibbs free energy calculation) and a substrate competition mechanism based on energy yield. It couples dynamic flux balance analysis (dFBA) with a reactive solute transport model (RTM), achieving dynamic iterative feedback between microbial metabolic processes and solute transport processes. This method characterizes the energy allocation and competition strategies of microorganisms at the metabolic mechanism level, can automatically simulate the switching of metabolic pathways, and the model parameters have physical meaning (such as thermodynamic parameters) and strong mobility, significantly improving the simulation accuracy and predictive ability of microbial-driven biogeochemical processes in complex groundwater environments. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the method flow of the present invention; Figure 2 This is a schematic diagram of the hardware device operation according to an embodiment of the present invention. Detailed Implementation

[0013] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0014] Before formally describing the present invention, a general description of the solution of the present invention will be given first to facilitate understanding.

[0015] Example 1 Please refer to Figure 1 The present invention provides a method for simulating groundwater reactive solute transport based on microbial metabolic constraints, comprising the following steps: S1. Collect microbial gene information in the target groundwater environment, and construct genome-scale metabolic models of at least two functional bacterial communities based on the microbial gene information to obtain a set of multi-community metabolic models. As one embodiment, step S1 is implemented as follows: S11. When the target microorganism is a known microorganism, obtain the complete genome information of the target microorganism from a public gene database, identify metabolism-related genes through gene function annotation, and construct the metabolic response network of each microbial community. Specifically, depending on the research subject, the collection of microbial genetic information falls into two categories: Scenario 1: Known microbial system When the research subject is a known microorganism, the complete genome information of the target microorganism is obtained from a public gene database. Commonly used public gene databases include, but are not limited to: NCBI (National Center for Biotechnology Information) GenBank database, JGI (Joint Genome Institute) IMG / M database, EBI ENA database, etc.

[0016] For example, if the research target is nitrate-reducing bacteria commonly found in groundwater environments, Pseudomonas spp. can be selected. Pseudomonas Representative strains, such as Pseudomonas stutzeri The complete genome sequence was downloaded from the NCBI database (e.g., GenBank accession number CP002881.1). Subsequently, gene function annotation tools (e.g., Prokka, RAST, KEGG Automatic Annotation Server, etc.) were used to annotate the genome sequence, identifying metabolism-related genes, including but not limited to: genes involved in carbon metabolism (e.g., the gltA gene encoding citrate synthase), genes involved in nitrogen metabolism (e.g., the narG, narH, and narI genes encoding nitrate reductase), and genes involved in energy metabolism (e.g., the atpA-H gene encoding ATP synthase).

[0017] Scenario 2: Unknown microorganisms or complex environmental microbial communities When the research subject is an unknown microorganism or a microbial community in a complex environment, metagenomic sequencing technology is used to obtain the microbial genetic information. The specific steps are as follows: 1) Collect microbial samples from the target groundwater environment and perform metagenomic sequencing using high-throughput sequencing technology (such as the Illumina NovaSeq platform) to obtain raw reads.

[0018] 2) Quality control of raw sequencing reads: Use tools such as Trimmomatic or FastP to remove low-quality reads, adapter sequences, and reads with an excessively high proportion of N to obtain high-quality sequencing data.

[0019] 3) Assemble high-quality sequences: Use metagenomic assembly tools such as MEGAHIT or metaSPAdes to assemble sequences into longer genomic fragments (contigs).

[0020] 4) Obtaining individual microbial genomes through binning: Using binning tools such as MetaBAT, MaxBin, or CONCOCT, contigs are clustered into different genome bins based on characteristics such as GC content, coverage depth, and tetranucleotide frequency. Each genome bin corresponds to the genome set of a microbial species. The integrity and contamination level of each genome bin are assessed by checking single-copy marker genes (e.g., using the CheckM tool), and qualified microbial genomes are screened out.

[0021] 5) Perform functional annotation on the screened genome set: Use tools such as Prokka or DRAM to identify metabolism-related genes in each genome, including genes encoding enzymes, transport proteins, transcription factors, etc.

[0022] S12. When the target microorganism is an unknown microorganism or a complex environmental microbial community, the original sequencing reads are obtained through metagenomic sequencing, quality control and sequence assembly are performed, multiple microbial genome sets are obtained through binning, functional annotation is performed on each genome set, metabolism-related genes are identified, and metabolic response networks of each microbial community are constructed. Based on the metabolism-related genes identified in the above steps, genome-scale metabolic models of each functional bacterial community are constructed. The specific steps are as follows: 1) Establish a set of metabolic reactions based on the correspondence between genes, enzymes, and reactions. Specifically, using gene annotation information, through databases such as KEGG, MetaCyc, or BiGG, genes are mapped to corresponding enzymes (EC numbers), and then enzymes are mapped to specific metabolic reactions. For example, if the nitrate reductase gene narG is annotated in a microbial genome, then the corresponding metabolic reaction "nitrate + 2H+" is mapped to this gene. + +2e- → Nitrite + H2O is incorporated into the metabolic reaction set of this bacterial community.

[0023] 2) Construct metabolic response networks in the KBase platform (Department of Energy Systems Biology KnowledgeBase) or the RAVEN toolkit. The KBase platform provides a complete set of metabolic model building tools. Users can import functional annotation results into KBase and use its automatic modeling function to generate genome-scale metabolic models.

[0024] S13. Export the constructed metabolic response networks of each microbial community as metabolic model files in SBML format to form a multi-microbial community metabolic model set. Specifically, export the model file: export the constructed metabolic model as a model file in SBML (Systems Biology Markup Language) format. This model file contains the following information: a list of metabolic reactions (including reaction equations, reaction directions, and stoichiometric coefficients), a list of metabolites, gene-reaction associations, and upper and lower limits of reaction flux (default setting is -1000 to 1000 mmol / gDW / h, which can be dynamically adjusted according to environmental conditions).

[0025] Repeat the above steps to construct genome-scale metabolic models for at least two functional bacterial communities, forming a multi-community metabolic model set. For example, for nitrate-contaminated groundwater systems, nitrate-reducing bacteria (such as...) can be constructed. Pseudomonas stutzer ) and sulfate-reducing bacteria (such as Desulfovibrio vulgaris A metabolic model of two bacterial communities was used to simulate the competition between the two communities for electron donors.

[0026] As one embodiment, step S1 is implemented as follows: Example: Construction of metabolic models for known microbial systems To construct Pseudomonas ( Pseudomonas stutzeri Taking a metabolic model as an example: 1) Genome Acquisition: Downloaded from the NCBI GenBank database Pseudomonas stutzeri The complete genome sequence of A1501 (GenBank accession number: CP002881.1, size approximately 4.6 Mbp).

[0027] 2) Gene Function Annotation: Automated annotation is performed using the RAST (Rapid Annotation using Subsystem Technology) tool. RAST can identify coding sequences in the genome and annotate genes to specific functional subsystems based on sequence similarity alignment (BLAST) results. Annotation results include: Metabolism-related genes: such as the gltA gene (locus tag: PST_0001) which encodes citrate synthase, and the narG gene (PST_1234) which encodes nitrate reductase.

[0028] Transporter protein genes: such as the actP gene that encodes the acetate transporter.

[0029] Energy metabolism genes: such as the atpA-H gene cluster encoding ATP synthase.

[0030] 3) Metabolic Reaction Network Construction: The RAST annotation results are imported into the KBase platform. Using KBase's "BuildMetabolic Model" function, a metabolic reaction network is automatically generated based on the gene-enzyme-reaction mapping relationship. This network contains approximately 1200 metabolic reactions and 800 metabolites.

[0031] 4) Model Validation and Correction: The integrity of the model is verified by checking whether it can simulate growth under minimal culture conditions. If some essential reactions are found to be missing (e.g., unable to synthesize a certain amino acid), the missing reactions are automatically filled in using KBase's "GapFilling" function.

[0032] 5) Export model file: Export the validated model as an SBML format file (e.g., "Pseudomonas_stutzeri.xml"), which contains complete stoichiometric matrix, reaction list, gene-reaction associations and other information.

[0033] Example: Construction of a metabolic model for an unknown microbial system Taking the construction of a metabolic model of complex microbial communities in groundwater as an example: 1) Sample collection and sequencing: Water samples were collected from the target groundwater contaminated site, DNA was extracted, and metagenomic sequencing was performed using the Illumina NovaSeq 6000 platform to obtain approximately 100 Gb of raw sequencing data.

[0034] 2) Quality control: Trimmomatic software was used for quality control to remove low-quality reads (Phred score <20), connector sequences, and reads with a length of less than 100 bp, resulting in approximately 85 Gb of high-quality clean data.

[0035] 3) Sequence assembly: Metagenome assembly was performed using metaSPAdes, with k-mer sizes set to 21, 33, 55, 77, 99, and 121, resulting in a total contig length of approximately 120 Mbp and an N50 of 5.2 kb.

[0036] 4) Binning: Contigs were binned using MetaBAT 2, and clustered based on sequence GC content, coverage depth, and tetranucleotide frequency characteristics to obtain 45 genome bins. CheckM was used to evaluate the integrity of each genome bin, and 12 high-quality genome bins with integrity >70% and contamination <5% were selected.

[0037] 5) Functional annotation: The selected genome bins were functionally annotated using DRAM (Distributed Refinery for Annotated Metagenomes). DRAM integrates multiple databases such as KEGG, Pfam, and dbCAN, and can identify metabolism-related genes.

[0038] Example of annotation results: Box 1 (presumably Pseudomonas): contains nitrate reductase genes such as narG, narH, and narI, as well as denitrification-related gene clusters.

[0039] Box 2 (presumably of *Desulfovibrio*): contains sulfate reductase genes such as dsrA and dsrB.

[0040] 6) Metabolic model construction: For each high-quality genome box, a metabolic model was independently constructed in the KBase platform, and finally, metabolic models of 12 bacterial communities were obtained, forming a multi-community metabolic model set.

[0041] S2. For the key metabolic models in the multi-microbial community metabolic model set, introduce thermodynamic feasibility constraints, calculate the Gibbs free energy change of the key metabolic reaction, and determine the feasibility of the metabolic flux based on the Gibbs free energy change to obtain the feasible region of metabolic flux after thermodynamic constraints. It should be noted that the introduction of thermodynamic feasibility constraints in step S2 further includes: S21. Obtain the concentrations of reactants and products in the environment, and calculate the real-time Gibbs free energy change ΔG of the key metabolic reaction, where ΔG = ΔG' + RT ln(Q), ΔG' is the standard Gibbs free energy change, R is the ideal gas constant, T is the absolute temperature, and Q is the reaction entropy. Specifically, for key metabolic responses in the metabolic model, the actual Gibbs free energy change ΔG under current environmental conditions is calculated. The formula for calculating ΔG is as follows: ΔG = ΔG' + R·T· ln(Q) in: ΔG': Standard Gibbs free energy change (kJ / mol), the standard value under the conditions of pH=7, temperature 25°C, and ionic strength 0.1M. This value can be obtained from databases such as the eQuilibrator database, the NIST thermodynamics database, or literature. For example, the reaction "nitrate + 2H+" + + 2e - → The standard Gibbs free energy of nitrite + H2O becomes ΔG°' = -163.0 kJ / mol.

[0042] R: Ideal gas constant, with a value of 8.314 J / (mol·K).

[0043] T: Absolute temperature (K), determined based on the actual temperature of the groundwater. For example, if the groundwater temperature is 15°C, then T = 288.15 K.

[0044] Q: Reaction quotient is defined as the ratio of the product of product concentrations to the product of reactant concentrations, with each concentration expressed as an activity. For the reaction "aA + bB → cC + dD", the reaction quotient is: Q = [C] c [D] d / ([A] a [B] b ) Where [ ] represents concentration.

[0045] When calculating Q, the required ambient solute concentration is obtained from the reactive solute transport model at the current time step. For example, assuming the current ambient conditions are nitrate concentration of 2.0 mmol / L, nitrite concentration of 0.5 mmol / L, and hydrogen ion concentration determined by pH (if pH=7, then [H+] mmol / L = 0.5 mmol / L), the concentration of hydrogen ions is determined by pH (if pH=7, then [H+] mmol / L = 0.5 mmol / L = 0.5 mmol / L). + ]=10 -7 (mol / L), electron activity is determined by redox potential. Substituting these concentrations into the formula for calculating Q yields the Q value.

[0046] S22. Set thermodynamic switching conditions: When ΔG>0, the reaction is determined to be a non-spontaneous reaction, and the flux of the metabolic reaction is forced to be 0; when ΔG is lower than the preset maintenance energy threshold, the reaction energy is determined to be insufficient to maintain the basic metabolism of microorganisms, and the flux of the metabolic reaction is forced to be 0. It should be noted that, based on the calculated ΔG value, the thermodynamic feasibility of key metabolic reactions is assessed, and thermodynamic switching conditions are set: 1) Non-spontaneous reaction constraint: When ΔG > 0, it indicates that the reaction is thermodynamically non-spontaneous and cannot proceed in the forward direction. In this case, the forward flux of the metabolic reaction is forced to be 0. For example, if a reaction has ΔG = +15.5 kJ / mol > 0 under current environmental conditions, then the upper limit of the metabolic flux for that reaction is set to 0, meaning that the reaction cannot occur under the current conditions.

[0047] 2) Maintenance Energy Constraint: Microorganisms need to maintain basic life activities during metabolism, including maintaining transmembrane ion gradients, repairing DNA, and synthesizing maintenance proteins. This energy requirement is called maintenance energy. When ΔG is lower than a preset maintenance energy threshold, it indicates that the energy released by the reaction is insufficient to support the maintenance needs of the microorganism. The maintenance energy threshold is usually set between -20 kJ / mol and -10 kJ / mol, and the specific value can be adjusted according to the type of microorganism and environmental conditions. For example, for typical chemoheterotrophic bacteria, the maintenance energy threshold can be set to -15 kJ / mol. When ΔG > -15 kJ / mol, the reaction flux is forced to be 0.

[0048] 3) Mathematical expression of thermodynamic switching: The above thermodynamic switching conditions can be expressed as linear constraints embedded in the solver of flux balance analysis. For the key metabolic reaction j, its flux v is set. j The following constraints must be met: If ΔG j If >0, then v j ≤ 0 (Only reverse reactions are allowed, but reverse reactions are usually not considered in the metabolic network, so the actual value is v) j = 0) If ΔG j <The sustaining energy threshold (e.g., -15 kJ / mol) then v j ≥ 0 (allowing positive reactions), and further limits its flux upper limit through other constraints. In practical solutions, a more common approach is to set both the upper and lower flux limits of a reaction to 0 for reactions where ΔG > 0 or ΔG > 0 sustaining energy thresholds, i.e., v j = 0. In this way, thermodynamically infeasible metabolic reactions are eliminated from the feasible solution space of the metabolic network.

[0049] S23. The thermodynamic switching conditions are embedded as linear constraints into the solver of flux balance analysis to obtain the feasible region of metabolic flux after thermodynamic constraints.

[0050] It should be noted that the aforementioned thermodynamic switching conditions are used as linear constraints to update the original flux feasible region for each microbial community metabolic model. The original flux feasible region is determined by the mass conservation constraint (S... v = 0) and reaction flux upper and lower limits constraints (v min ≤ v ≤ v maxThe new flux feasible region is composed of ( ). After thermodynamic constraints, the new flux feasible region is: {v | S v = 0, v min ' ≤ v ≤ v max ', and satisfies the thermodynamic switching condition} Among them, v min 'and v max 'These are the upper and lower limits of flux adjusted according to thermodynamic switching conditions.'

[0051] As one embodiment, the thermodynamic constraint process in step S2 can be described in detail through the following specific examples: Example: Thermodynamic constraints of nitrate reduction reaction Suppose we simulate a groundwater column experiment with the following initial conditions: nitrate concentration 2.0 mmol / L, nitrite concentration 0.1 mmol / L, pH=7.0, and temperature 15°C (T=288.15 K).

[0052] 1) Calculate ΔG for the nitrate reduction reaction. Reaction equation: NO3 - + 2H + + 2e - → NO2 - + H2O According to the eQuilibrator database, the standard Gibbs free energy change for this reaction is ΔG°' = -163.0 kJ / mol (at pH=7 and 25°C).

[0053] Calculate the reaction entropy Q: Q = [NO2] - ] / ([NO3 - ][H + ]²[e - ]²) Assuming that electron activity is determined by redox potential, and given that the groundwater environment is under moderately reducing conditions, let the electron activity [e] be... - ]=10 -6 M.

[0054] Substitute the value: [NO3] - ] = 2.0 × 10 - ³ M; [NO2] - ] = 1.0 × 10 -4 M; [H] + ] = 10 -7 M; [e] - ] =10 -6 M; but: Q = (1.0×10-4 ) / (2.0×10 - ³ × (10 -7 )² × (10 -6 )²)= (1.0×10 -4 ) / (2.0×10 - ³ × 10 - ¹ 4 × 10 - ¹²) = (1.0×10 -4 ) / (2.0×10 - ² 9 = 5.0 × 10² 4 ;; Calculate RT ln(Q): RT = 8.314 × 288.15 = 2395 J / mol = 2.395 kJ / mol; ln(Q) = ln(5.0×10² 4 ) = ln(5.0) + 24×ln(10) = 1.609 + 24×2.303 =1.609 + 55.272 = 56.881; RT ln(Q) = 2.395 × 56.881 = 136.2 kJ / mol; Therefore: ΔG = ΔG°' + RT ln(Q) = -163.0 + 136.2 = -26.8 kJ / mol; 2) Determination of thermodynamic switch The sustaining energy threshold was set at -15 kJ / mol. The calculated ΔG = -26.8 kJ / mol < -15 kJ / mol indicates that the reaction is thermodynamically feasible and the released energy is sufficient to support the sustaining needs of the microorganisms. Therefore, the forward flux of the reaction is permissible.

[0055] If environmental conditions change, for example, if nitrite accumulates to 5.0 mmol / L, the calculated ΔG = -15.2 kJ / mol, which is close to the sustaining energy threshold; if nitrite further accumulates to 10.0 mmol / L, ΔG may be greater than -15.0 kJ / mol. At this point, the thermodynamic switch will shut down the reaction, forcing the nitrate reduction flux to 0, simulating the automatic shutdown of the nitrate reduction pathway.

[0056] 3) Constraint Embedding The above thermodynamic switching conditions are transformed into linear constraints. Taking the Pseudomonas metabolic model as an example, this model includes the nitrate reduction reaction (reaction ID: R).nar The original flux range is [-1000, 1000] mmol / gDW / h. If the current time step calculates ΔG = -26.8 kJ / mol, the reaction is feasible, and the original flux range is maintained. If ΔG = -12.0 kJ / mol > -15.0 kJ / mol, then R... nar The flux upper limit is set to 0, i.e., v_nar ≤ 0. Since the model usually does not allow reverse reactions, in actual solutions, v nar = 0. This constraint is implemented using the changeRxnBounds function in the COBRA Toolbox.

[0057] S3. Within the feasible region of metabolic flux after the thermodynamic constraints, a substrate competition mechanism between different bacterial communities is established based on energy yield. The upper limit of substrate uptake flux is allocated according to the energy yield of each bacterial community using the same substrate. The metabolic flux distribution of each bacterial community under dynamic environmental conditions is calculated through dynamic flux balance analysis. It should be noted that the establishment of a substrate competition mechanism between different bacterial communities based on energy yield in step S3 further includes: S31. For each substrate, calculate the energy yield that each functional bacterial community can obtain through respiratory metabolism using the substrate, wherein the energy yield is the amount of ATP generated per mole of substrate. Specifically, for each substrate, the energy yield obtained by each functional bacterial community through respiratory metabolism using that substrate was calculated. Energy yield is defined as the amount of ATP generated by the complete oxidation of one mole of substrate (unit: mol ATP / mol substrate).

[0058] Energy yield calculations are based on ATP production responses in metabolic models. A typical method is to calculate the theoretical maximum ATP yield through flux balance analysis: setting microbial community growth as the objective function, the amount of ATP produced per mole of substrate consumed under optimal growth conditions is calculated. The specific formula is as follows: Y_ATP = ATP_generation_rate / substrate_uptake_rate in: ATP_generation_rate: The sum of fluxes of all ATP-generating reactions in the metabolic model, including substrate-level phosphorylation (such as pyruvate kinase reactions) and oxidative phosphorylation (ATP synthesis reactions coupled through the electron transport chain).

[0059] substrate_uptake_rate: Substrate uptake flux.

[0060] For example, for Pseudomonas ( Pseudomonas stutzeriThe denitrification reaction, using acetic acid as an electron donor and nitrate as an electron acceptor, can have its energy yield calculated using a metabolic model. With an acetic acid uptake flux of 1 mmol / g DW / h to maximize bacterial growth, the ATP production flux was obtained by solving for the free radical basis (FBA). The calculated Y_ATP was approximately 8-12 mol ATP / mol acetic acid, with the specific value depending on the detail of the metabolic model and parameter settings.

[0061] For example, for desulfurization vibrio ( Desulfovibrio vulgaris Using acetic acid as an electron donor and sulfate as an electron acceptor in sulfate reduction has a low energy yield, approximately 2-4 mol ATP / mol acetic acid. This is because the electron transport chain in sulfate reduction is short, resulting in low proton pumping efficiency.

[0062] S32. Rank the functional bacterial groups competing for the same substrate according to their energy yield from high to low, and the bacterial group with the highest energy yield shall have the priority to take up the substrate. Specifically, for multiple bacterial communities competing for the same substrate, they are ranked from highest to lowest energy yield. The bacterial community with the higher energy yield obtains greater energy benefits from the substrate and holds a dominant position in resource competition.

[0063] For example, in a system where acetic acid is the sole electron donor, *Pseudomonas* (Y_ATP≈10 mol ATP / mol acetic acid) and *Desulfovibrio* (Y_ATP≈3 mol ATP / mol acetic acid) compete for acetic acid. Based on energy yield, *Pseudomonas* has a higher energy yield and is therefore in a dominant position in the competition.

[0064] S33. Based on the relative ratio of energy yield of each bacterial community, dynamically adjust the upper limit of uptake flux of each bacterial community to the same substrate. The upper limit of uptake flux of bacterial community with high energy yield is set to the higher value of bacterial community with low energy yield, and the upper limit of uptake flux of bacterial community with low energy yield is compressed and lower than that of bacterial community with high energy yield. Specifically, the upper limit of uptake flux of each bacterial community for the same substrate is dynamically adjusted based on the relative ratio of energy yield of each community. The allocation rules are as follows: Let K be the set of bacterial communities competing for the same substrate, and let Y be the energy yield of each community. k (k∈K). Define the energy yield weighting factor: w k = Y k / Σ {j∈K} Y j The upper limit of substrate uptake flux for each microbial community is allocated according to a weighting factor: v max,k = v total × wk × α in: v total The total available uptake rate of the substrate in the environment (determined by substrate concentration and total microbial biomass). α: Competition intensity adjustment coefficient, 0 < α ≤ 1, which can be calibrated according to experimental data. Usually, α = 1 (allocated entirely according to energy yield) or α = 0.8 (considering competition resistance due to non-energy factors).

[0065] For example, assuming the initial concentration of acetic acid in the environment is 1 mM, and the total available uptake rate is 0.5 mmol / L / h, and the energy yield of *Pseudomonas* is 10, while that of *Desulfovibrio* is 3, then: w pseudo = 10 / (10+3) ≈ 0.769 w desulfo = 3 / (10+3) ≈ 0.231 The upper limit of acetic acid uptake flux for Pseudomonas is set at 0.5 × 0.769 = 0.3845 mmol / L / h, and the upper limit of acetic acid uptake flux for Desulfovibrio is set at 0.5 × 0.231 = 0.1155 mmol / L / h.

[0066] As environmental conditions change (such as nitrate depletion), the energy productivity of each bacterial community may change, and the upper limit of uptake flux will also be dynamically adjusted accordingly.

[0067] S34. Using the adjusted upper limit of uptake flux as the boundary constraint for dynamic flux balance analysis, solve for the metabolic flux distribution of each microbial community.

[0068] Specifically, based on the aforementioned thermodynamic and substrate competition constraints, dynamic flux balance analysis is used to calculate the metabolic flux distribution of each bacterial community under dynamic environmental conditions. The core idea of ​​dFBA is to solve a linear programming problem at each time step to obtain the optimal metabolic flux distribution under the current environmental conditions, then update the environmental concentration, and proceed to the next time step.

[0069] The specific implementation steps are as follows: 1) Set the time step: Set the total simulation time T and the discrete time step Δt. The selection of the time step needs to balance computational efficiency and numerical stability. Usually, Δt is set to 0.1-1 hours, and can be appropriately extended for slow processes.

[0070] 2) Initialize environmental concentrations: Set the initial concentration distribution of each solute (electron donor, electron acceptor, metabolite) and the initial biomass of each bacterial community.

[0071] 3) Execute the following sub-steps within each time step: a. Constructing constraints: For each bacterial community, construct the following set of constraints: Mass conservation constraint: S k · v k = 0; Flux upper and lower limits constraints: v min ,k ≤ v k ≤ v max,k ; Thermodynamic constraint: Based on the ΔG value calculated in step S2, the flux of infeasible reactions is set to 0; Substrate competition constraint: Limit the substrate uptake flux according to the uptake flux limit allocated in step S33; b. Define the objective function: Typically, the objective function is to maximize the biomass growth of each bacterial community. For community k, the objective function is: max Z k = μ k · X k Where μ k X is the specific growth rate (determined by biomass synthesis flux). k This represents the biomass of bacterial community k.

[0072] c. Solving the linear programming problem: Use the flux balance analysis solver in COBRA Toolbox (MATLAB environment) or COBRApy (Python environment) (such as using the linprog function or optimization solvers such as Gurobi and CPLEX) to solve the above linear programming problem and obtain the metabolic flux vector v of each bacterial community. k .

[0073] d. Calculate the reaction rate: Convert the substrate uptake flux and product formation flux in the metabolic flux vector into the macroscopic reaction rate. For example, if the uptake flux of a certain bacterial community of substrate S is v... S (Unit: mmol / gDW / h), then the consumption rate of substrate S by this bacterial community is v. S × X k (Unit: mmol / L / h), where X k The biomass concentration of bacterial community k is expressed in gDW / L.

[0074] e. Biomass renewal: based on specific growth rate μ k Updated microbial biomass: X k (t+Δt) = X k (t) + μ k × X k (t) × Δt 4) Proceed to the next time step: Transfer the calculated reaction rate and updated biomass to the reactive solute transport model (step S4) to update the environmental solute concentration, and then proceed to the next time step, repeating step 3).

[0075] As one embodiment, the substrate competition process in step S3 can be described in detail through the following specific examples: Example: Competitive allocation of acetic acid Consider a groundwater contamination system containing two functional bacterial communities: Group A: Pseudomonas ( Pseudomonas stutzeri It can utilize acetic acid as an electron donor and nitrate as an electron acceptor for denitrification.

[0076] Group B: Desulfurization Vibrio ( Desulfovibrio vulgaris It can use acetic acid as an electron donor and sulfate as an electron acceptor to reduce sulfate.

[0077] Acetic acid (initial concentration 5.0 mmol / L) in an environment where two bacterial communities compete.

[0078] 1) Calculate the energy yield of each bacterial community. Based on genome-scale metabolic models of various bacterial communities, flux balance analysis was used to calculate energy yield. The acetic acid uptake flux was set at 1 mmol / g DW / h, and the ATP production flux was calculated using biomass maximization as the objective function.

[0079] For Pseudomonas, metabolic model calculations show that for every 1 mmol of acetic acid consumed, approximately 10.2 mmol of ATP is generated (Y_ATP = 10.2 mol ATP / mol acetic acid).

[0080] For Desulfovibrio, metabolic model calculations show that for every 1 mmol of acetic acid consumed, approximately 3.5 mmol of ATP is generated (Y_ATP = 3.5 mol ATP / mol acetic acid).

[0081] 2) Energy Yield Ranking Sorted by energy yield from highest to lowest: Pseudomonas: 10.2 mol ATP / mol acetic acid (ranked 1st); Desulfuric vibrio: 3.5 mol ATP / mol acetic acid (ranked 2nd); The energy yield of Pseudomonas is about 2.9 times that of Desulfovibrio, giving it a clear advantage in the competition for acetic acid.

[0082] 3) Calculate the upper limit allocation of intake flux Assuming the current time step has an acetic acid concentration of 5.0 mmol / L in the environment, the total available uptake rate is determined by the acetic acid concentration and the total microbial biomass. Let the upper limit of the total acetic acid uptake flux be 0.5 mmol / L / h (this value is determined based on the acetic acid concentration and the transport capacity of the microorganisms).

[0083] Calculate the energy yield weighting factor: w_pseudo = 10.2 / (10.2 + 3.5) = 10.2 / 13.7 = 0.745; w_desulfo = 3.5 / (10.2 + 3.5) = 3.5 / 13.7 = 0.255; If we set the competition intensity adjustment coefficient α = 1.0, then: The upper limit of acetic acid uptake flux for Pseudomonas aeruginosa = 0.5 × 0.745 = 0.3725 mmol / L / h; Upper limit of acetic acid uptake flux by *Desulfovibrio* = 0.5 × 0.255 = 0.1275 mmol / L / h; 4) Constraint Implementation in dFBA In the COBRA Toolbox, competitive constraints are achieved by setting an upper limit on substrate uptake flux for each microbial community metabolic model.

[0084] Pseudomonas model: model_pseudo = changeRxnBounds(model_pseudo, 'EX_ac_e',-0.3725, 'l'); where the uptake flux is negative to set the lower limit.

[0085] Desulfurization Vibrio model: model_desulfo = changeRxnBounds(model_desulfo, 'EX_ac_e', -0.1275, 'l'); 5) Dynamic adjustment As time progresses, environmental conditions change. For example, when nitrates are depleted, Pseudomonas can no longer utilize nitrates as electron acceptors, and its energy yield drops significantly (to about 2.5 mol ATP / mol acetic acid). At this point, a recalculation of energy yield ranking and upper limit allocation of uptake flux reveals that Desulfovibrio may obtain more acetic acid, mimicking the process of bacterial succession.

[0086] It should be noted that the dynamic flux balance analysis described in step S3 further includes: Within each time step, the objective function is to maximize the biomass of each microbial community. The constraints are the steady-state constraints of the metabolic network stoichiometry matrix, the upper and lower limits of the reaction flux, and the thermodynamic and substrate competition constraints introduced in steps S2 and S3. The metabolic flux distribution of each microbial community is solved by linear programming.

[0087] As one embodiment, the dFBA solution process in step S3 can be described in detail through the following specific examples: Example: Solving the dFBA of a Pseudomonas metabolic model The metabolic model of Pseudomonas aeruginosa includes the following elements: number of metabolic reactions: 1250, number of metabolites: 950, and stoichiometric matrix S: 950 rows × 1250 columns.

[0088] 1) Constructing constraints The following is a partial code example of using the COBRA Toolbox to construct constraints for a linear programming problem in the MATLAB environment: % Load metabolic model model = readCbModel('Pseudomonas_stutzeri.xml'); % Set mass conservation constraints (automatically included in the model) % Set upper and lower limits for flux model.lb = -1000 * ones(1250, 1); % Lower limit (intake is negative); model.ub = 1000 * ones(1250, 1); % Upper limit (generates positive values); % Set substrate uptake constraints model = changeRxnBounds(model, 'EX_ac_e', -0.5, 'l'); % Upper limit of acetic acid intake: 0.5 mmol / g DW / h model = changeRxnBounds(model, 'EX_nitrate_e', -0.4, 'l');% Upper limit of nitrate intake: 0.4 mmol / g DW / h; % Set thermodynamic constraints (based on ΔG calculation results) model = changeRxnBounds(model, 'R_nar', 0, 'u');% If nitrate reduction is thermodynamically infeasible; model = changeRxnBounds(model, 'R_nar', 0, 'l'); 2) Define the objective function The objective function is set to maximize biomass. Metabolic models typically include a "biomass reaction," which synthesizes biomass by consuming various precursors (amino acids, nucleotides, lipids, etc.) in a fixed proportion. The objective function is to maximize the flux (v_biomass) of the biomass reaction. A partial code example is shown below: % Set the target function model.c = zeros(1250, 1); biomass_rxn_id = find(strcmp(model.rxns, 'BIOMASS')); model.c(biomass_rxn_id) = 1; % Maximize biomass 3) Solve the linear programming problem The solveLP function is used to solve linear programming problems. A partial code example is shown below: % Solve solution = solveLP(model, 'max'); % Extraction Results v_biomass = solution.x(biomass_rxn_id); % Specific growth rate v_ac = solution.x(find(strcmp(model.rxns, 'EX_ac_e')));% Acetic acid uptake flux v_nitrate = solution.x(find(strcmp(model.rxns, 'EX_nitrate_e')));% Nitrate uptake flux 4) Multi-microbial community joint solution For multi-community systems, a linear programming problem needs to be solved for each community separately, taking into account the competitive constraints between communities. The solution yields the metabolic flux distribution and specific growth rate for each community.

[0089] 5) Dynamic updates At the next time step, the substrate uptake constraint is updated based on the updated environmental concentration, and the above solution process is repeated.

[0090] S4. The metabolic flux distribution is converted into solute reaction rate and coupled into a reactive solute transport model. By iteratively updating the environmental solute concentration and microbial metabolic flux, the microbial-driven biogeochemical processes in the groundwater system are simulated.

[0091] It should be noted that step S4, which involves iteratively updating the environmental solute concentration and microbial metabolic flux, further includes: S41. In the current time step, convert the metabolic flux distribution of each microbial community calculated in step S3 into substrate consumption rate and product generation rate. Specifically, the metabolic flux distribution of each bacterial community calculated in step S3 is converted into the reaction rate term required by the reactive solute transport model.

[0092] For the i-th solute (such as nitrate, nitrite, acetic acid, etc.), the overall reaction rate R i It equals the sum of the net production rates of the solute by each bacterial community: R i = Σ {k=1}^{N} (v i,k × X k ) in: v i,k : In the metabolic flux vector of bacterial community k, the net flux of solute i (unit: mmol / g DW / h). If v i,k A positive value indicates that bacterial community k produces solute i; a negative value indicates that bacterial community k consumes solute i.

[0093] X k Biomass concentration of bacterial community k (unit: gDW / L).

[0094] N: Total bacterial count.

[0095] For example, suppose the metabolic flux distribution of *Pseudomonas* shows: nitrate uptake flux is -0.5 mmol / gDW / h (the negative sign indicates consumption), nitrite formation flux is +0.5 mmol / gDW / h, and biomass is 0.2 gDW / L; *Desulfovibrio* has no direct effect on nitrate. Then the overall nitrate reaction rate is: R nitrate = (-0.5) × 0.2 = -0.1 mmol / L / h; the total reaction rate of nitrite is: R nitrate = (+0.5) × 0.2 = 0.1 mmol / L / h.

[0096] S42. Substitute the substrate consumption rate and product formation rate as reaction terms into the reactive solute transport equation, and solve the environmental solute concentration distribution at the next time step using the PFLOTRAN computing platform. Specifically, the above reaction rate R i Substituting the reactive solute transport equations, the environmental solute concentration distribution at the next time step is obtained by solving the equations using the PFLOTRAN computing platform.

[0097] The reactive solute transport equation is described by the convection-dispersion-reaction equation: C i / t = ·(D· C i ) - ·(u· C i ) + R i in: C i : Concentration of the i-th solute (unit: mmol / L or mg / L); D: Dispersion coefficient tensor (unit: m² / s), D = α L ·|u| + D m α L For longitudinal dispersion, D m The molecular diffusion coefficient; u: Groundwater velocity vector (unit: m / s), obtained by solving the groundwater flow equation; R i : Reaction term (unit: mmol / L / s), calculated from step S41.

[0098] PFLOTRAN is a high-performance parallel computing platform for simulating reactive solute transport in two-dimensional and three-dimensional heterogeneous porous media, supporting multi-component reactive solute transport simulations. The solution process is as follows: 1) Mesh Discretization: Based on the geological characteristics of the study area, a three-dimensional stratigraphic model is established and meshed. The mesh size must ensure that it can capture key processes (such as the movement of the reaction front), and the mesh is usually densified in areas where geological conditions change drastically.

[0099] 2) Initial conditions setting: Set the spatial distribution of each solute concentration at the initial time, as well as the initial distribution of microbial biomass.

[0100] 3) Boundary condition setting: Set boundary conditions, including Dirichlet boundary (constant concentration boundary), Neumann boundary (constant flux boundary), or Cauchy boundary (mixed boundary), to simulate the exchange between groundwater and the external environment.

[0101] 4) Time Discreteness and Iterative Solution: An implicit time discrepancy method is used to transform the partial differential equations into a system of algebraic equations. Within each time step, the nonlinear equation system is solved using the Newton-Raphson iterative method. During the iteration process, the reaction rate R_i at the current concentration is substituted into the equations as a source term to obtain the new concentration distribution C. i new.

[0102] S43. Feed back the environmental solute concentration of the next time step to steps S2 and S3, and update the Gibbs free energy calculation parameters of each metabolic reaction and the upper limit of substrate uptake flux of each bacterial community. Specifically, the environmental solute concentration distribution obtained from the PFLOTRAN solution at the next time step is fed back to steps S2 and S3 to update the following parameters: In step S22, ΔG is calculated to determine the required concentrations of reactants and products; In step S33, the total available uptake rate required to allocate the upper limit of substrate uptake flux is determined. In step S34, the substrate concentration constraint required for dFBA solution is obtained.

[0103] For example, if PFLOTRAN calculates that the nitrate concentration decreases from 2.0 mmol / L to 1.5 mmol / L at the next time step, this concentration value is passed back to the dFBA calculation module to update the upper limit of nitrate uptake flux for Pseudomonas.

[0104] S44. Repeat steps S41 to S43 until the preset simulation termination time is reached, thereby realizing the dynamic iterative coupling of microbial metabolic processes and solute transport processes.

[0105] Specifically, steps S41 to S43 are repeated until the preset simulation termination time is reached. Through this iterative coupling method, dynamic feedback between microbial metabolic processes and solute transport processes is achieved: microbial metabolism changes environmental concentrations, and environmental concentrations, in turn, constrain microbial metabolic activities.

[0106] As one embodiment, the iterative coupling process in step S4 can be described in detail through the following specific examples: Example: Iterative simulation of an experiment on a groundwater column contaminated with nitrates.

[0107] A 1 m long sand column was set up for the experiment. Initial conditions: the column was filled with water, the initial concentration of nitrate was 2.0 mmol / L, the initial concentration of acetic acid was 5.0 mmol / L, and there was no nitrite. Water free of nitrate and acetic acid was injected into the left end of the column, and the right end was allowed to drain freely.

[0108] Time step settings: The total simulation time is set to 48 hours, and the time step Δt = 0.1 hours (6 minutes). In the initial stage of the reaction (the first 12 hours), due to the large concentration gradient and vigorous reaction, a smaller step size is used; when the reaction stabilizes in the later stage, the step size can be adaptively increased to 1 hour.

[0109] Iterative process: 1) Hour 0 (Initial Time Step): Initial concentration distributions were obtained from PFLOTRAN: nitrate 2.0 mmol / L, acetic acid 5.0 mmol / L, nitrite 0.

[0110] The concentration is passed to the dFBA calculation module.

[0111] 2) Hours 0 to 0.1 (first time step): dFBA calculation: Based on the current concentration, calculate the metabolic flux distribution of each bacterial community. The metabolic fluxes for Pseudomonas are shown as follows: acetic acid uptake flux 0.5 mmol / gDW / h, nitrate uptake flux 0.4 mmol / gDW / h, and nitrite formation flux 0.4 mmol / gDW / h.

[0112] Converting metabolic flux to reaction rates: With a Pseudomonas biomass of 0.1 gDW / L, the acetic acid consumption rate is 0.05 mmol / L / h, the nitrate consumption rate is 0.04 mmol / L / h, and the nitrite production rate is 0.04 mmol / L / h.

[0113] The reaction rate is passed to PFLOTRAN.

[0114] PFLOTRAN Solution: The implicit finite element method was used to solve the convection-dispersion-reaction equation, obtaining the concentration distribution after 0.1 hours. Calculation results: The nitrate concentration at the left end of the column decreased to 1.96 mmol / L, the acetic acid concentration decreased to 4.95 mmol / L, and the nitrite concentration increased to 0.04 mmol / L.

[0115] 3) Hours 0.1 to 0.2 (second time step): Obtain the updated concentration distribution from PFLOTRAN.

[0116] Repeat the above dFBA calculation process, but at this time the nitrite concentration increases, which may affect the thermodynamic constraint (the nitrite reduction reaction becomes energy-favorable).

[0117] The updated reaction rate is then passed to PFLOTRAN to solve for the concentration at the next time step.

[0118] 4) Continue iterating until 48 hours have passed: The coupling process described above is repeated at each time step.

[0119] When the nitrate concentration drops below a certain threshold, the thermodynamic constraint automatically shuts off the nitrate reduction pathway, simulating the termination of the denitrification process.

[0120] If sulfate is present, the metabolic activity of desulfovibrio may be activated after nitrate is depleted, mimicking a switch in metabolic pathways.

[0121] Convergence control: To ensure the numerical stability of the coupled computation, an internal iteration is set up within each time step. Specifically, in the Newton iteration of PFLOTRAN, the reaction rate is recalculated using dFBA after each concentration update until the concentration change is less than a preset tolerance (e.g., 1×10⁻⁶). -6 (mmol / L) or reach the maximum number of iterations (e.g., 10 times).

[0122] Results output: After the simulation, the distribution of nitrate, nitrite, and acetic acid concentrations at each time point and spatial location, as well as the biomass evolution curves of each bacterial community, will be output. These results can be used for comparison and verification with experimental observation data.

[0123] Example 2: Please see Figure 2 , Figure 2 This is a schematic diagram of the hardware device in operation according to an embodiment of the present invention. The hardware device specifically includes: a groundwater reactive solute transport simulation device 401 based on microbial metabolic constraints, a processor 402, and a storage medium 403.

[0124] A groundwater reactive solute transport simulation device 401 based on microbial metabolic constraints: The groundwater reactive solute transport simulation device 401 based on microbial metabolic constraints realizes the groundwater reactive solute transport simulation method based on microbial metabolic constraints.

[0125] Processor 402: The processor 402 loads and executes the instructions and data in the storage medium 403 to implement the groundwater reactive solute transport simulation method based on microbial metabolic constraints.

[0126] Storage medium 403: The storage medium 403 stores instructions and data; the storage medium 403 is used to implement the above-mentioned method for simulating groundwater reactive solute transport based on microbial metabolic constraints.

[0127] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for simulating reactive solute transport in groundwater based on microbial metabolic constraints, characterized in that: Includes the following steps: S1. Collect microbial gene information in the target groundwater environment, and construct genome-scale metabolic models of at least two functional bacterial communities based on the microbial gene information to obtain a set of multi-community metabolic models. S2. For the key metabolic models in the multi-microbial community metabolic model set, introduce thermodynamic feasibility constraints, calculate the Gibbs free energy change of the key metabolic reaction, and determine the feasibility of the metabolic flux based on the Gibbs free energy change to obtain the feasible region of metabolic flux after thermodynamic constraints. S3. Within the feasible region of metabolic flux after the thermodynamic constraints, a substrate competition mechanism between different bacterial communities is established based on energy yield. The upper limit of substrate uptake flux is allocated according to the energy yield of each bacterial community using the same substrate. The metabolic flux distribution of each bacterial community under dynamic environmental conditions is calculated through dynamic flux balance analysis. S4. The metabolic flux distribution is converted into solute reaction rate and coupled into a reactive solute transport model. By iteratively updating the environmental solute concentration and microbial metabolic flux, the microbial-driven biogeochemical processes in the groundwater system are simulated.

2. The method for simulating groundwater reactive solute transport based on microbial metabolic constraints as described in claim 1, characterized in that: The introduction of thermodynamic feasibility constraints in step S2 further includes: S21. Obtain the concentrations of reactants and products in the environment, and calculate the real-time Gibbs free energy change ΔG of the key metabolic reaction, where ΔG = ΔG' + R·T·ln(Q), ΔG' is the standard Gibbs free energy change, R is the ideal gas constant, T is the absolute temperature, and Q is the reaction entropy. S22. Set thermodynamic switching conditions: When ΔG>0, the reaction is determined to be a non-spontaneous reaction, and the flux of the metabolic reaction is forced to be 0; when ΔG is lower than the preset maintenance energy threshold, the reaction energy is determined to be insufficient to maintain the basic metabolism of microorganisms, and the flux of the metabolic reaction is forced to be 0. S23. The thermodynamic switching conditions are embedded as linear constraints into the solver of flux balance analysis to obtain the feasible region of metabolic flux after thermodynamic constraints.

3. The method for simulating groundwater reactive solute transport based on microbial metabolic constraints as described in claim 1, characterized in that, Step S3, which establishes a substrate competition mechanism among different bacterial communities based on energy yield, further includes: S31. For each substrate, calculate the energy yield that each functional bacterial community can obtain through respiratory metabolism using the substrate, wherein the energy yield is the amount of ATP generated per mole of substrate. S32. Rank the functional bacterial groups competing for the same substrate according to their energy yield from high to low, and the bacterial group with the highest energy yield shall have the priority to take up the substrate. S33. Based on the relative ratio of energy yield of each bacterial community, dynamically adjust the upper limit of uptake flux of each bacterial community to the same substrate. The upper limit of uptake flux of bacterial community with high energy yield is set to the higher value of bacterial community with low energy yield, and the upper limit of uptake flux of bacterial community with low energy yield is compressed and lower than that of bacterial community with high energy yield. S34. Using the adjusted upper limit of uptake flux as the boundary constraint for dynamic flux balance analysis, solve for the metabolic flux distribution of each microbial community.

4. The method for simulating groundwater reactive solute transport based on microbial metabolic constraints as described in claim 1, characterized in that: Step S4, which iteratively updates the environmental solute concentration and microbial metabolic flux, further includes: S41. In the current time step, convert the metabolic flux distribution of each microbial community calculated in step S3 into substrate consumption rate and product generation rate. S42. Substitute the substrate consumption rate and product formation rate as reaction terms into the reactive solute transport equation, and solve the environmental solute concentration distribution at the next time step using the PFLOTRAN computing platform. S43. Feed back the environmental solute concentration of the next time step to steps S2 and S3, and update the Gibbs free energy calculation parameters of each metabolic reaction and the upper limit of substrate uptake flux of each bacterial community. S44. Repeat steps S41 to S43 until the preset simulation termination time is reached, thereby realizing the dynamic iterative coupling of microbial metabolic processes and solute transport processes.

5. The method for simulating groundwater reactive solute transport based on microbial metabolic constraints as described in claim 1, characterized in that: The construction of a genome-scale metabolic model of at least two functional bacterial communities described in step S1 further includes: S11. When the target microorganism is a known microorganism, obtain the complete genome information of the target microorganism from a public gene database, identify metabolism-related genes through gene function annotation, and construct the metabolic response network of each microbial community. S12. When the target microorganism is an unknown microorganism or a complex environmental microbial community, the original sequencing reads are obtained through metagenomic sequencing, quality control and sequence assembly are performed, multiple microbial genome sets are obtained through binning, functional annotation is performed on each genome set, metabolism-related genes are identified, and metabolic response networks of each microbial community are constructed. S13. Export the constructed metabolic response networks of each microbial community as metabolic model files in SBML format to form a multi-microbial community metabolic model set.

6. The method for simulating groundwater reactive solute transport based on microbial metabolic constraints as described in claim 1, characterized in that: The dynamic flux balance analysis described in step S3 further includes: Within each time step, an objective function is set, such as maximizing the biomass of each bacterial community. The metabolic flux distribution of each bacterial community is solved by linear programming, with the steady-state constraints of the stoichiometric matrix of the metabolic network, the upper and lower limits of the reaction flux, as well as the thermodynamic constraints and substrate competition constraints introduced in steps S2 and S3 as constraints.

7. The method for simulating groundwater reactive solute transport based on microbial metabolic constraints as described in claim 1, characterized in that, The reactive solute transport model described in step S4 uses a convection-dispersion-reaction equation to describe the solute transport process: C i / t = ·(D· C i ) - ·(u· C i ) + R i Where C i Let be the concentration of the i-th solute, D be the dispersion coefficient tensor, u be the groundwater velocity vector, and R be the groundwater velocity vector. i The reaction term is obtained by converting the metabolic flux distribution of each bacterial community calculated in step S3.

8. A storage medium, characterized in that: The storage medium stores instructions and data to implement the groundwater reactive solute transport simulation method based on microbial metabolic constraints as described in any one of claims 1 to 7.

9. A groundwater reactive solute transport simulation device based on microbial metabolic constraints, characterized in that: include: A processor and a storage medium; the processor loads and executes instructions and data in the storage medium to implement the groundwater reactive solute transport simulation method based on microbial metabolic constraints as described in any one of claims 1 to 7.

Images