The genome of S. oneidensis MR-1 was obtained from GeneBank (GCA_000146165.2) and annotated using RAST, version 2.0 (Overbeek et al., 2005). The draft model was constructed using ModelSEED, version 2.0 (Overbeek et al., 2014). Annotation information from KEGG (Kanehisa et al., 2016), UniProt (Consortium, 2007), DeepEC (Ryu et al., 2019) and iMR1_799 was used to extended the model. For genes that were not included in the draft model but were annotated with EC numbers in KEGG, UniProt, DeepEC, or iMR1_799, reaction information used to supplement GPRs was extracted from the ModelSEED reaction database according to the corresponding EC numbers.

The biomass composition was adapted from the experimental data of S. oneidensis MR-1, which is consistent with that used in iSO783 (Pinchuk et al., 2010). To simulate biomass synthesis, the cell growth and non-growth associated ATP maintenance values (GAM and NGAM, respectively) were assumed to be identical to those used in iSO783. Gapfilling was carried out under the conditions of an l-lactate uptake rate of 4.06 mmol/gDW/h, while exchange reactions of other basic nutrient components (including H₂O, O₂, NH₃/NH₄ ⁺, sulfite, phosphate and metal ions) were not restricted. The gapfilling algorithm developed based on weight-added pFBA was as follows: Min ∑ ( a b s ( v i ) + a b s ( y j ) ) (1) W ⋅ S ⋅ v + U ⋅ y = 0 (2) 0 < v t arg e t ≤ 10 (3) v l b ,   i ≤ v i ≤ v u b ,   i (4) y l b ,   j ≤ y j ≤ y u b ,   j (5) where S is the stoichiometric matrix representing the model to be gap-filled, v is the vector of fluxes through reactions in S , U is the stoichiometric matrix representing the universal ModelSEED database from which reactions are activated to fill gaps, y is the vector of fluxes through reactions in U , v t a r g e t is the flux of demand reaction for which a precursor could not be synthesized, v l b ,   i and v u b ,   i are respectively the lower and upper bounds of fluxes through reactions in the original model, W is the weight for reactions in the original model (set to 1,000), while y l b ,   j and y u b ,   j are respectively the lower and upper bounds of fluxes through reactions in ModelSEED.

The gap filling process consists of three steps. Firstly, identifying biomass compositions that cannot be synthesized. The demand reaction of each biomass composition was added and set as the optimization target for FBA calculation. If the result is 0, it means that the biomass composition cannot be synthesized, and there are gaps in its synthesis pathway. Then, the gaps for each composition were filled using the weight-added pFBA gapfilling algorithm (Eqs 1–5). The algorithm is implemented by adding 1,000 weights to all reactions in the Model_02 or Model_03 to minimize the number of gapfilling reactions introduced from ModelSEED. Finally, all the reactions used to fill the gaps are added into the Model_02 or Model_03 to enable the synthesis of all the biomass compositions.

The wild-type strain S. oneidensis MR-1 used in this study was obtained from our laboratory. The strains were stored at −80°C and revived by growing at 30 °Con LB agar plates containing (per liter): 10 g tryptone, 5 g yeast extract, and 10 g NaCl. The Biolog phenotype microarray was used to evaluate the utilization capabilities of different carbon, nitrogen, sulfur and phosphorus sources in MR-1. The PM1, PM2A, PM3B, PM4A (A-E), PM4A (F-H) plates consist of 190 carbon sources, 95 nitrogen sources, 59 phosphorus sources and 35 sulfur sources, respectively. All these plates were inoculated with cell suspensions at 100 μL/well. The plates were placed at 30°C and monitored for 72 h with readings taken at 15 min intervals. The kinetic information was recorded and quantified using OmniLog OL_FM_12 kinetic software (Biolog, United States) followed by data analysis. In the event of substrate utilization, photographic readings of color intensity resulted from dye reduction were represented in OmniLog units (OU). The substrates that were 1.5-fold higher or 0.8-fold lower than the control were selected as positive or negative results, respectively.

Gapfilling was conducted for those substrates that could be used in the Biolog experiment but not in the model. Two categories were found: 1) Substrates that lack the exchange reactions (From extracellular to intracellular). These gaps were filled manually; 2) Substrates that lack utilization pathway. These gaps were filled using the gap-filling method described above.

The bounds of exchange reactions were manually set to minimum medium conditions, with unrestricted H₂O, O₂, NH₃/NH₄ ⁺, sulfite, phosphate and metal ions. The uptake rate of the chosen carbon source was set to 4.06 mmol/gDW/h, while rates for other carbon sources were set to zero. For the simulation of terminal electron acceptors, ATPM (“rxn00062_c0”) was set as the objective function with an l-lactate uptake rate of 4.06 mmol/gDW/h. All the exchange reactions for terminal electron acceptors were turned off, but those for the chosen one ware not restricted. The metabolic model testing suite (MEMOTE) was used to evaluate the consistency and annotation of the model (Lieven et al., 2020). The validation of essential genes consists of two steps. Firstly, for each gene, the associated reactions are removed from the model; Then, biomass was set as objective function with an l-lactate uptake rate of 4.06 mmol/gDW/h. If the result is 0, it means that the gene is essential for growth.

Flux balance analysis (FBA) (Orth et al., 2010) is used for optimizing a pre-defined objective function in the specified metabolic constraints. Unless stated otherwise, maximization of the biomass production rate was considered as the objective function. Parsimonious flux balance analysis (pFBA) (Lewis et al., 2010) was chosen for analyzing the product biosynthesis pathway. These methods were conducted using the COBRApy (0.18.1) (Ebrahim et al., 2013). The optimization solvers GLPK and Gurobi were used for linear and quadratic programming. Flux scanning based on the enforced objective flux (FSEOF) algorithm was used to identify the metabolic modification targets (Choi et al., 2010).

Coulombic efficiency (CE) is commonly used to quantify the performance of MFCs and is defined as the ratio of electrons actually transferred to the anode to that theoretically present in the starting carbon source (Mao and Verwoerd, 2014). For electrons exported to the anode, four reactions were added into the model to represent the electron production process. The “Electrode_Demand” reaction was set as objective to represent the generation of all electrons. Thus, CE could be represented by: CE % = C o u t p u t C s u b s t r a t e × 100 % = F ( m m o l / g D W / h ) × 2 × 100 s u b s t r a t e   u p t a k e   r a t e   ( m m o l / g D W / h ) × N × % (6) where F is the flux of ‘Electrode_Demand’, N is the theoretical maximum number of electrons donated by the carbon source. It can be calculated from the reduction degree of the substrate (Dugar and Stephanopoulos, 2011). For example, N for lactate is C₃H₆O₃ = 3*4 + 6*1 + 3*(-2) = 12.

The theoretical maximum yield was calculated according to the ratio of reduction degrees of substrate and product. Theoretical maximum yield ( mol/mol ) = C s u b s t r a t e C p r o d u c t = N M (7) where N is the reduction degree per mole of substrate, M is the reduction degree per mole of product.

To reconstruct the genome-scale metabolic network of S. oneidensis MR-1, the genome of MR-1 was annotated by RAST, then uploaded into ModelSEED to generate a draft metabolic model. Charges and formulas of metabolites provided by ModelSEED were added into the draft model. The temporary gene IDs (e.g., 211586.77. peg.1131) in the draft model were transformed into specific gene IDs of S. oneidensis MR-1 (e.g., SO_1198) by local Blastp which was utilized for homology matching under the criteria: e-value < 1E-30; identity = 100%. The results of Blastp are listed in Supplementary Datasheet S1 (Supplementary Table S1). Furthermore, all the exchange reaction in the draft model was not constraint. Thus, the bounds of exchange reactions were manually set to minimum medium conditions, while l-lactate (4.06 mmol/gDW/h) was used as sole carbon and energy source. The draft model contained 866 genes, 1,436 metabolites and 1,524 reactions (Figure 1), which was more than the published models iSO783 and iMR1_799. For instance, dihydrofolate reductase was annotated with SO_3646 and SO_0772 in the draft model, but it was only annotated with SO_3646 in iMR1_799.

However, the draft model could only cover 20.6% (866/4196) of the annotated ORFs. Some genes, which have been reported/annotated with clear metabolic functions, were not annotated by RAST. For example, SO_4706 was annotated as malate synthase (EC: 2.3.3.9) in KEGG, but it was missed in the draft model. In order to compensate for these genes and their corresponding reactions, manual extension of the draft model was performed based on the public biochemical databases KEGG and UniProt, in combination with the published model iMR1_799. Moreover, further protein annotation from DeepEC (Ryu et al., 2019), which is a deep learning-based computational tool for prediction of EC numbers, was also used to extend the draft model. A total of 265 genes, 217 reactions and 148 metabolites were extracted to extend the draft model and generate the Model_01. All these genes and reactions can be tracked in Supplementary Datasheet S2 (Supplementary Table S2).

As shown in Figure 1, the ATP production rate (up to the upper bound of the boundary 1,000) simulated by Model_01 were incorrect. Therefore, Model_01 had to be revised. Analysis of the ATP synthesis pathway revealed that there was infinite energy-generation caused by incorrect reaction directions in the Model_01. The wrong direction of some reactions led to the infinite generation of ATP or proton driving force, resulting in energy-generating cycles (Figure 2). Therefore, the MetaCyc database (Caspi et al., 2019) and iML1515 model (Monk et al., 2017) were used as references to correct the reaction directions in Model_01. The directions of 92 reactions were modified, and energy-generating cycles were avoided (Supplementary Datasheet S1, Supplementary Table S3).

After these modifications, the optimal generation rate of ATP was changed from 1,000 to 31.3 mmol/gDW/h. The P/O of the aerobic respiratory chain in Model_01 was 1.5, which is lower than the 2.5 from iSO783 based on literature (Toledo-Cuevas et al., 1998). Thus, all aerobic respiration chains were curated based on iSO783. Consequently, a 54.5 mmol/gDW/h ATP generation rate was obtained. Furthermore, six anaerobic respiratory chains were modified and 11 anaerobic respiratory chains were added (Supplementary Datasheet S1, Supplementary Table S3). After these corrections, Model_02 was obtained (Figure 1). A total of 16 terminal electron acceptors of S. oneidensis MR-1, with known electron transport pathways confirmed by experimental data, were chosen for in silico simulations (Table 1). If the flux of energy generation is greater than 0, the terminal electron acceptor was regarded as a positive result in silico. The agreement of growth between the simulated results of iLJ1162 and in vivo results achieved 100% (16/16) of all tested electron acceptors. Nevertheless, 62.5% (10/16) and 56.25% (9/16) concordance was obtained by iMR1_799 and iSO783, respectively. This can be explained by the key reactions missing in these models. For instance, the electron transfer reaction for the conversion of AsO⁴⁻ was missing in iMR1_799 and iSO783, but was successfully annotated to be associated with genes SO_4591, SO_1777, SO_1778 and SO_1779 in iLJ1162. Notably, when using metal ions such as Fe³⁺ and Co³⁺ as terminal electron acceptors, S. oneidensis MR-1 obtained the highest energy synthesis efficiency in silico. This may be due to the reduction of these electron acceptors to generate the proton driving force and resulting energy production.

However, the simulated growth rate of Model_02 was still zero, while that the reported experimentally measured growth rate of S. oneidensis MR-1 cultured in minimal medium was 0.085 h⁻¹, and the corresponding l-lactate consumption rate was 4.06 mmol/gDW/h under aerobic conditions (Pinchuk et al., 2010). This can be explained by the inability of Model_02 to synthesize a number of biomass precursors, such as lipopolysaccharide and glycogen. Therefore, the weight-added pFBA algorithm was used to fill the gaps in Model_02 (Materials and Methods). A total of 14 reactions (Supplementary Datasheet S1, Supplementary Table S4) were added into Model_02 to generate Model_03. Under the same l-lactate uptake rate, the predicted cell growth rate of iSO783 was 0.108 h⁻¹ while that of iMR1_799 was 0.192 h⁻¹. However, Model_03 predicted a rate of 0.105 h⁻¹, which was the closest to the actual in vivo value of 0.085 h⁻¹ (Supplementary Datasheet S1, Supplementary Figure S1).

In order to improve the substrates utilization capacity of Model_03, 190 carbon sources, 95 nitrogen sources, 59 phosphorus sources and 35 sulfur sources were systematically investigated in Biolog phenotype microarray experiments (Supplementary Datasheet S1, Supplementary Table S5). The Biolog results showed that 27 carbon sources, 32 nitrogen sources, 47 phosphorus sources and 14 sulfur sources can be used by S. oneidensis MR-1. However, 54 of these could not be utilized by Model_03 due to missing transport or utilization reactions (Supplementary Datasheet S1, Supplementary Table S6). After supplementing 187 reactions obtained by the weight-added pFBA algorithm, iLJ1162 was generated (Supplementary Datasheet S2). A total of 289 substrates were used to evaluate the substrate utilization abilities in different models. 94 true positive (substrates that can be used in both the in vivo experiments and simulations) and 157 true negative (substrates that cannot be used in both the in vivo experiments and simulations) results were obtained by iLJ1162 (Figure 3). This means iLJ1162 has an agreement of 86.9% with the experimental data, which is significantly higher than the 62.3 and 61.9% obtained by iMR1_799 and iSO783, respectively (Figures 3B,C). This can be attributed to the extensively extended GPRs associated information and gapfilling curation for these substrates (Supplementary Datasheet S1, Supplementary Table S6). However, due to current limited understanding of metabolism of some substrates, there are still some inconsistencies of in vivo results with those obtained by iLJ1162. For instance, the utilization pathways of Tween 20, 40, and 80 as well as gelatin, are still unknown, and these cannot be simulated in silico.

The final version of the model, iLJ1162, includes 1,162 genes, 1818 metabolites and 2084 reactions, which were divided into the three extracellular, periplasmic and cytosolic compartments. MEMOTE was used to assess the quality of iLJ1162 and the other two published GEMs of MR-1 (Supplementary Datasheet S3). In addition to the significant improvements in consistency as well as annotation of reactions and metabolites, iLJ1162 has clear annotations for genes and SBO-terms (Systems Biology Ontology) (Figure 4). These result in a clear advantage of iLJ1162 in the total score. All validation results confirmed the high reliability of iLJ1162, indicating that it can be used for further analysis.

Essential genes are an important characterization of model quality (Lee et al., 2014; Nouri et al., 2020). Therefore, iLJ1162 was used to analyze the essentiality of individual genes of MR-1 under aerobic conditions (Materials and Methods). From the report by Yang et al. (2014), we collected 781 unambiguous gene inactivation results whose gene IDs could be matched to the three MR-1 models. After comparison of in vivo results with in silico simulations, iLJ1162 achieved an accuracy of (681/781) 87.2%, significantly higher than (471/781) 60.3% of iSO783 and (477/781) 61.1% of iMR1_799 (Supplementary Datasheet S4, Supplementary Table S7). This can be majorly attributed to the significantly improved ORFs coverage in iLJ1162. However, there were still some inconsistent results, including 59 false negative results (essentiality for experiment, but non-essentiality for iLJ1162) and 35 false positive results (essentiality for iLJ1162, but non-essentiality for experiment). For these false negative results, it is probably caused by the cellular transcriptional regulation in MR-1. For example, SO_3415 was annotated as aspartate kinase and assessed as essential gene from the in vivo experiments. However, in the iLJ1162, there exist other three isoenzymes (SO_4055, SO_3427 and SO_3986), which might be strictly regulated and unexpressed in MR-1. As for false-positive results, it is likely caused by the incomplete annotation (coverage of ORFs, 27.8%) in iLJ1162. For example, SO_4234 was annotated as 3-isopropylmalate dehydratase and participates an essential role in l-leucine synthesis in all the three models. If the results of the experiment are accurate, MR-1 may have potential isoenzymes for 3-isopropylmalate dehydratase or alternative pathways for leucine synthesis.

In MFCs, electrons obtained through the oxidation of the substrates are transferred to electrodes via the electron transport chain (Figure 5). The electrode is actually one of terminal electron acceptors, and can be set as a pseudo-metabolite which is similar to metal ions such as Fe³⁺ (Mao and Verwoerd, 2014). Thus, four artificial reactions (Electrode_Sink, → Oxidated_Electrode; Electrode_Reduction_1, 2.0 Oxidated_Electrode + mql7 → 2.0 Reduced_Electrode + Menaquinone_7 + 2.0 H_plus__e0; Electrode_Reduction_2, 2.0 Oxidated_Electrode + methymenaquinol_7 → 2.0 H_plus__e0 + 2.0 Reduced_Electrode + methylmenaquinone_7; Electrode_Demand, Reduced Electrode →) were added into the model for simulation of electron production (Supplementary Datasheet S4, Supplementary Table S8). The “Electrode_Demand” reaction was set as the objective function to represent the generation of electrons.

To simulate the electron generation ability of MR-1 with different substrates, 15 available carbon sources experimentally validated using Biolog PM1-2 and five heterogenous carbon sources widely used as laboratory substrates were chosen for computational simulation. Taking both the biomass synthesis and electron generation into consideration, the biomass growth rate was set to 10% of the maximum theoretical value, then the maximization of electron generation was set as objective function. This is consistent with that MR-1 grows slowly in MFCs under anaerobic conditions (Logan et al., 2019).

In general, the simulated electron generation rate per C-mole substrate (Figure 6A) was positively correlated with the inherent electron number per C-mole substrate (Figure 6B). Glycerol was predicted to be the optimal substrate for electron generation with the highest rate (4.57 mmol per C-mol glycerol), which was ∼15% higher than that of the common laboratory common substrates l-lactate or N-Acetyl-d-glucosamine (NAG). Many substrates obtained high CE values over 95%, which means that the electrons of these substrates can practically to be fully transferred to the electrode by the metabolic network of MR-1 (Figure 6C). However, some substrates such as adenosine and inosine generated low C-mole electron generation rates despite their high inherent electron numbers (Figures 6A,B), which led to low CEs values from 40 to 70%. This is attributed to the production of intermediate metabolites during the metabolism of these substrates, such as hypoxanthine production for adenosine and inosine utilization, or succinate production for the utilization of the dipeptide glycyl-l-proline. However, these intermediate metabolites cannot be further metabolized by MR-1 under anaerobic conditions, and have to be transported outside of the cell, causing electron waste.

The enhancement of the regeneration of reducing equivalents, mostly in the form of NADH, has been reported to be one of the most successful metabolic engineering strategies for the improvement of MFCs (Li et al., 2018). In view of the important role of NADH regeneration in MFCs, the iLJ1162 model was used to predict potential modification targets for increasing the NADH regeneration capacity of S. oneidensis MR-1 using the FSEOF algorithm (Choi et al., 2010). An artificial reaction (NADH_G, NADH_c → NAD⁺_c + H⁺_c) was introduced into the model and used as the optimization objective to calculate the optimal pathways for NADH regeneration from NAD⁺. In addition, since there is low or no measurable activity of the pyruvate dehydrogenase (PDH) encoded by SO_0424 and SO_0425 in S. oneidensis MR-1 under anaerobic conditions (Yang et al., 2006), the reactions corresponding to PDH were temporarily turned off. In the simulation, the reactions with significant flux changes from optimal biomass synthesis to optimal NADH regeneration were chosen as potential targets. Two common laboratory substrates for Shewanella species, l-lactate and NAG were introduced as simulated substrates (Supplementary Datasheet S4, Supplementary Table S9). In the modification strategies obtained from the FSEOF simulation with l-lactate as substrate, 17 overexpression targets were mainly distributed in the serine cycle (Figures 7A,B). Among them, NAD⁺ dependent malate dehydrogenase (Mdh), pyruvate formate-lyase (Pfl) and NAD⁺ dependent formate dehydrogenase (Fdh) were consistent with previously reported experimental results (Li et al., 2018), which means that the modification strategies predicted by iLJ1162 are meaningful. Thus, other targets may be potentially effective and worth implementing to improve electron production in vivo. In the case of NAG (Figures 7A,C), the simulation results revealed that 31 potential targets, mainly distributed in the PP pathway (e.g., Rpe and RpiA) and serine bypass (e.g., SerA and Mdh).

An important approach for increasing EET rates based on high intracellular regeneration of NADH is the transfer of electrons from NADH to the respiratory chain. The first step of the respiratory chain is the dehydrogenation of NADH by NADH dehydrogenases (NDHs) to release electrons, which are important for energy supply and cell growth (Figure 8A). Three types of prokaryotic NADH dehydrogenases are generally recognized in the literature (Melo et al., 2004), including those with proton-pumping activity (type I), proton/sodium-pumping activity (type II) and sodium-pumping activity (type III). Notably, the genome of S. oneidensis MR-1 encodes all three types of NDHs (type I: SO_1009-SO_1021, type II: SO_1103-SO_1108 and SO_0902-SO_0907, type III: SO_3517). It has been reported that electrical current production or growth rate with NAG as the substrate under anaerobic conditions would be extremely inhibited by knocking out NDHs (Kasai et al., 2019; Madsen and Teravest, 2019). However, the simulation results of iLJ1162 indicated that there was no obvious decrease of electron production or growth rate. Further analysis revealed that there are two alternative pathways with NDHs function, including 1) LdhA (pyruvate reductase encoded by SO_0968, EC: 1.1.1.28) and Dld (lactate dehydrogenase by SO_1521, EC: 1.1.2.4) (Figure 8B); 2) ProC (pyrroline-5-carboxylate reductase encoded by SO_3354, EC: 1.5.1.2) and PutA (proline dehydrogenase encoded by SO_3774, EC: 1.5.5.2) (Figure 8C). The former has been experimentally verified (Kasai et al., 2019), while the latter was uncovered by iLJ1162 for the first time. However, the two pathways failed to play a role in electricity production or cell growth, probably due to the strict transcriptional regulation and low pools of intermediate metabolites (d-lactate or l-proline) (Kasai et al., 2019; Madsen and Teravest, 2019). Thus, overexpression of NDHs and/or activation of LdhA-Dld/ProC-PutA might be an effective strategy for channeling more electrons toward MQ to improve EET rate.

In MES systems, electrons from the cathode are transferred to the microbial cells and used for the synthesis of value-added chemicals (Logan et al., 2019). The MES process of S. oneidensis MR-1 can be summarized as NADH generation by receiving extracellular electrons and protons using NAD⁺ (Yang et al., 2015; Rowe et al., 2017). Thus, the reaction of “NADH_R” (2.0 e⁻ + H_p + NAD_c → NADH_c) and ‘Sink_e⁻’ (→ e⁻) were added to represent the electron receiving process. Six high value products were chosen for simulation in MESs, including 3-hydroxypropanoate (3HP), isopropanol (2ppoh), 3-hydroxybutanoate (3HB), acetoin (Actn), succinate (Suc), and 2,3-butanediol (23btdl). All these products require no more than two heterologous reactions in MR-1. The cell growth was omitted to simplify the complexity of biosynthesis pathways for different products, which is also consistent with that EAB almost did not grow in MESs (Chen et al., 2018; Tefft and Teravest, 2019) The demand reactions of these products are set as objective function respectively.

The simulation results of these six products in MESs using l-lactate as sole carbon source under anaerobic conditions are shown in Table 2. Without inward electron (no extra electron supply), the simulated yields were lower than their theoretical maximum yields calculated from reducing degrees. However, with inward electrons, all products achieved or exceeded the theoretical maximum yields calculated from the reduction degrees. In the simulation of MESs, CO₂ is infinitely supplied, which is consistent with the actual MESs aerated with CO₂ to ensure a strict anaerobic condition (Chen et al., 2018). Thus, with extra electron supply (inward electrons), it is theoretically feasible to fix environmental CO₂, and consequently break mass balance merely based on substrate. As one of the C4 platform chemicals, 3-hydroxybutanoate (3HB) was chosen as an example for MESs simulation using iLJ1162. As shown in Figure 9, a simulated yield of 0.43 mol/mol was obtained using l-lactate as the sole electron source, which is only 64% of the theoretical maximum yield calculated from reducing degrees. This can be explained by the fact that the key intermediate metabolite formate cannot be further converted to 3HB by the metabolic network of MR-1, resulting in the loss of one third of total carbon in l-lactate. Moreover, to obtain NADPH for 3HB synthesis, a part of carbon fluxes was directed into the serine bypass and PP pathway, and oxidized to release electrons for NADPH generation.

With inward electrons from an electrode, the simulated yield was significantly increased to 1 mol/mol l-lactate, which was 57% higher than the its theoretical maximum yield. This can be explained by the significantly increased fluxes of the threonine bypass (Figure 9) (Lin et al., 2015). The threonine bypass is able to not only utilize additional reducing equivalents to fix extracellular CO₂ (4.06 mmol/gDW/h, Figure 9), but also re-assimilate formate released by PFL (4.06 mmol/gDW/h, Figure 9), and generate the key precursor acetyl-CoA with a higher yield than that obtained merely via the carbon-losing PFL pathway. As a result, 3HB can be generated at a rate of 4.06 mmol/gDW/h in silico, which represents the yield of 1 mol/mol l-lactate. It is worth mentioning that the 3HB yield would not be infinitely improved in MESs of MR-1, which might be owing to the lacking of key intermediate metabolite formate. Finally, these results indicated that metabolic regulation of the threonine bypass will be an important aspect of MR-1 for enhancing the yield of chemicals in future MES.

Since glycerol turns out to be the best substrate for election generation (Figure 6), we also simulated the product synthesis using glycerol in MESs (Supplementary Datasheet S4, Supplementary Table S10). Three major findings were obtained: 1) Without inward electrons, the product yields using glycerol are higher that using l-lactate, which is consistent with that in electron generation. 2) With inward electrons, the product yields using glycerol or l-lactate are equivalent, which means that extra electron supplying might eliminate the difference caused by different reducing degree in substrates. 3) The requirement of inward electrons for the same product is significantly lowered when using glycerol rather than l-lactate. This indicates that glycerol might still be a better substrate in MESs, which lowers the demand for the electron uptake rate of S. oneidensis strains.