Model

×

Section 1: Wet Lab

To further optimize perchlorate reduction efficiency in E. coli and enhance metabolite exchange within our closed-loop, two-microbe co-culture system, Team Cornell employed two computational approaches to guide our experimental conditions and predictive analysis:

  1. Michaelis-Menten Kinetics: This approach was used to model the perchlorate reductase reduction dynamics in the transformed E. coli by analyzing its enzyme-substrate interactions and time-course concentrations across the pcrABCD and cld biological pathways. This can be used to determine specific rate-limiting steps and identify how to optimize specific enzyme and substrate concentrations for experimental design.
  2. Michaelis-Menten
    Equation 1.1. Michaelis-Menten Equation
    MM Enzyme Kinetics
    Equation 1.2. Michaelis-Menten Enzyme Kinetics
  3. Co-Culture Logic Gate: Using an AND-gate framework and a series of Boolean logic operations, we can model the interaction between E. coli and Synechococcus PCC 7002 to better understand how each microbe responds to different environmental inputs and contributes to the overall metabolite cycle in the co-culture. In this system, E. coli acts as the upstream processor by reducing perchlorate into chloride, while Synechococcus functions as the downstream strain that uptakes the chloride and produces sucrose to support E. coli growth.
  4. A. Upstream Process - E. coli Perchlorate Reductase (Pcr) and Chlorite dismutase (Cld)

    Using the AND Gate framework, we can treat E. coli as the upstream processor: PcrABCD (perchlorate reductase complex) reduces perchlorate to chlorite, and Cld (chlorite dismutase) converts chlorite to CO2 and Cl-. The two enzymes must act simultaneously to complex the detoxification process and feed into the Synechococcus system. Using a standard Hill's coefficient of n = 2, we can standardize the inducers to model similarly to a moderately switch-like behavior that allows for expression to rise sharply around a specific threshold dissociated constant (Kd) for each enzyme while minimizing the effects of leaky basal expression. As seen in Figure 1.1, this strategy also allows us to ensure each individual enzyme output is symmetric and co-limiting.

    AND Gate Behavior
    Figure 1.1. Figure 1.1. AND-gate behavior of Pcr and Cld module. Final Detox output (nM) at steady state: Both Off ≈ 0 nM, pcrABCD, pcrABD Only ≈ 20 nM (Half-response), Cld Only ≈ 20 nM (Half-response), Both On ≈ 40 nM.
    Detoxification Output
    Figure 1.2. Figure 1.2. Steady-state detoxification output heatmap. With the blue stripes/areas representing zero input confirming AND dependence; rapid transition at red plateau regions to indicate sharp activation response near Kd and max flux when box inducer are at a high concentration
    3D AND Gate
    Figure 1.3. Figure 1.3. The 3D surface plot provides another visualization of the detoxification output at varying inducer concentrations and enzyme activation. With the “steep wall" area indicating that the gate is thresholded and saturating: sharp activation around Kd, then near-max flux across most of the interior. The symmetry indicates Pcr and Cld are co-limiting and balanced.

    At very low Kd (~0.1 nM), the promoters for pcrABCD and cld are very sensitive to their inducers. Even the tiniest amount of inducer can cause strong enzyme expression, resulting in high detoxification (accumulation of CO2 and Cl-) output from the perchlorate reductase pathway. As a result, the system behaves more like an OR gate since either one enzyme complex can drive a substantial response. However, at the intermediate Kd range (~3-10 nM), both the pcrABCD and cld must reach a sufficient threshold of induction before the system reaches the maximal detoxification output, allowing for the most distinct AND-like response when both inducers are present. At very high Kd (>30 nM), the enzymatic system becomes less responsive, as a high concentration (strong induction) is needed from both inputs to achieve activation; however, the actual detoxification output decreases as the enzyme expression saturates out under maximal induction.

    Sensitivity Analysis
    Figure 1.4. Sensitively Analysis of varying Kd values to evaluate detoxification output of Cl- and CO2

    Both the detoxification output of the upstream E. coli system feeds into two time-resolved fluxes that become the sources of input for the Synechococcus module. With the chloride production rate from the chlorite dismutase feeding into the chlorite capacity of the Synechococcus module and the CO2 production rate from the perchlorate reductase feeding into the available inorganic carbon (DIC) rate. With these two inputs, we're able to examine the Synechococcus organic carbon fixation rate via its available biomass and also sucrose formation.

    B. Downstream Process - Halorhodopsin (Chloride Uptake) and Sucrose Fixation

    As E. coli reduces perchlorate and releases CO2 and Cl- into the co culture medium, Synechococcus takes up those outputs and converts the inorganic carbon into fixed carbon such as sucrose. In this closed loop, Cl⁻ acts as a cofactor or quota that can throttle photosynthesis if undersupplied, while DIC (dissolved inorganic carbon = CO₂ + HCO₃⁻) is the actual substrate for carbon fixation. Synechococcus can import chloride via halorhodopsin (HR), a retinal binding microbial rhodopsin originally described in haloarchaea and expressed here heterologously, which uses light energy to pump Cl⁻ inward through a photoisomerization cycle that moves chloride from outside to inside the cell.

    For our co culture, expressing HR in Synechococcus is a practical way to boost effective chloride uptake under light, reducing the chance that a chloride quota limits photosynthetic carbon fixation, especially when external chloride is low. If we do not model DIC explicitly, Synechococcus would keep fixing carbon even after the medium runs out of CO2 or HCO3-, which is not realistic. The equations used below are:

    1. fL(I) = (I) /( KI + I), describes the light limitations on photosynthesis
      1. I = light intensity
      2. KI = half-saturation constant
    2. vfixation =Vfixation,max fL(I) * (DIC)/ (KCi + DIC), describes the max carbon fixation scaled by light and inorganic carbon availability
      1. Vfixation, max = max fixation rate per biomass (mmol)
      2. DIC = dissolved inorganic carbon (CO2 + HCO3-)
      3. KCi = half saturation constant for DIC
    3. μ0=Yx/c​vfix−msyn, growth rate from fixed carbon minus maintenance
      1. Yx/c = biomass yield from fixed carbon
      2. msyn = maintenance drain
    4. Chloride Need
      1. Demand of Chloride: rcl, required = qcl μ0
        1. qcl = Cl needed per unit biomass growth
        2. rcl, required = Demand of chloride needed
      2. Capacity of Chloride: vCl, capacity = VCl, max*(Clmedium)/ (KCi+Clmedium)
        1. VCl, max = max uptake per biomass
        2. Clmedium = concentration of chloride in medium
    5. Actual Chloride Uptake: vCl = min(vCl, capacity , rcl, required)
    6. Penalty when CI is low:
      1. Min (1 , vCl, / (rcl, required)
    7. Mass balances (bioreactor)
      1. dClmedium/dt = Jc1E.Coli - vCl Xsynechoccocus
        1. Jc1E.Coli = Cl produced by E.coli in Cld step
        2. Xsynechoccocus = Synechoccocus Biomass Concentration
      2. DIC Concentration and Pool
        1. dDIC/dt =JCO2in - ((vfixation,effXsynechoccocus )/ (YC, fixation)
          1. JCO2in = CO2 supplied from E.Coli
          2. YC, fixation) = mmol of DIC per mmol C fixed)
      3. dXsynechoccocus /dt = μXsynechoccocus → Biomass of Synechococcus
      4. dsucrose/dt = fsucrose vfixation,eff Xsynechoccocus - ksucrose[Sucrose]
        1. fsucrose = fraction of fixed C secreted as sucrose
        2. ksucrose = consumption/decay rate

    The graphs below show the Chloride concentration, DIC concentration, Sucrose concentration and Synechococcus biomass over the course of 10 hours.

    Synechococcus
    Figure 2. Synechococcus over ten hours

    The co-culture is carbon-limited, not chloride-limited. In the Synechococcus culture, DIC fell from ~1 mM to ~0 within ~20-40 min, so biomass and sucrose became pinned to the CO₂ feed. On the E. coli side, the AND-gate heatmaps show detox output peaking only when both inducers are ON at intermediate Kd (≈3-10 nM); low Kd leaks (OR-like), high Kd gives weak absolute output.

  5. Flux Balance Analysis (FBA): A constraint-based modeling approach used to simulate the metabolic exchange within the co-culture system as a complex network, allowing us to predict steady-state flux distributions based on mass balance and defined reaction constraints. This method helps us optimize key metabolic conditions—such as perchlorate turnover rates, chloride uptake efficiency, and sucrose fixation—and identify potential pathway knockouts that could enhance the production of critical enzymes involved in the system.

Enzyme Kinetics of Differing Perchlorate Reductase

Michaelis-Menten Assumptions:

  1. Steady-state approximation: In our perchlorate reduction model, we assume that the rate of PcrABCD-substrate complex formation (enzyme binding to perchlorate or chlorate) quickly reaches a steady state where complex formation and breakdown occur at equivalent rates. This allows us to approximate the catalytic rate as constant over time.
  2. Free ligand assumption:We assume the substrate concentration in this case perchlorate and its intermediates products (chlorate and chlorite) is significantly higher than the total enzyme concentration (PcrABCD or Cld). This allows us to assume the reaction rate depends primarily on the amount of substrate concentration that we have rather than the enzyme concentration.
  3. Rapid equilibrium assumption:We assume that the enzyme-substrate binding and catalytic conversion steps proceed primarily in the forward direction (toward perchlorate reduction and oxygen/chloride generation), and that product rebinding is negligible.

Co-culture Conditions and Metabolite Linking between E. coli and Synnechococcus PCC. 7002

Modeling bacterial consortia using logic computing

The cooperative co-culture system between E. coli and Synechococcus PCC 7002 can be further effectively modeled using logic computing frameworks designed for bacterial consortia. Specifically, the perchlorate reduction mechanism performed by the E. coli and the chloride uptake and carbon fixation performed by the Synechococcus can be modeled using a XOR gate design. As described in Figure 1 , the XOR gate employs Boolean logic to model the mechanism of the two model organisms in which each strain performs a complementary function towards an unified biological function.

AND Gate
Figure 3. AND logic gate only produces true (1) output response when every input is true.

This formalized logic gate framework allows our team to simplify the complex mechanism into simpler biologically-implementable units (AND, OR and NOT Gates) and ordinary differential equations (ODEs) that can be used to model the two microbial cooperative behaviors across a series of assumptions and mathematical interpretations. In our co-culture system, the E. coli's perchlorate reductase complex protein (pcrABCD) and chlorite demutase (cld) will serve as the input responses while the Synechococcus acts as more of a downstream processing microbe that consumes the chloride and recycles the fixed carbon back to E.coli in the form of sucrose. This system design seeks to emulate the quorum-sensing mechanism exhibited by many microbial species in the environment. In the case of a co-culture system, both microbial species must exhibit quorum sensing and symbiotic behaviors to accurately stimulate the closed microbial metabolic system.

A series of ODE equations are used to simulate the gene expression dynamics, metabolite exchange, and system-wide response underlying varying environmental conditions. These parameters include diffusion constants to model the exchange of chloride and sucrose between the two microbes, degradation rates to simulate the consumption of metabolite, and both species' growth dynamic and Michaelis-Menten constants to predict the enzyme-substrate catalytic reactions in the E. coli.

This approach allows us to minimize the wiring complexity in the biological system while also creating a predictive model that not only provides us with a mechanistic understanding of the complex molecular and metabolic dynamics, but also guides the optimization of experimental conditions. Through stimulating time-course behavior of essential metabolites (e.g. perchlorate, chloride, sucrose), this approach allows us to better anticipate gene expression level, enzymatic behaviors and metabolite exchange rates that are key in developing a scalable application such as Martian regolith remediation on a greater scale.

Structural Modeling of pcrAB Construct

Perchlorate reductase (PcrABCD) is the periplasmic enzyme complex that reduces perchlorate (ClO4-) to chlorite (ClO2-), before which chlorite dismutase (Cld) catalyzes the conversion of chlorite to Cl- and O2 [1]. Genes encoding this metabolism are located on the Perchlorate Reduction Genomic Island (PRI), which is horizontally transferred and tightly conserved in strains like Azospira and Dechloromonas aromatica. Within the perchlorate reduction pathway, pcrA encodes the catalytic a-subunit, electron transfer beta-subunit within pcrB. pcrC and pcrD are required coenzymes: pcrC a periplasmic tetraheme c-type cytochrome which is responsible for electron transfer from the membrane quinone pool to pcrB while pcrD is a necessary coenzyme for the assembly of the pcrAB complex.

As such, we chose to model the pcrAB complex given its catalytic function with Michaelis-Menten kinetics. PcrAB both act in concert during perchlorate reduction, further motivating the decision to model them together. To better study the relationship between pcrA and pcrB, a point mutation was modelled at residue number 461. This point mutation replaces tryptophan (W) with glutamic acid (E), which functionally deactivates pcrA.

Wt vs mutation
Figure 4. Range of enzyme kinetics for the wild type Azospira (blue) compared to Azospira (functionally inactive PcrA) (red).

After plotting the Michaelis-Menten curves, a significant drop in enzyme affinity and rate occurs at lower concentrations in the Azospira strain with a mutation. Despite a higher apparent kcat at saturating substrate, the mutated Azospira has a ~26x drop in catalytic efficiency (kcat/Km) which indicates severely reduced performance at environmental ClO4- (Figure 4). This suggests impaired long-range electron transfer or substrate gating when the PcrAB wiring is disrupted. After plotting the Michaelis-Menten curves, the loss in catalytic efficiency observed in the mutated Azospira strain indicates that disruption of the PcrAB complex severely impairs substrate binding and electron transfer. Thus, the W461E mutation found here can serve as a negative control, and establishes that both PcrA and B are required for effective perchlorate reduction.

To further elucidate the structural basis for this relationship between PcrA and PcrB, the structure of the perchlorate reductase complex was modelled to visualize how the catalytic (α) and electron-transfer (β) subunits cooperate to facilitate Mo-centered reduction.

Normal pcr
Figure 5. Normal view of perchlorate reductase.
Mutated pcr
Figure 6. Mutated view: Effect of a single point mutation near molybdenum

A key mutation is the effect of changing the residue #461 from tryptophan to glutamic acid [1], visualized in pink on the above figure. This mutation occurs near the molybdenum (Mo) cofactor pocket and forms part of an aromatic gate that stabilizes substrate entry and controls product retention. The negative charge of the glutamic acid repels electrons, which prevents electrons from being coordinated by the molybdenum atom during perchlorate reduction. Perchlorate Reductase is part of the Type II DMSO Reductase Protein family, which is part of the larger DMSO Reductase Protein Family that includes proteins like nitrate reductase. All proteins of this family Type II DMSO Reductase Protein Family require molybdenum (Mo) that is buried in the catalytic (α) subunit inside a bis-molybdopterin guanine dinucleotide (bis-MGD) cofactor pocket. Initially, Mo has an oxidation state of +4 with 2 bonds to sulfur, 1 double bond to oxygen and 1 amino acid side chain bond (typically Asp in Type II DMSO Reductase Protein Family), Mo is able to strip one oxygen from the substrate (in this case perchlorate) to produce a form of oxidation state +6 with all of its previous coordinate bonds plus an additional double bond to oxygen. The Mo (VI) is able to reversibly form water with 2 additional electrons from the Fe-S centers and 2 protons. In the presence of atmospheric oxygen, Mo will double with 3 oxygens instead of the normal 2. This “kills” the enzyme in a sense as due to oxygen's high electronegativity it is able to pull the electron density away from the Mo, reducing its effective energy state and reactivity (Figure 6). The Mo becomes electron poor and unlikely to further react in redox reactions where its role is to act as a temporary oxygen holder before being moved elsewhere.

VEP
Figure 7. Vacuum electrostatic potential map of the PcrAB catalytic pocket surrounding the molybdenum (Mo) cofactor. Mutated above, normal below.

To further visualize the effect of the point mutation at residue 461, the electrostatic surface potential (ESP) was calculated using the Adaptive Poisson-Boltzmann Solver (APBS) under vacuum conditions, meaning the surrounding environment was treated as empty space without solvent or ionic effects.

In the unmutated PcrAB complex, the MO pocket exhibits a predominantly neutral to slightly positive potential (blue) surrounding the substrate entry channel, while the deeper regions near the cofactor and coordinating aspartate residue show localized negative potential (red) arising from oxygen and carboxylate groups. These charge distributions collectively form a balanced electrostatic landscape that is favorable for stabilizing the incoming perchlorate (ClO4-) substrate.

Following structural and electrostatic modeling, wet lab experiments were designed first to express pcrAB alongside necessary cofactors pcrC and D. Following this, it would be possible to validate the proposed role of pcrA and pcrB in catalysis and electron transfer within the perchlorate reductase complex. The overarching aim of these experiments is to experimentally confirm that both subunits are necessary for efficient perchlorate reduction, as suggested by our structural data and kinetic modeling. Another perturbation to the system is the effect of oxygen, as it competes for electrons that cannot be shuttled to perchlorate and potentially damages sensitive perchlorate reducing enzymes.

Effect of O2 on Perchlorate Reduction via Stochastic Simulation

Biological redox systems, especially those operating under microaerobic conditions, often exhibit fluctuations that cannot be captured by purely deterministic models. The perchlorate reductase (PcrAB) complex operates in a periplasmic environment where stochastic changes in oxygen concentration and electron donor availability can significantly affect catalytic turnover. While our deterministic Michaelis-Menten analysis describes the mean behavior of PcrAB, it does not capture the inherent randomness of molecular interactions or the probabilistic nature of Mo redox cycling under O₂ stress [3]. Since cld catalyzes the reaction of chlorate to oxygen and chlorides, we hope to demonstrate a negative feedback loop in effect of perchlorate reduction is observed.

To further investigate how oxygen exposure influences PcrAB catalysis, we developed a stochastic simulation of the perchlorate reduction reaction network. This approach accounts for biological noise and simulates discrete reaction events over time rather than relying on continuous rate laws and structural analysis.

To account for stochasticity, we reformulated our kinetic rate equations into propensity functions and implemented the Gillespie algorithm, or stochastic simulation algorithm (SSA). This method captures the probabilistic evolution of enzyme-substrate interactions in time, allowing us to explore transient behaviors such as partial inhibition, random reactivation, and oxygen-dependent quenching of the Mo cofactor.

In this model, we fit the effect of O2 to gradually reduce the propensity for the catalyzation of perchlorate (ClO4-). We incorporated the effect of the inhibitory effect on perchlorate reduction, which was clearly demonstrated in our simulation. As O2 levels increased, perchlorate levels began to build up, clearly demonstrating the negative feedback from O2 production catalyzed by cld.

Stochastic
Figure 8. Stochastic Simulation Modelling Perchlorate Buildup Depending on O2 In-Flow from Low to high O2.

Conclusion

Structural Analysis

Reinforced Wet Lab's decision to study the expression and enzyme kinetics of pcrAB, which requires necessary coenzymes pcrC for electron transfer from the quinone pool and pcrD which allows for assembly.

Effect of O2 on Perchlorate Reductase

References

[1] Ji, Weiyue, et al. “A Formalized Design Process for Bacterial Consortia That Perform Logic Computing.” PLoS ONE, vol. 8, no. 2, 28 Feb. 2013, pp. e57482-e57482, https://doi.org/10.1371/journal.pone.0057482. Accessed 30 Jan. 2025.

[2] Melnyk, Ryan A., and John D. Coates. “The Perchlorate Reduction Genomic Island: Mechanisms and Pathways of Evolution by Horizontal Gene Transfer.” BMC Genomics, vol. 16, no. 1, 26 Oct. 2015, https://doi.org/10.1186/s12864-015-2011-5.

[3] https://pubs.acs.org/doi/10.1021/acscatal.7b01749

Section 2: Product Development Modeling

Thermal Modeling and Material Suitability under Martian Environment

The Product Development modelling is divided into two major components.

The first focuses on future scale-up planning for Martian missions. Although the ARC bioreactor was originally designed based on current laboratory-scale requirements and the logistical constraints of Mars transportation, we must anticipate the potential bioreactor demands that ARC could face in future applications. A key consideration is how heat exchange between the bioreactor and its surrounding environment may affect overall performance. On Mars, this factor becomes particularly critical, as it directly determines the choice of construction materials, modifications to the operational workflow, and ultimately, how these findings can be integrated into large-scale use cases. In our case, large-scale Martian soil treatment.

During the scale-up analysis, several key variables were identified as determinants of bioreactor performance:

  1. External placement feasibility:Given that human habitat space on Mars will be a highly limited resource, any scaled-up deployment may require locating bioreactors outside pressurized habitats. This introduces direct dependence on the convective heat transfer coefficient, which characterizes how efficiently heat is exchanged between the reactor surface and the external atmosphere.
  2. Target reactor volume:The required culture volume is closely tied to the volume of soil to be processed. As there are currently no empirical data from interplanetary missions to define realistic soil throughput demands, we conservatively estimate a reactor capacity of 10 L—sufficient for effective perchlorate degradation and soil remediation.
  3. Material selection:The chosen material must satisfy all preceding constraints—thermal compatibility with Martian environmental conditions, feasible mass and density for transport, and long-term structural stability. For these reasons, our current design validation is conducted at a laboratory scale, which substantially reduces prototyping cost, launch payload, and total transport volume while enabling meaningful thermal and material characterization for eventual extrapolation to Martian conditions.

To quantitatively evaluate how these environmental and material constraints influence reactor performance under Martian conditions, we introduced Biot number (Bi) modelling as a thermal characterization framework.

The Biot number, defined as Bi = hL/k, represents the ratio between external convective heat transfer (governed by the Martian atmosphere, with a low convective coefficient h) and internal conductive heat transfer (determined by the material's thermal conductivity k and reactor dimension L).

By computing Bi values across different materials (borosilicate glass, stainless steel, AL-6XN, carbon steel, and aluminum alloy) and Martian convective conditions (0.6-9 W m-2 K-1), we can identify which reactor configurations maintain sufficient thermal uniformity. Convective heat transfer coefficients under Martian surface conditions were taken from Soria-Salinas et al. (2015), Convective Heat Transfer Measurements at the Martian Surface, which reports h ≈ 0.6-9 W m-2 K-1 for wind velocities of 0.5-4 m/s at T ≈ 250 K and P ≈ 720 Pa. This modelling serves as a quantitative link between the environmental assumptions of future Martian operations and the practical material and design decisions for the ARC bioreactor family.

These values span more than two orders of magnitude, allowing us to compare materials ranging from highly insulative (glass) to highly conductive (aluminum). In this context, lower k values indicate poorer heat conduction and thus higher temperature gradients within the reactor, while higher k values promote rapid heat transfer and a more uniform internal temperature distribution. This wide range provides a clear physical contrast when analyzing how different reactor materials perform under Martian convective heat transfer conditions.

To further investigate the scalability of the ARC bioreactor system for future Martian soil remediation missions, the modelling framework was extended from the baseline 1 L laboratory configuration to a 10 L conceptual design. This up-scaling reflects a realistic operational scale for in-situ perchlorate degradation, where larger culture volumes would be required to achieve meaningful soil treatment throughput under mission-level constraints.

In this scenario, the reactor's characteristic length (L)was proportionally increased from 0.10-0.17 m (for 1 L) to 0.35-0.50 m, representing the expected geometric range for a 10 L vessel of similar aspect ratio. An increase in L therefore yields a proportional increase in Bi, implying greater potential for internal temperature gradients at larger volumes. This analysis enables quantitative assessment of whether each material can maintain Bi < 0.1, the threshold for near-isothermal performance, even at a 10 L scale. Through this comparative framework, we can determine which structural materials remain thermally viable for scaled-up ARC systems deployed in Martian environmental conditions.

To translate the conceptual thermal framework into quantitative results, we developed a numerical simulation of Biot number variation under Martian environmental conditions. By iterating across a range of realistic h, k, and L values, we can visualize how different materials behave thermally when scaled to both laboratory and mission-scale reactors.

borosilicate
Figure 9.1. Biot Number (Borosilicate Glass)
SS304
Figure 9.2. Biot Number (Stainless Steel)
stainless alloy
Figure 9.3. Biot Number (Stainless Alloy)
carbon steel
Figure 9.4. Biot Number (Carbon Steel)
anodized aluminum
Figure 9.5. Biot Number (Anodized Aluminum)

The simulation was implemented in MATLAB, using parameter sweeps to generate a 3D surface plot of Biot number as a function of reactor height and material index (h/k).

This approach allows us to identify which material-geometry combinations achieve Bi < 0.1, the threshold for near-isothermal operation, and which configurations fall into higher-Bi regimes where temperature gradients may arise.

To further quantify this relationship, the model also allows the user to input a target Biot number (e.g., Bi = 0.05), enabling reverse identification of the corresponding h, k, and L values that achieve this thermal equilibrium condition.

This workflow bridges the conceptual understanding of Martian convection with practical design parameters for the ARC bioreactor family.

The five plots respectively illustrate how the Biot number varies across different materials under Martian environmental conditions (h = 0.6-9 W m-2 K-1). The x-axis represents the reactor height L (0.10-0.17 m), the y-axis denotes the convective heat-transfer coefficient h, and the z-axis shows the calculated Biot number (Bi = h x L / k). The color gradient from blue to yellow indicates increasing Biot values. A smaller Biot number corresponds to a more uniform internal temperature distribution and higher thermal-transfer efficiency. In general, Bi > 0.1 implies a pronounced internal temperature gradient, whereas Bi < 0.1 indicates an effective isothermal system.

For our bioreactor application, maintaining an isothermal chamber is advantageous—it ensures that all regions of the culture experience consistent reaction conditions. This thermal stability is particularly beneficial for maintaining co-culture balance and sustaining biofilm formation over long-term operation.

Overall, borosilicate glass exhibits the lowest thermal conductivity among the tested materials. While this provides strong thermal insulation, it also limits internal heat diffusion, potentially creating localized hot or cold zones. However, at the laboratory scale of the ARC reactor (approximately 1 L), such local gradients are minimal, justifying our continued use of a glass beaker as the primary culture vessel. The other metallic materials demonstrate significantly higher conductivities, enabling more even temperature profiles—an important feature for larger-scale systems or precise temperature control.

Among all tested materials, anodized aluminum alloy (k ≈ 170 W m-1 K-1) stands out. Its exceptionally high conductivity results in nearly uniform temperature distribution (Bi < 0.01) even under low Martian convection coefficients. This characteristic makes it ideal for off-Earth applications, where external environmental fluctuations could otherwise destabilize internal conditions.

In the ARC reactor design, these modeling insights directly informed our material selection. Stainless steel was initially considered but ruled out due to machining complexity and higher mass. In contrast, aluminum alloy offered excellent manufacturability via CNC milling, low weight, and superior heat-transfer capability. Consequently, we selected CNC-machined anodized aluminum for the top structure—a choice that aligns with both our engineering modeling results and the practical requirements of a lightweight, thermally stable Martian bioreactor.


Modeling Population of Cells in Media

Cell Concentration

The growth of cells in solution can be modeled with a whole host of equations: the Monod, Tessier, Contois, and Moser equations, among others. For the purposes of simplicity, a modified Monod equation will be utilized that allows for substrate inhibition of growth at higher concentrations.

Equation 1
Equation 2.1. The first factor is the modification. The second factor is the general Monod equation. Cp is the product concentration, rg is the growth rate, Cp* is the product concentration when cellular metabolism goes to zero, μmax is the maximum specific growth reaction rate, Cs is the substrate concentration, Ks is the Monod constant, and n is a constant that can be found through experimentation.

The death of cells can be modeled by the following equation.

Equation 1
Equation 2.2. Equation modeling cell death rate rd. Ct is the concentration of a toxic material, and kd and kt are both constants for the death rate and toxic material buildup respectively. Ctkt can be zero if the term is irrelevant to the system of interest.

Substrate and Product Concentration

Since the Monod equation takes into account the product and substrate concentration as well, two other mass balances for the substrate and product should and can be determined in a similar manner to the cell mass balance. The equations are such:

Equation 2.5
Equation 2.5. Equation modeling the change in substrate concentration over time in the growth phase. Ys/c is mass of substrate consumed for growth/mass of new cells, rsm is the substrate consumption rate for the maintenance of the cells (no growth). The constant ms/w is the ratio of mass of substrate consumed for maintenance divided by the dry weight of the cells and time.
Equation 2.6
Equation 2.6 Equation modeling the change in substrate concentration over time in the stationary phase.
Equation 2.7
Equation 2.7. Equation modeling the change in product concentration over time in the growth phase. Yp/c is the mass of a product divided by the mass of cells formed in the growth phase. The constant rp is the rate of production formation.
Equation 2.8
Equation 2.8. Equation modeling the change in product concentration over time in the stationary phase. Yp/s is the mass of a product divided by the mass of substrate needed by cells to form the product.

The above equations can be combined to model the growth of a bacterium at any point in its growth or stationary phase.

Modeling Accumulation of Biomass on Carrier Lattice

The following is from a paper authored by Picioreanu et al. 1997 in which they aimed to model the formation of biomass on a gel bead. The table was included to demonstrate all the various parameters that were taken into account in growing the biomass. Many of these parameters are of relevance to the first part of the task: modeling the growth of bacteria in the solution.

Modeling the growth of biofilm on a lattice, a more geometrically complex structure than a gel bead, will no doubt have more factors associated with it. For purposes of practicality, using a complex mathematical model may not be in Cornell iGEM's best interest. While the growth of bacteria in the solution is easier to modulate, the growth of biofilm on the carrier may lead to more variable results depending on the carrier's structure.

Table 2.1
Table 2.1. Table credits go to Picioreanu et al. 1997.

Here is where the creative element of the task comes into play. We can map the values of bacteria concentration in the solution to the numerical value of the cell count that adheres to the biofilm carrier over a given period. This is not a beautiful, mathematical solution, but what this calibration will do is allow the experimenter to optimize biofilm formation precisely by simply changing the cell concentration in the solution. We can assume, for the purposes of experimentation, that they are related by a constant and a function that is related to the ratio of product concentration and substrate concentration:

Equation 2.9
Equation 2.9.

The fundamental guess here is that the mapping between cell count of the biofilm and the concentration of bacteria in the solution are related by a function whose input is the ratio between the product concentration and the remaining substrate concentration.

Computation of Biofilm Formation

Using the equations highlighted above, the mapping function is as follows:

Equation 2.10
Equation 2.10.

Where Nbf is the attached biofilm (cells) per unit area, λ is the attachment coefficient, and f(Cp/Cs) is the heuristic function representing the product to substrate ratio affecting the attachment. If the following model had feedback inhibition, the heuristic function could be defined as follows: Equation 2.11

Equation 2.11.

This model shows the decrease in biofilm formation as there is an accumulation of biomass. The parameter n represents the Hill Coefficient. As biofilm attachment is non-linear, this is needed to control how sharply biofilm formation decreases as the ratio of product to substrate increases. If n is closer to 1, this means that the function follows Michaelis-Menten inhibition, if it is larger it follows more of a sinusoidal behavior, and if it is less it is less sensitive to biofilm build up. This coefficient can be found from experimental data of controlled substrate and product concentrations over time. To simulate this model, MATLAB will be used.

Limitations

After a series of multiple tests, it was determined that the rate of formation of the biofilm is quite a complex matter, and it is difficult to quantify biofilm formation. A mathematical model that can model the rate of biofilm growth would certainly be beautiful but may be inconsistent. The Product Development subteam members will experiment with open minds and may look for other equations that better model biofilm formation as more experimental evidence comes to light.

Data Collection and Optimization of the Reactor Chambers

A methodology for finding the efficiency of each of the modules was defined in the Hardware section. Thus, depending on the experimental results, it is possible to find the power output levels where the ratio of valuable output to energy input is maximized. Operating in these conditions is ideal. As of now, testing has not begun, but once results are obtained, an outline of optimal operating conditions can be created. Mass and thermal balances can also be quantified, but again, without data, no concrete conclusions can be made. This is the next large priority for the PD subteam.

Ansys Fluent

Ansys Fluent is a computational software that yields insights into the nature of fluid flow in various systems. This includes shear forces, flow vectors, contours, etc. The Product Development subteam is currently working on modeling the fluid flow in our bioreactor using Ansys Fluent as a way to find inefficiencies in the system and to better understand the various fail-points in the setup.

Heat Transfer

The carrier used in the design is glass wool, which is an insulator. Glass wool works well as a carrier as it has high surface area which allows for easy attachment, but it has low effective thermal conductivity, therefore it impedes heat transfer. Therefore, maintaining 37°C could be a challenge and could have temperature gradients.

To check whether the wool would be isothermal or if there would be internal gradients, the Biot number serves as an estimate.

The Biot number, NBi, is a heat transfer calculation that is commonly used in heat transfer calculations. This number is the ratio of internal conductive resistance and external convective resistance and it is dimensionless.

The equation below is used to calculate the heat Biot number.

The convective heat-transfer coefficient between the fluid and glass wool in this case would be in the format of assisted convection, as the convection occurs by the heat pad. The convective heat transfer coefficient varies depending on the fluid, and in forced convection of water and liquids the range of this value is 10 - 1000 Wᐧm-2ᐧk-1 and it is assumed that this value is 200 Wᐧm-2ᐧk-1. The length of the glass wool will be experimentally determined, depending on the mass, volume, and length of the wool.

The thermal conductivity (keff) of typical glass wool is 0.0343 Wᐧm-2ᐧk-1.

If this number is above 1, the convection heat rate dominates conduction heat transfer. The goal is to have a consistent temperature within the reactor, and therefore this is good. Thus, a conclusion is that it will take longer for the wool to reach the temperature of 37°C as compared to a system fully composed of fluids. Some solutions to this challenge would be to make the glass wool thinner, as this would reduce the characteristic length of the glass wool.

Mission to Mars...