Wet Lab Modelling

Executive Summary

Oncurex is an industrially scalable platform to increase the supply of ursolic acid for breast cancer treatment. The wet lab subteam planned to use Saccharomyces cerevisiae for all experimentation. To direct our approach, Team Cornell used two computational models to inform our wet lab experimentation:

1. Michaelis-Menten Kinetics: Analyzing substrate-level interactions to pinpoint methods of enzymatic intervention for wet lab optimization
2. Graph-Theory Based Network Topology Analysis: Considering global metabolic perspectives to identify hubs, bottlenecks, flux imbalances, and system robustness for optimization of ursolic acid production

Our modelling approaches were designed to accomplish the following:

Our approaches revealed that the MdOSC1 variant of alpha amyrin synthase, sourced from domesticated apples, is more efficient than the loquat variant. Our stochastic simulations showed that manipulating individual enzymes is not the most efficient strategy as the specificity required for these substrate-level approaches would not be feasible to develop during the iGEM timeline. However, our network topology analysis of the greater metabolic pathway highlighted a few practical wet lab approaches. After simulating gene knockouts, analyzing degree distributions and clustering coefficients, quantifying system robustness, and augmenting with visual analyses, we identified the best methodology to optimize UA production. In this method, we knock out the voltage-dependent channel POR2, which transports acetyl-CoA from the cytoplasm into the mitochondrial matrix. We summarize our key activities below:

Introduction

Team Cornell 2024 is looking to develop Oncurex, a scalable platform to increase the supply of ursolic acid treatment. Ursolic acid (UA) is a triterpenoid commonly found in apples and loquats that has garnered attention due to its anti-cancer activity.

Though therapeutically promising, UA production is unsustainable as it relies on extraction from fruit, and the chemical synthesis of ursolic acid is low-yield. Current research shows that Saccharomyces cerevisiae can produce 2,3-oxidosqualene, a precursor to UA. S. cerevisiae has all the machinery needed to produce UA except for one enzyme: ɑ-amyrin synthase. S. cerevisiae can be genetically engineered to express ɑ-amyrin synthase. However, simply expressing a single enzyme will not drastically enhance ursolic acid production. For Oncurex to be a competitive scale-up company, we must develop a method to further optimize production. Our wet lab subteam developed two models as we explored the available methods.

All wet lab model development is based on the following metabolic pathway. This pathway combines research from Jia et al. [1], Jordá and Puig [2], and Jin et al. [3].

Figure 1. Metabolic map of Saccharomyces cerevisiae encompassing glycolysis, gluconeogenesis, pentose phosphate pathway, citric acid cycle, mevalonate pathway, β-oxidation, amino acid synthesis/degradation, nucleotide synthesis, and ursolic acid synthesis. Compartmentalization of pathways represented by blue (cytoplasm) and green (mitochondria)

Substrate-Level Michaelis-Menten Kinetics

Theoretical Development

Michaelis-Menten kinetics can be used to gain a basic understanding of the activity of enzymes in our systems. Dalwadi et al. utilized Michaelis-Menten kinetics to develop and solve a model for the mevalonate pathway [4]. Following this, we assume Michaelis-Menten kinetics for all reactions considered in wet lab modelling. This is a mathematical representation of the lock and key model of enzyme mechanics.

Figure 2. Theoretical development of Michaelis-Menten kinetics. The differential equation dP/dt represents the change in product formation based on enzyme and substrate concentration, and the maximum reaction rate. This model is consistent with the lock and key enzyme mechanics analogy.

Michaelis-Menten kinetics relies on three assumptions:

1. Steady-state approximation: The rate of enzyme-substrate complex formation are equivalent to the rates of enzyme-substrate breakdown. In reality, metabolism is dependent on environmental/physical factors
2. Free ligand assumption: The concentration of substrate is higher than that of enzyme. In reality, this is quite accurate to metabolism
3. Rapid equilibrium assumption: Reactions only move forward. In reality, enzymes may catalyze reverse reactions, which can be mathematically accounted for by multiplying by -1

We decided to model kinetics as a series of differential equations in MATLAB due to its extensive library of numerical solvers. We used ode45 to solve all equations. MATLAB also has high quality figures and an easy to use IDE (interactive design environment).

Kinetics of ɑ-Amyrin Synthase

We first utilize Michaelis-Menten Kinetics to characterize important enzymes located in the pathway right before the synthesis of ursolic acid occurs, dubbed the metabolic neighborhood:

Figure 3. Metabolic neighborhood of ursolic acid production with defined enzyme names.Lanosterol formation competes with ɑ-amyrin as both use the substrate squalene epoxide (2,3-oxidosqualene). The final product in the neighborhood, synthesized by ERG24, is compared to the production of UA.

Our focus is on increasing the efficiency of UA production, and given that we are exogenously expressing ɑ-amyrin synthase (AAS), we took the opportunity to identify the most efficient variant. We considered both variants of OSC1 synthases from Malus x domestica and Eriobotrya japonica, which are domesticated apples and loquats, respectively. MdOSC1 is the gene that codes for AAS.

The code for this model considers all reactions shown in Figure 3, though there are some caveats between Figure 1 and Figure 3.Our model for the OSC2 branch in FIgure 1 is informed by Yu et al.., who found that MdOSC1 produced α-amyrin, β-amyrin, and lupeol at a 86:13:1 ratio [5]. Thus, we modelled α-amyrin production as a function of the production of all three compounds mentioned.

We also consider the conversion of squalene to 2,3-oxidosqualene (squalene epoxide) implicitly. With each time step, there is a small amount added to the initial concentration, simulating replenishment of the substrate pool.

Figure 4. Considerations made in the first kinetic model of ursolic acid production. All modelling choices made in consultation with literature or with stakeholder feedback (Ex: Dr. Ben Cosgrove). In reality, the kinetic balance is influenced by external factors and non-diminishing substrate concentrations.

This model further considers the reaction catalyzed by CYP450 and CPR as Michaelis-Menten. Though studies differ on the kinetics of this cytochrome, a study performed by Shumyantseva et al. showed that the steady-state assumption holds for this complex reaction [7]. However, CYP450 has not been characterized with α-amyrin, so to source kinetic parameters, we looked at a reaction between CYP450 and the steroid hormone pregnenolone. We make this argument based on the thermodynamics of enzyme-substrate complex formation, where the functional groups need to maintain a specific orientation to minimize binding energy.

Lastly, with parameters established for the direct pathway of UA synthesis, we also include the effects of competing pathways. Lanosterol synthase can catalyze 2,3-oxidosqualene into lanosterol, directly competing with the activity of α-amyrin synthase. Similar to the ursolic acid pathway, an enzyme in the P450 family catalyzes a reaction which the Lanosterol pathway, where the lanosterol 14-alpha-demethylase catalyzes a reaction with lanosterol to form 4,4-Dimethylcholesta-8, 14, 24-trienol (which will be referred to as trienol). With this view, the lanosterol synthase and lanosterol 14-alpha-demethylase are analogous to the alpha amyrin synthase and amyrin c-28 oxidase enzymes in the UA pathway. Unlike the exogenous pathway which ends after two reactions, lanosterol forms ergosterol, the yeast equivalent of cholesterol.

Figure 5. Comparison of OSC1 from domestic apples (Malus x domestica) and loquats (Eriobotrya japonica). kcat value for MdOSC1 is ~3x greater than EjAS. Similarly, the Km value is double that of EjAS, indicating that on a kinetics argument, MdOSC1 is more efficient.

Based on Figure 5, we conclude that MdOSC1 and EjAS have opposing behaviors. EjAS’s high Km value indicates that the enzyme has a much lower affinity for 2,3-oxidosqualene. Thus, the initial reaction rates will be slower as the system reaches steady-state over 4 minutes. However, the MdOSC1 variant reaches steady-state in just 2 minutes, a 50% reduction in time which is significant on an industrial scale. Furthermore, the low Km and high kcat values indicate that this enzyme has been evolutionarily optimized and requires no further modification (through methods such as rational protein design). Therefore, the MdOSC1 variant is more efficient and more optimal for wet lab experimentation.

Stochasticity of MdOSC1 Reactions

Like many real-world systems, biological systems are highly complex and difficult to analyze using pure deterministic methods. Biology is inherently random and complex, so we must pursue alternative methods to understand how systems are affected by this noise. To expand on the deterministic solution established in the previous section, we modelled the metabolic neighborhood as stochastic.

Figure 6. Stochastic models account for biologically randomness by using discretized, or step-by-step, iteration. In each iteration, the reaction that occurs is dictated by the propensity of said reaction. Propensity is dependent on kcat and the number of molecules involved in a reaction. By using propensity functions correctly, you capture the randomness of a system evolving over time.

To account for stochasticity, we converted our Michaelis-Menten equations into propensity functions. Then, we adapted the Gillespie algorithm, also known as the stochastic simulation algorithm. This model is a powerful tool for enzyme kinetics because it provides a detailed time course of a system’s evolution. It can also accommodate complex reaction networks including multi-substrate reactions, allosteric inhibition, and reversible reactions (the key is in how you define the propensity functions). Despite these algorithms being computationally intensive, we followed the below workflow to create a stochastic solver.

Figure 7. Pseudocode used to implement stochastic solver algorithm. Reaction that occurs is chosen via random number generator. Time (tau) represents the change in time between the previous step and the next. It is chosen by a second random number generator plugged into a log(time) function. The stochastic solver ends once a solution converges or if a maximum amount of steps has been reached

This workflow provides the steps required to create one complete iteration of the Gillespie algorithm with n steps and n-1 tau values. We can take this one step forward by adapting the Gillespie algorithm into a Monte Carlo simulation. Monte Carlo Methods are a broad class of algorithms that rely on repeated random sampling when generating numerical results. In other words, a Monte Carlo simulation involves repeatedly using the Gillespie algorithm to generate many simulations, each with many samples. This probabilistic view rightfully models enzyme kinetics as a random process filled with uncertainties and biological noise We used the kinetic parameters for MdOSC1, CYP450 (amyrin C-28 oxidase), CYP450 (lanosterol 14-a-demethylase), and lanosterol synthase from Figure 5. We assume that all metabolites have a starting concentration of zero, except for the starting metabolite of 2,3 Oxidosqualene, while also prescribing enzyme concentrations to be constant. Since 2,3 Oxidosqualene levels will be replenished, our model takes into consideration that some of this substrate will enter the system at each time step. Given these various assumptions, the main takeaways of this stochastic model will be the relative final concentrations of metabolites. We also plot many outcomes of this simulation to demonstrate the wide range of possible final concentrations.

Figure 8. Results of stochastic modelling of ursolic acid production against the competing lanosterol pathway. Yellow = a-amyrin, green= lanosterol, black = average of each data set (bottom = a-amyrin, top = lanosterol). Lanosterol production in a majority of cases is higher than that of a-amyrin.

The simulation shows the propensity for lanosterol to be formed over α-Amyrin and any resulting ursolic acid in a 5:3 ratio. As such, if knocked out, we predict that ursolic acid production could increase two or three fold with an increase in available 2,3 oxidosqualene. However, on closer examination of the lanosterol pathway (shown in Figure 1), we see that it is the precursor to ergosterol which functions similarly to cholesterol in animal cells. Ergosterol is essential to yeast survival, and there have been studies that show that knocking out lanosterol synthase can severely suppress cell growth [7] or be fatal [8]. Thus, we rule out lanosterol synthase knockout. However, we conclude that knockdown of lanosterol synthase could be effective depending on how S. cerevisiae responds to the lower levels of the triterpenoid.

Sensitivity Analysis of Metabolic Neighborhood

Expanding on the Monte-Carlo method from the previous section, we expand it into a sensitivity analysis, also known as a “what-if analysis.” In this analysis, we manipulate the initial conditions and rate constants to test which parameters have the most significant impact on the output. Our method of choice to conduct our sensitivity analysis is parameter sweeps. This is where we run simulations with different ranges of parameters to explore how the system behaves under various conditions. Note that a parameter sweep can be likened to a Monte Carlo simulation. That is, in each iteration of a Monte Carlo simulation, we change the enzymatic parameters, allowing us to generate large data sets.

Figure 9. The goal of this sensitivity analysis is to provide further proof that MdOSC1 is an optimal candidate for high efficiency production of ursolic acid. To do this, we use parameter sweeps of kcat and Km values to identify how the ratio of ursolic acid/trienol changes. These parameter sweeps simulate methods that directly edit enzyme function such as site-directed mutagenesis.

Though our stochastic simulation describes the activity of 4 enzymes, they can be analyzed in pairs (Figure 10a/b). Since the alpha amyrin synthase and lanosterol synthase both compete directly for 2,3 oxidosqualene as substrate, their activities can be quantified by the final ratio of alpha amyrin and lanosterol, their respective products. This provides a direct comparison between the two competing pathways. Similarly, the pair of CYP450 enzymes in Lanosterol 14-alpha-demethylase and Amyrin C-28 Oxidase can be analyzed together since they catalyze reactions with the products of alpha amyrin synthase and lanosterol synthase (Figure 11a/b). To quantify their activity, the ratios of the products ursolic acid and trienol are tracked.

Figure 10a. Effect of Changing Km and Kcat Values for a-amyrin synthase on the UA Producing Pathway. The ratio of ursolic acid/trienol is represented by the ratio a-amyrin/lanosterol, the names corresponding to the pathways that catalyze the formation of the product. Ursolic acid concentration is high when alpha amyrin synthase has a high \(K_\text{cat}\) and low \(K_\text{m value}\)

Figure 10b. Effect of Changing Km and Kcat Values for lanosterol synthase on the UA Producing Pathway. Ursolic acid concentration is high when lanosterol synthase has a low \(K_\text{cat}\) and low \(K_\text{m value}\)

Figure 11a. Effect of Changing Km and Kcat Values for alpha C-28 oxidase (CYP450) on the UA Producing Pathway. Ursolic acid concentration is high when this enzyme has a high \(K_\text{cat}\) and low-medium \(K_\text{m value}\)

Figure 11b. Effect of Changing Km and Kcat Values for lanosterol 14-a-demethylase (CYP450) on the UA Producing Pathway. Ursolic acid concentration is high when this enzyme low \(K_\text{cat}\) and \(K_\text{m value}\)

Similar to the previous pair of enzymes, the activity of alpha c-28 oxidase and lanosterol 14-alpha-demethylase are dictated by their kcat value. When the kcat value of alpha c-28 oxidase is increased, the ratio of ursolic increases and when the kcat value of Alpha C-28 Oxidase , with the opposite relationship for holding for lanosterol 14-alpha-demethylase. Just like the lanosterol synthase, the effect of lanosterol 14-alpha-demethylase depends solely on its kcat value. Upon further scrutinization of the ratios, the ursolic acid-trienol ratio is incredibly high for the naturally occurring Alpha C-28 Oxidase enzyme in yeast at a ratio of 3.5. This demonstrates the comparatively low rate of reaction for the lanosterol pathway, proving that lanosterol 14-alpha-demethylase is one of the rate limiting reactions in the ergosterol pathway [2].

Conclusion - Substrate-Level Michaelis-Menten Kinetics

Michaelis-Menten modelling of potential alpha amyrin synthase strains were used to determine which strain should be expressed in yeast, settling on the MdOSC1 and its related strains as the best candidate given their high expression and high enzymatic activity. With a strain determined, it was possible to run stochastic simulations to determine the activity of the aforementioned 4 enzymes in the same system. The competitive inhibition between the alpha amyrin synthase and lanosterol synthase was clear, with the alpha amyrin to lanosterol ratio of 0.6. Upon further scrutinization of the ergosterol pathway, it is clear that wet lab will have to strike a delicate balance between ursolic acid production and lanosterol production, since lanosterol is the precursor for ergosterol, a metabolite necessary for the survival of yeast.

It was then possible to run sensitivity analyses where stochastic simulations were repeatedly run with changes in parameters each time, giving us a better understanding of the activity of these enzymes. In the future, It is possible for wetlab to design a repression or promoter system depending on if yeast cells are not viable (too much ursolic is being produced) or healthy (could increase expression of ursolic acid producing enzymes). Therefore, by combining deterministic and stochastic modelling approaches, we were able to gain a deeper understanding of the complex interactions between enzymes and substrates in this metabolic pathway.

Besides the limitations of using the Michaelis-Menten model to describe these enzymes, this model of the synthetic pathway only takes into account a small portion of the entire metabolic pathway of yeast. There are likely key inhibitors and NADPH/NADP+ balance necessary to fully understand the activity of the synthetic ursolic acid pathway within the ergosterol pathway.

While our analysis suggests that manipulating enzymes involved in the MVA pathway, specifically alpha C-28 oxidase, could potentially increase UA production, enzymatic manipulation is not the most practical approach. For instance, altering the Km or kcat values of enzymes through techniques like site-directed mutagenesis can be challenging and time-consuming. Additionally, such modifications could have unintended consequences on other metabolic pathways or cellular functions, leading to undesirable side effects. Furthermore, the complex regulatory mechanisms that govern the mevalonate pathway may make it difficult to achieve sustained and predictable changes in enzyme activity through direct manipulation. Therefore, while enzymatic manipulation offers a theoretical possibility, it may not be a viable strategy for increasing ursolic acid production in practice.

Figure 12. Enzyme kinetics is one part of global metabolism. Analyzing kinetics under the Michaelis -Menten assumption yielded useful information, but unfortunately, it is not very actionable. Thus, we consider global metabolism, which accounts for more than just kinetics

To this end, we decided to step back from enzyme kinetics and consider a global metabolic perspective using network topology analysis. We believe this will help identify points of intervention that can be exploited to optimize UA production.

S. cerevisiae Global Metabolism Analysis

The above Michaelis-Menten modelling, and attached Monte Carlo simulation, provide a view of reaction kinetics. However, like we said in the conclusion of the previous section, it is difficult to manipulate enzymatic properties with such a high level of specificity that the deterministic and stochastic models suggested.

Therefore, we consider the larger metabolome of Saccharomyces cerevisiae. Since analyzing the reaction kinetics of multiple converging pathways would be computationally intensive, we instead decided to interrogate metabolic topology. Network topology analysis is a method of modelling the structure of a network to describe how information flows. As metabolism is a complex series of interconnected reactions, we don’t need to understand the kinetics of every reaction. However, we can benefit from identifying bottlenecks/key components (called hubs), understanding system robustness and simulating knockouts by adjusting the modelled topology.

To this end, we decided to pursue an integrated approach where we performed a network topology analysis of relevant pathways that connect to our metabolic neighborhood. Figure 1 contains a collection of metabolic pathways that can affect ursolic acid production, either directly or indirectly. We decided to focus on the following pathways:

Figure 13. List of all the metabolic pathways represented in Figure 1. Note that in Figure 1, no inhibition or cofactors/co-enzymes were shown due to a limitation in our model development. The crossing-over between these pathways affects the overall production of ursolic acid.

Theoretical Development

Before we get into network topology analysis, we must first understand graph theory, a branch of mathematics that studies graphs, which are networks that model pairwise relationships between objects. The following schematic represents a basic graph where there are many possible connections between the different circles:

Figure 14. A schematic of a graph where each line represents a pairwise relationship between two circles. Note that each circle is a separate entity and the only connections between circles occurs through a line connection.

As mentioned before, network topology analysis looks at how information flows between the elements of a system. Graph theory operates under the following assumptions and scope:

1. Abstraction: Representing complex pathways as a system of nodes and edges. In reality, details such as cooperativity, allosteric inhibition, and other complex behavior occurs, information that is lost in this representation
2. Discretized elements: Each metabolite is represented as a discrete node with discrete edges. In reality, metabolite concentrations are continuously dictated by reaction rates
3. Pairwise relationships: One reaction (edge) converts one substrate to one product. In reality, many reactions involve multiple substrates and products. Though edges can represent substrate-product relationships, complex mechanisms are not represented such as cooperativity/ultrasensitvity
4. Static Structure: Graph does not change over time. In reality, metabolic networks are affected by cell signaling, changes in gene expression, etc
5. Finite: Graphs represent a finite size of nodes and edges. In reality, S. cerevisiae's metabolism is much larger than is represented in this graph

Figure 15. Visualization of the loss of information between chemical reactions and graph theory representation. This reaction involving pyruvate dehydrogenase is a complex, yet essential reaction, to aerobic respiration. In graph theory representations, cofactors/coenzymes cannot be involved. Nor can complex mechanisms and sub-reactions that occur within the PDH complex. Though some information is lost, the relationship between pyruvate and acetyl-CoA provides information on metabolite flow, which is useful when looking for genetic knockout candidates.

In a metabolic engineering application, each graph consists of a set of nodes (points which represent metabolites) and edges (connections between nodes that represent the flow of information, or in this case, a reaction). An important consideration is that all reactions do not occur with the same likelihood. Genetic and environmental factors heavily influence what reactions occur and when. However, we cannot model these effectively in a graph because they can change at non-discrete times. Since a graph is static by nature (assumption #4, illustrated in Figure 15), we need a static metric. Therefore, we weighted edges based on their Gibbs Free Energy value (ΔG). As a measure of spontaneity, a more negative ΔG value indicates a stronger likelihood that a specific reaction will occur.

As an aside, we were presented with an alternative modelling approach to the greater metabolome. During an interview with Dr. Ben Cosgrove of Cornell Biomedical Engineering, he suggested we look into flux balance analysis (FBA), which is a computational method that predicts metabolic fluxes under steady-state conditions. It is based on the assumption that the system is at metabolic steady-state, meaning that the rates of production and consumption of each metabolite are balanced. However, with the introduction of ΔG weights, our network topology analysis uses basic principles of flux balance analysis.

Figure 16. Our approach to modelling the yeast metabolome incorporates elements of flux balance analysis within standard graph theory development. Our model is capable of predicting metabolic flux since we constrain our model using ΔG weights. This allows for us to make targeted calculations.

To summarize, we incorporated elements of flux balance analysis into our graph theory approach. We constrained our model using ΔG weights which provide an objective measure of the likelihood a reaction occurs without needing to represent the complexity involved with genetic and environmental influences.

Figure 17. Yeast-GEM is a repository of relevant data regarding Saccharomyces cerevisiae. The data it contains can be used for many applications, listed above

Graph Initialization

Under a specific temperature/pH condition, the ΔG value is constant which further satisfies the requirement that the weight must be static. Knowing that S. cerevisiae’s optimal temperature is 37°C, we calculated the ΔG at this temperature and assumed pH 6.

We sourced the values from yeast-GEM, or yeast genome-scale metabolic model [9]. It features detailed mappings of metabolic pathways and interactions and can be used for a wide variety of systems biology projects. The information of interest to us were the simulated thermodynamic coefficients. These values were determined by linking reactions based on the identified metabolites and enzymes. Then, the model is calculated and validated against experimental data.

The values we sourced are in the following figure:

Figure 18. Table documenting the enzyme names, substrates, products, and delta G values used in our model. The OSC1, OSC2, and CYP450/CPR values were not found as OSC1/2 are not endogenous to yeast.

Utilizing the information from the above table, we built a custom representation of the portion of the metabolome we are considering using the Jung Graph Builder in Java. This library provides a flexible framework for building complex graphs, performing graph algorithms, and rendering the results.

Figure 19. Jung Graph representation of the metabolic map. Blue circles are nodes that represent reactants. Grey lines are edges, which represent reactions. The thickness of the edge correlates with the ΔG value. The more spontaneous a reaction, the thicker the edge appears visually.

Referencing Figure 19, we will define some key terms:

• Undirected graph: flow of information between nodes is bidirectional
• Directed graph: information flows in one direction only between nodes
• Weighted graph: a value that represents the strength/cost of an edge
• Degree: the number of edges a node is connected to
• Path: a sequence of edges that connects 2 directly OR indirectly nodes

Below we define the strengths/uses of network topology analysis:

1. Centrality: Identifying metabolites that play a crucial role based on degree. Identifies hubs and bottlenecks which indicate that these compounds have important roles in the overall pathway we are modelling
2. Shortest Path Analysis: Finding the most efficient routes for metabolic conversions. This approach is the basis for all graph-based analyses
3. Simulating gene knockouts: Removing nodes and edges to understand how paths change when looking at ways to optimize UA production. We can do this by manually removing them in the adjacency matrix
4. Robustness: Analyzing the ability of a metabolic network to withstand external changes. This can be done by simulating KOs, adding new metabolites or increasing the temperature of the system (which changes the ΔG values), etc.

Simulating KOs and Quantifying Shortest Path

Specific information, like the degree distribution of a network, can also be extracted from graphs. Degree distribution plots, like the one in Figure 20, provide information about the connectivity of the network. Plots with a binomial distribution indicate that the nodes have similar degrees, while distributions with a long tail (power law distribution) indicate the presence of hubs, which are nodes with a much higher degree than the other nodes in the network. This information is useful for optimizing a metabolic pathway, as changes made to a hub likely would have impactful changes on the overall network.

Figure 20. Degree distribution plot of graph represented in Figure 19. Many nodes have a degree of 2 which indicates that the graph has a high amount of linearity present (one formation reaction, one consumption reaction). However, the presence of hubs indicates that the degree distribution follows a power law distribution.

The nodes with the greatest degrees are acetyl-CoA and pyruvate, both with six neighbors. There are two ways to interpret the presence of these hubs. First, hubs can make the network more resilient to perturbations as the provide alternative pathways for metabolites to flow. However, hubs also represent high importance. This means that the hub is central to the function of the metabolic pathway.

Figure 21. Zoomed in view of Figure 1, focusing on pyruvate as a hub. Pyruvate’s metabolic neighborhood reflects the compound’s centrality in the metabolic pathway modelled.

Thus, there are two approaches here: leave the hub alone as it is central to the pathway or manipulate reactions near the hub and quantify changes in metabolic flux. Considering that the two 6-degree hubs are pyruvate and acetyl-CoA, both of which are central to many metabolic pathways, targeting reactions near it is a more viable approach. Directly manipulating pyruvate kinase/pyruvate dehydrogenase or modulating the concentration of acetyl-CoA directly is a high-risk strategy and does not guarantee success. Therefore, we explored KOs of the edges (reactions) connected to these two hubs.

To simulate a knockout, the edges connected to acetyl-CoA and pyruvate were deleted from the graph representation from Figure 19. We compared the distance of the shortest path, as well as the overall ΔG value, to quantify whether there is a difference in efficiency between a normal version and a KO.

Figure 22. Modifying Figure 14 into a weighted graph. Node 1 is the start and Node 6 is the destination. Thicker edges indicate these edges are easier to traverse (lower ΔG value in the biochemistry analog).

Knowing that our graph is weighted, we used Dijkstra's algorithm to calculate the shortest path. This algorithm identifies the shortest path via the 4 rules stated in Figure 21. Starting at the 1st node (or what was identified as the first node), the algorithm will systematically traverse edges adjacent to the current node. Each edge is associated with a score (ΔG value in the biochemistry analog). However, following rule #1, only positive weight values could be used. Our negative weight with the greatest magnitude had a Gibbs free energy value of -448.19 kJ/mol, so we added 449 kJ/mol to each weight.

To calculate the true ΔG, we back-calculated by subtracting (449 kJ/mol) * (# of nodes in path) from the adjusted distance. From there, we can make comparisons regarding the effect of each KO.

In the ΔG score calculation, we made the following assumptions for the ΔG value of enzymes/transporters we could not find ΔG values for in yeast-GEM:

• Oxidosqualene cyclase = 10 kJ/mol
• CYP + CPR = 10 kJ/mol
• Phosphoserine aminotransferase = 10 kJ/mol
• Acyl-CoA Synthetase = 10 kJ/mol
• Citrate transport = 10 kJ/mol
• Pxa1p, Pxa2p = 10 kJ/mol
• POR2 = 2.5 kJ/mol

The values for all enzymes/transporters except for POR2 were made based on the observation that the average edge weight was around 10 kJ/mol. Without knowing the true extent of the effect of these edges, we assume that their effect is similar to the rest of the metabolome modelled. POR2, which codes for a transporter that moves acetyl-CoA into the mitochondrial matrix, has a lower ΔG because S. cerevisiae will be cultured in a controlled bioreactor. Thus, it should always undergo aerobic respiration because there should always be a supply of O2. Based on this, the POR2 transporter will be more active as it shuttles acetyl-CoA into the mitochondria for the citric acid cycle. To reflect this preference, POR2's ΔG is 2.5 kJ/mol.

Figure 23a. Pyruvate Node of metabolic graph representation. 3 nodes produce pyruvate (pyruvate kinase and malic enzyme, serine dehydratase) while 3 nodes consume pyruvate (pyruvate dehydrogenase, pyruvate decarboxylase, pyruvate carboxylase).

Figure 23b. Acetyl-CoA node of metabolic graph representation. 3 edges produce acetyl-CoA (citrate synthase, pyruvate dehydrogenase, acetyl-CoA synthetase) and 3 edges consume it (ERG10, POR2 transporter, acetyl-CoA C-acetyltransferase)

Looking at Figures 23a and 23b, we considered knockouts of all edges that connect to pyruvate and acetyl-CoA; the results of the pyruvate KOs are listed below:

Figure 24. Table comparing changes in path length and distance (kJ/mol) based on knockouts of neighboring enzymes. Pyruvate dehydrogenase (PDH) increases path length in almost all cases. Stark increases in adjusted distance for 4 out of 5 starting nodes in PDH KOs. Knockout of malic enzyme, pyruvate carboxylase, and pyruvate decarboxylase have no effect on shortest path

Similar to pyruvate, all edges connected to acetyl-CoA were knocked out (edge removed in graph), and then the path lengths and distances were calculated. Pyruvate dehydrogenase was no longer under consideration as a knockout given that its knockout increases the path distances from glucose, ribulose-5-phosphate, serine, and glycine. Knockout of ERG10 blocks all five starting nodes from 2,3-Oxidosqualene. Knocking out citrate synthase blocks the path from fatty acid to 2,3-Oxidosqualene. Neither knockout of Acetyl-CoA synthetase, acetyl-CoA carboxylase, nor POR2 transporter changes any of the shortest path distances.

Thus, the following enzymes were left under consideration because they do not increase the path distance from any of the five starting nodes to 2,3-Oxidosqualene:

• Malic enzyme
• Pyruvate carboxylase
• Pyruvate decarboxylase
• Acetyl-CoA synthetase
• Acetyl-CoA carboxylase
• POR2 transporter

However, this list can be further specified. We know that these enzymes do not affect path length and have minimal effect on distance/adjusted distance. To gain a more comprehensive understanding of the metabolic network and identify potential targets for intervention, we combined our quantitative measurements with a qualitative assessment of the graph.

Recall the two arguments made about the degree distribution: hubs can make the network more resilient to perturbations or their centrality is too important to change. The nodes pyruvate and acetyl-CoA synthase are integral to the network, and malic enzyme and pyruvate carboxylase catalyze the formation of pyruvate, so it makes little sense to knock them out.

Figure 25. Zoomed in view of Figure 1, focusing on acetyl-CoA as a hub. Acetyl-CoA’s metabolic neighborhood reflects the compound’s centrality in the metabolic pathway modelled.

While pyruvate decarboxylase does reduce the amount of pyruvate, the pyruvate decarboxylase edge leads to acetyl-CoA after two nodes, so its knockout likely would not have a large impact. Similarly, acetyl-CoA synthetase catalyzes the formation of Acetyl-CoA, so its knockout would be unhelpful. Acetyl-CoA carboxylase and POR2 (Porin 2) transporter both decrease the amount of acetyl-CoA. However, POR2 transporter is notable because it is a bottleneck. That is, the transport of acetyl-CoA into the mitochondria is a necessary step that leads to the citric acid cycle and oxidative phosphorylation.

Figure 26. Effect of POR2 knockout on acetyl-CoA metabolic flux. POR1 can partially compensate for POR2 knockout while the remaining acetyl-CoA is redirected to cytosolic pathways, such as the MVA pathway

POR1 and POR2 are isoforms of voltage-dependent anion-selective channel (VDAC), also known as eukaryotic mitochondrial porin [10]. POR1 is the more abundant isoform, indicating a more important role, presumably one that extends beyond the scope of Oncurex. However, given that POR2 has similar electrophysiological features to POR1, the only difference is that POR2 expression is lower [10]. Knocking out POR2 will redirect metabolic flux from the mitochondria into the cytosol, where the acetyl-CoA will then be available for use in the cytosolic mevalonate pathway. Therefore, exploiting this bottleneck and knocking out POR2 is a good strategy to optimize UA production.

Conclusion - S. Cerevisiae Global Metabolism Analysis

This model presents a comprehensive analysis of the metabolic pathways involved in ursolic acid production in Saccharomyces cerevisiae.

Through the integration of graph theory and flux balance analysis, we have developed a robust model that accurately represents the complex interactions within the yeast metabolome. The use of ΔG weights as a constraint has allowed us to evaluate the likelihood of reactions occurring without the need for detailed kinetic information.

Our analysis revealed that the POR2 transporter is a promising target for optimizing ursolic acid production. By knocking out POR2, we can redirect metabolic flux towards the cytosolic mevalonate pathway, leading to increased ursolic acid synthesis. Our model confirms this as we identified POR2 as a bottleneck in yeast metabolism. We also note that knocking out the POR2 transporter does not increase the path length or adjusted score yielded by Dijkstra’s algorithm. Thus, though the reaction may not be faster, it will occur more due the redirection of acetyl-CoA from the mitochondria to the cytosol.

While this study provides valuable insights, further research is needed to validate our findings experimentally. Additionally, exploring other potential targets and investigating the effects of various environmental conditions on ursolic acid production could further enhance our understanding of this metabolic pathway.

Product Development

Executive Summary

Oncurex is an industrially scalable platform to increase the supply of ursolic acid for breast cancer treatment. The product development subteam planned to develop a continuous stirred-tank reactor (CSTR) to culture the genetically modified yeast made by the wet lab subteam. To direct our approach, Team Cornell used computational and mathematical models to understand the bioreactor system:

1. Modeling the Continuous Stirred-Tank Reactor System: Mathematically characterizing the bioreactor in a lab-scale setting
2. Modeling Ursolic Acid Detection in the Biosensor: Designing and analyzing how the biosensor detecting Ursolic Acid
3. Industrial Scale-Up: Designing an industrially-scaled version of the bioreactor system that characteristically matches the lab-scale version

Introduction

Our Product Development team designed a continuous stirred-tank bioreactor (CSTR) containing genetically engineered yeast to efficiently produce Ursolic Acid. In addition, an integrated biosensor was designed to track Ursolic Acid yield in real time. To understand the system of a CSTR while building it, modeling was utilized to guide our approach to scale-up. By having a comprehensive understanding of our CSTR, we can assess the viability of the model in lab-scale and industry-scale. We first state the simplifying assumptions we made for the models, such as well-mixed reactor contents and steady-state conditions. We then implement Michaelis-Menten kinetics to model the rate of ursolic acid production and rate of lanosterol production. We use the Monod equation to model the growth rate of yeast, and combine this information with yeast death rate and inflow and outflow of yeast from the CSTR to generate a model of the total change in yeast concentration. Our steady-state assumption means that we can set the rate of change of the concentration of yeast in the CSTR to 0. We next analyze how various assumptions have affected our model. We used a software called Visimix to model our CSTR and generate a Reynolds number to determine whether our reactor would support laminar or turbulent flow. During testing, we obtained a higher-than-expected Reynolds number indicating turbulent flow, so we analyze possible causes of this discrepancy. We calculated a Biot number to view heat transfer conditions in our bioreactor and conclude that convection heat transfer dominates conduction heat transfer. This Biot number will be essential in future scale-up.

Modeling the Continuous Stirred-Tank Reactor System

Visimix Lab-Scale Computational Modeling for CSTR Mixing

Visimix is a specialized software for modeling mixing processes, to develop a computational model of the CSTR. When constructing the model, several parameters were defined and calculated to ensure that the reactor functioned optimally for our application.

We measured and set the geometry and volume of the reactor as a cylindrical tank with a diameter of 15 cm and media volume of 1 liter. Additionally, to simplify the model, the density and viscous properties of the media were set to match water.2 The impeller position from the bottom of the tank was also measured and set to ensure effective mixing.

Figure 27: 2D model of the CSTR in Visimix with a diameter set to 15 cm and media volume of 1 liter.

The impeller size and speed were critical parameters to set for the CSTR. We wanted to keep the speed at an optimal velocity where fluid flow would stay laminar while considering homogeneity within the media and promoting cell growth. By keeping the flow laminar, we would be able to precisely control the other parameters in the system, increase the efficiency of the chemical reactions, and reduce the risk of shear damage on yeast cells. As mixing has a large effect on fluid flow, a dimensionless number called the Stirred Tank Reynold’s number was calculated where Ni is the stirrer speed, Di is the impeller diameter, ρ is density of the fluid, and μ is viscosity.

This number models the fluid flow of our system and determines if our system is laminar or turbulent. We decided that a Reynold’s number less than 2000 would be ideally implemented so that steady state within the reactor would be maintained while also having enough circulation in the reactor to promote cell growth. Furthermore, Reynold's numbers under 4,000 would indicate laminar flow.

Multiple speeds of the impeller were tested and we decided that 50 RPM with a propeller diameter at 7.5 cm would have a RCF equal to 0.21 g. This would be the most optimal in ensuring low speed operation that would support laminar flow and bring ordered flow patterns. The power consumption of the impeller was calculated and measured to be 2.5 Watts which is optimal for low speed stirring and laminar flow. As this is a continuous reactor, we had to determine the appropriate feed and discharge flow rates to maintain steady-state operation. The flow rate of media into the reactor is also an influencer in laminar or turbulent flow so we tested and determined that a feed and discharge rate of 20 mL/min would ensure a slow and steady supply of nutrients and cells into the system.

Initially, we aimed to achieve laminar flow in the reactor to maintain precise control and prevent shear damage to the yeast cells. However, during testing, we unexpectedly calculated a Reynolds number of ~33,000, indicating turbulent flow instead of the intended laminar regime. The numbers utilized to calculate this number are listed below:

This discrepancy from our Visimix model, which predicted a Reynolds number below 2000, can be attributed to several factors. First, there may have been inaccuracies in the input parameters, such as the actual viscosity of the media, which may have been higher than the values used in the model. To keep Reynold's number at a reasonable number, the impeller speed during operation was raised from 0.21 RCF to 0.42 RCF. Additionally, the reactor's geometry, particularly the impeller size and placement, might have caused uneven flow patterns that led to higher turbulence. Another major assumption that may have contributed to this discrepancy would be the presence of the multiple sensors and oxygen diffuser that disrupt laminar flow. Lastly, external factors, such as temperature fluctuations or media composition, could have affected the fluid properties, increasing the flow velocity and contributing to the higher Reynolds number.

Through Visimix, we calculated a Reynold's number of 1930, ensuring a stable, low-shear environment for growing genetically modified yeast. However, the test run discussed previously suggests that further adjustments are needed to our model to reflect what is occurring in real life. While laminar flow is ideal, turbulent flow is not necessarily bad either as many industrial bioprocesses can be turbulent as well . This taught us that in future renditions of the Visimix model, turbulent conditions can be assessed and optimized for yeast growth to ensure optimal Ursolic Acid production and tracking.

Heat Transfer within the CSTR using Biot Number

While building the bioreactor, we recognized that heat transfer is an integral characteristic to understanding what is occurring within our system. For optimal yeast growth and Ursolic Acid production, it was determined the bioreactor had to be kept at 37°C. However, due to the nature of the system, we immediately became curious about the heat transfer between the inside and outside of the main reservoir. To model this, the Biot number was used with h representing the convective heat transfer coefficient in \(W/(m^2K)\), k representing thermal conductivity of the main reservoir in \(W/(mK)\), and L representing the height in meters of the reactor.

The Biot number is a dimensionless number that models the ratio of conduction and convection within a model, allowing us to view heat transfer conditions within the bioreactor. Assuming conduction and convection heat transfer are significantly larger than radiation, we can use this number to model the heat transfer of the system. It was also assumed that the reactor is filled and therefore heat is being transferred at a consistent rate. Using these assumptions, a Biot number was calculated to be 22.78 using the values listed below.

With a number above 1, this tells us that convection heat transfer dominates conduction heat transfer within our bioreactor. As we want to maintain a consistent temperature within the media in the reactor, this is a good sign to see as heat will be mainly transferred within the media rather than being lost through the walls. However, any loss of heat can be fixed using a sous vide machine within our system to provide heat transfer from a water bath external of the reactor to the inside. With this Biot number, we can also begin the industrial-scaling of our bioreactor. The Scale-Up and Future Models section of this page explains how this works!

Biological Reactions within the CSTR

As mentioned within Wet Lab modeling, there will be multiple reactions occurring inside the bioreactor, some leading to the production of Ursolic Acid, and some that create other metabolic by-products. For these reactions, we can follow a general structure for developing a model:

1. Identifying the key biological reactions and pathways that occur in the bioreactor, such as the production of UA, as well as other reactions that could possibly compete with that production.
2. Defining the reaction kinetics for these biological reactions, as well as for the growth of yeast using Michaelis-Menten and Monod models.
3. Develop mass balance equations given the information from the kinetic models.
4. Integrate or account for environmental factors into the models created.

Before developing a model, it is important to state some assumptions that will have an impact on the results of the model:

1. The reactor contents are well mixed, so that the concentrations of species, yeast, and environmental conditions are constant throughout the CSTR
2. A constant temperature, pH, and dissolved oxygen for aerobic metabolism and other pathways to proceed normally: optimal conditions will be discussed in a later section
3. The reactor is working under steady-state conditions so that there is no accumulation of any given species
4. Substrate and product inhibition is negligible, which is important especially because we aim to incorporate an excess concentration of glucose in the CSTR
5. The biomass yield coefficients for each species are constants, so that there is a constant and stable relationship between substrate consumption and product formation

The first step in developing models for the bioreactor is identifying the key biological reaction and pathways that occur in the bioreactor. These have already been defined in earlier sections of the document. As a summary, the most important pathway to model out is the actual ursolic acid production pathway, followed by some key competing pathways such as for the ergosterol pathway through the production of lanosterol. Additionally, even the metabolism for yeast will impact the rates of UA production, because in addition to glucose metabolism, lipid synthesis/oxidation will also have an impact because acetyl-CoA is also integrated into the pathway. For these reactions, the general growth of the yeast can thus be modeled for simplicity, as the general metabolic pathways directly relate to the growth of the yeast in the controlled CSTR.

Given these reactions, we can implement Michaelis-Menten kinetics using the same motivation described in the section titled “Substrate-Level Michaelis-Menten Kinetics”, subsection titled “Step 1.”:

Ursolic Acid Production:

• dU/dt refers to the rate at which ursolic acid is being produced
• K_m,U is the Michaelis-Menten constant for UA production
• [S] refers to the concentration of the substrate, which in this case is glucose
• X is the yeast biomass concentration in the bioreactor, which also depends on the growth of the yeast itself

Lanosterol Pathway:

The specific step that is important to account for this model is that it is governed by the enzyme Lanosterol Synthase, which converts 2,3-oxidosqualene to lanosterol, which happens to be the final step in the pathway while also directly competing with the UA production pathway. The latter is necessary for yeast survival, thus ruling out the possibility for knocking out this competing pathway, and so it must be accounted for as it uses up an intermediate used to synthesize ursolic acid.

• dL/dt refers to the rate at which lanosterol is being produced
• K_m,L is the Michaelis-Menten constant for lanesterol production (11 for lanosterol synthase)

Growth of Yeast:

Using the Monod equation, the growth rate of yeast can be modeled to be:

• Mu refers to the growth rate of the yeast
• Mu_max refers to the maximum specific growth rate of yeast, which will have to be determined by testing on our specific bioreactor
• K_s is the Monod constant, referring to the substrate concentration when Mu is at half of Mu_max

Now that the kinetics of these key reactions have been defined, they can be integrated into mass balance equations to determine the change in concentration of each species in the CSTR. A relationship between concentration of yeast and our product, Ursolic Acid can be graphed as well. These will generally take the form:

• dC/dt refers to the change in concentration of the species in the bioreactor
• dP/dt refers to the rate of production of the species in the bioreactor
• D is the dilution rate of the bioreactor, which can be further defined as Fout/V, where Fout is the outflow rate and V is the volume of the bioreactor

Figure 28. Graph of model equation for Ursolic Acid production, modeling the relationship between concentration of yeast and concentration of ursolic acid.

Growth of Yeast:

The change in yeast concentration, because yeast is continuously being fed into the CSTR, must take into account three different factors: the growth of the yeast, as described by the Monod equation, the inflow and outflow of yeast to and from the CSTR, as well as the death rate of the yeast that is in the reactor:

• Fin*Xin is the flow rate from the yeast reservoir into the main CSTR multiplied by the concentration of yeast in that reservoir.
• kd is the specific death rate of the yeast.

• µmax is the maximum growth rate, which is found to be 0.65 hr^-1
• K is found to be 0.108 g/L
• S is the concentration of substrate (glucose) assumed to be at 20.0g/L (2% wt)
• Y is the yield coefficient arbitrarily set to 0.5 (reasonable assumption)
• D is th edilution rate adjusted down to 0.2hr^-1
• X is the concentration of yeast with units g/L

Figure 29. Graph of Biomass and Substrate Concentrations over Time. As time passes, the biomass concentration increases as substrate concentration decreases in the CSTR, and when the substrate concentration has reached a negligible point than growth of biomass stops. After industrial scale-up, the parameters of the model will change and the model will have to be readjusted.

The reactor contents are well mixed

This assumption was made assuming that the agitator in the CSTR was able to create a uniform concentration of species across the reactor itself. However, the model itself does not completely account for that. Without perfect mixing, concentration gradients could develop, as well as unaccountable differences in yields, which could be compounded by the subsequent creation of regions with suboptimal conditions for UA production. This could result in decreased yield and performance inconsistencies.

A constant temperature, pH, and dissolved oxygen

Having a constant temperature, pH, and dissolved oxygen are important because they determine the enzyme kinetics and the growth rate of the yeast itself. Oxygen is especially important as well for maintaining aerobic metabolism of the yeast—as compared to fermentation—maximizing energy efficiency. This is not accounted for in the model, and deviations from these optimal conditions can decrease yield, favor competing pathways, and also decrease the UA yield itself.

The reactor is working under steady-state conditions

The reactor is not at steady-state, there are many other factors during startup and shutdown that could significantly affect the production of species and the growth of the yeast, essentially independent of the equations derived in this section. The yield and other parameters affected by growth and production will be altered—likely be lowered—which can even affect other times of optimal production.

Substrate and product inhibition is negligible

The assumption held that substrate and product concentration remained within optimal ranges, avoiding inhibition effects on both product formation and yeast growth. Although this is not explicitly accounted for in the model, effects such as having too much glucose may suppress UA synthesis by promoting ethanol production through fermentation or increasing the lanosterol pathway etc. However, concentrations for products can be managed through shifting the specifications of the reactor itself.

The biomass yield coefficients for each species are constants

If the yield is variable, then the terms for each of the model equations will be complicated significantly, which could impact the behavior of the yeast in a way that is not accounted for by the model, ultimately leading to a less predictable—and most likely lower—yield.

Modeling Ursolic Acid Detection within the Biosensor

For the biosensor, we aim to use Ohm's law and basic biosensor principles to quantify the amount of ursolic acid. There are three basic components for any electrochemical biosensor: an analyte (what you are measuring), a biorecognition element that binds the analyte, and a transducer that relates the binding of the analyte into a measurable change in voltage. Our analyte, in this case, is of course UA. For our biorecognition element, Ding et al [5]. found that α-glucosidase is inhibited by UA. Our idea is to affix α-glucosidase to a column via an azide click reaction [6] for affinity chromatography. The electrode with the α-glucosidase will use a glassy carbon electrode (GCE) with gold nanoparticles (GNPs) as a transducing element [7].

Figure 29: Overview of Biosensor Components

Our model will affect our development cycle by first providing us with rough estimations of required important quantities like glucose and media given the volume of our main bioreactor tank and the intended extraction rate. We can generate an initial estimate to help us with the initial setup of our pumps to ensure that the reactor never empties and maintains a constant volume during preliminary tests. We can further adjust our model to fit once a working bioreactor is built by making tight adjustments within our model and going through multiple iterations to fully optimize the reactor's conditions for efficient production of Ursolic Acid.

Mathematical Signal Modeling

To measure the efficacy of our sensor, we can develop a mathematical response model. To start, we consider a modification of the Hill equation suited for our scenario of inhibition (the modifications being the replacement of the dissociation constant Kd with IC50, the half-maximal inhibitory concentration):

Where %Response represents the maximum biosensor response (unitless), Rmax represents the maximum response, [G] represents the concentration of glucose, n is the Hill coefficient, assumed to be one for simplicity, and IC50representing the half-maximal inhibitory concentration of UA, found to be approximately 5.08µg. Because we are assuming that the inhibition of α-glucosidase is positively correlated to the amount of UA present in solution, we convert our %Response to a positive signal from the sensor:

We are, of course, making several assumptions to get to this point:

• Ursolic Acid acts as a non-competitive inhibitor, and its binding is reversible.
• Ursolic Acid is the only inhibitor, competitive or non-competitive, present in solution.
• There are no activators of α-glucosidase present in solution.
• The only substrate binding to α-glucosidase is UA.
• α-glucosidase is not affected by other conditions (e.g. temperature, pH) of the outflow.
• A Hill Coefficient of 1, which represents the degree of cooperativity, makes the assumption that affinity is independent of concurrent substrate binding.

\(R_\text{max}\) introduces a difficult problem in this scenario, however, because it is usually calculated experimentally. As such, one of our future steps would be to calculate Rmax so that we can develop an accurate response curve for our model. We would do this using Surface Plasmon Resonance (SPR), but were unable to do so for this project due to time and budgetary restraints. Surface plasmon resonance is an optical phenomenon that occurs when light strikes a metal surface at a specific angle, causing electrons to oscillate and create a rapidly diminishing wave that is highly sensitive to changes in the refractive index near the surface, allowing for label-free, real-time detection of molecular interactions.

Glucosidase Binding Affinity for Ursolic Acid Detection

However, there are other ways to model the biosensor, starting with the biosensor’s biorecognition element. To make our sensor return a measurable signal, we would utilize the noncompetitive characteristics of UA on glucosidase [5]. During start of the reactor, a glucose concentration baseline, using a known glucose concentration that is inputted into the reactor, can be obtained and used to normalize sensor readings once UA starts being produced. In our model, we assume that all glucose binding sites are filled and remain bound to α-glucosidase while UA binds to separate binding sites. This conformational change caused by the UA binding would create a signal between the reference and counting electrodes that can be read and evaluated against the baseline to give accurate UA concentration readings. This would be accomplished with the use of commercially-available glucose sensors.

Wang et al. determined the binding affinity of UA and α-glucosidase at three different binding sites on α-glucosidase with molecular docking scores of −43.938 kcal mol−1, −31.2816 kcal mol−1, and −12.3694 kcal mol−1. These negative scores indicate a strong binding affinity between the two. This can also be visually seen using MolModa, a molecular docking model that we used to confirm the binding of UA and α-glucosidase that Wang et al. reported.

Figure 30: UA Bound to α-Glucosidase

Figure 31: Closer View of UA Bound to α-Glucosidase

These binding affinity scores become useful when considering the theoretical operation of our biosensor. Wang et al. used a linearized Michaelis-Menton equation (shown below) to do so, and we would do the same.

Where v is the reaction velocity in [mol/time], k_m is the dissociation coefficient, k_i is the IC50 coefficient, [I] is the inhibitor concentration (UA), [S] is the substrate concentration (glucose), and v_max is the maximum reaction velocity. This would give an expected binding curve as such:

Figure 32: Tabulation of Data for Figure 31

Figure 33: Expected Plot Based on Wang et. al

What this equation tells us is that the presence of an inhibitor increases 1/v thereby decreasing, through an inverse proportionality, the activity of the enzyme. With this and our previous discussion on response percentage, we can construct an approximation for the change in voltage we expect to see with the binding of Ursolic Acid.

Figure 34: Expected Voltage Strength Trend

These models are valid with the following assumptions:

1. The Hill Coefficient remains at 1, meaning binding UA does not change the affinity of additional UA being bound. This is a common assumption to simplify the math by ignoring transient change in bound inhibitors.
2. The electrode's voltage is directly proportional to signal strength. This was also assumed to simplify the math, as with a Hill Coefficient of 1 we are assuming that the binding will only increase linearly.
3. The voltage would act in the same direction of existing current (not opposing flow). A voltage acting in opposition to the existing current would accomplish the same thing but would add unnecessary complexity to the calculations.
4. All enzymes already have bound glucose when measurements begin. This simplifies the calculations by removing the extraneous “priming” time where we are waiting for glucose to bind.
5. Glucose concentration ([S]) in the outflow is constant. Our reactor is supposed to keep concentrations in the main chamber constant and given that we are already using the well-stirred assumption in our bioreactor modeling, the outflow should also have a relatively constant concentration.

We now have a method to detect Ursolic Acid concentration in the biosensor but there are a few factors in vitro that might affect the efficacy of our model. The Hill Equation with n=1 provides a useful starting point for developing a basic model of ursolic acid's interaction with α-glucosidase. However, experimental data is crucial to confirm whether α-glucosidase exhibits cooperative binding or not, which would influence the accuracy of this model. Furthermore, molecular docking using MolModa indicates that ursolic acid has a high affinity for α-glucosidase, suggesting that the sensor can respond to multiple binding conformations. This versatility in binding increases the likelihood of both true and false positives, particularly if molecules structurally similar to ursolic acid bind to α-glucosidase. Finally, the linearized Michaelis-Menten plot presented by Wang et al. offers valuable insights into the relationships we should anticipate and examine during experimental testing. Specifically, we expect to see the non-competitive inhibition of α-glucosidase by UA and a corresponding linear decrease in voltage from the electrode. Finally, we expect to see a linear change in circuit current as a result of the voltage change as the former is linearly proportional to the latter via Ohm’s Law. In the real world, we expect these trends to be affected by factors in the outflow that would affect the activity of α-glucosidase such as pH and temperature, likely affecting our readings but not so much as to make the trend nonlinear.

These relationships can serve as important benchmarks for both validating our model and understanding the kinetics of the ursolic acid and α-glucosidase interaction. What our model has told us is that there exists a linear relationship between α-glucosidase activity and UA concentration and that we can quantify this relationship by reading a proportional change in voltage within the biosensor.

Scale-Up and Future Considerations

Considering that we are an industrial-scale project, our modeling would not be complete without modeling what the industrial scale of the bioreactor would look like. Luckily, we calculated two dimensionless numbers to help with this in our previous models: Biot Number and Reynolds Number. When scaling up a reactor from lab-scale to industry, the dimensionless numbers must remain constant to ensure the model’s physical phenomena reflect that of the scale-up. This is because dimensionless numbers are inherent correlations of the physical conditions of a system. Since they have no dimension, they can be used to describe the same system at any size as long as the other parameters within the system are scaled to keep the overall dimensionless number constant.

Scaling up Heat Transfer Conditions:

The Biot number is the dimensionless number that deals with heat transfer by giving the ratio of convection to conduction. [3] In the “Heat Transfer within the CSTR” "section of modeling, this number was calculated to be 22.7. However, this number was calculated using a thick glass material. While glass can work for lab-scale, it is not a viable material in industry as the industry standard uses a type of steel. However, switching material causes the thermal conduction and convective heat transfer coefficients to change; therefore, we constructed a 3D plot on MATLAB comparing reactor height (L), to the material index (h/k), and the Biot number (Bi).

Figure 35: 3D Plot of Biot Number vs Reactor Height and Material Index (h/k)

An additional code was also produced to input the h and k values of 3 industry standard steels (stainless steel, 6L-AXN stainless steel alloy, and carbon steel 1.5%) and the original values of our model (glass). Inputting this data, the MATLAB code was set to input our preferred Biot number and scan the 3D plot to find the ideal material and the ideal length that would keep the Biot number the same. The results matched to the AL-6XN stainless steel alloy material with a resulting Biot number of ~22.7 as shown below.

Figure 36: Resulting Biot Number of Industrial Scaled Model

By building the industrial bioreactor with 6L-AXN stainless steel alloy, the heat transfer principles will match those of our models and therefore achieve accurate scale-up on the heat transfer side. For other iGEM teams to model their own scale-up of their bioreactor using Biot number, the MATLAB code is linked in a repository on the bottom of the modeling page!

Now that we have the material and height of our industry-scale bioreactor, we can use it to calculate the radius as well. Assuming a cylindrical geometry, an industrial scale of 1000 liters, and having a height of 15.54 meters from Biot number calculations, a simple V=r2h can find the radius. Inputting our numbers, the radius is 4.52 m as shown below. With the exterior scaled, we shifted our focus to what was flowing in and out of the rector.

Scaling up Fluid Flow Conditions:

The Stirred Tank Reynolds number is the dimensionless number that deals with fluid flow by comparing fluid properties within the reactor to determine if the system has turbulent or laminar flow. It is also used to compare a fluid’s inertial force to its viscous force [2]. Since the fluid within the lab scale and industrial scale of the reactor will remain the same, the scaling up of the reactor fluid flow-wise is simple; only the impeller speed (Ni)and diameter (Di) need to be altered. Through a ratio, we can determine what the propeller size of the 9.05 m wide reactor should be using numbers from our Visimix model. Through the ratio, the propeller diameter should be 6.46 m as shown below.

Using that diameter, it can be inputted into the Reynolds number equation to determine the propeller speed and maintain the same fluid flow conditions. If Reynolds number were to be changed to be less turbulent, these calculations can be easily repeated to adjust the propeller diameter and speed and meet the new conditions. Using the equation, the propeller speed necessary to scale-up fluid flow conditions is 0.001 rad/s as shown below.

However, this is incredibly slow and is not industrially feasible. Due to this, we decided to shift our industrial scale model to change from 1 large propeller to many small ones instead. This way, we can keep the propeller size small in the industrial scale but raise the propeller speed to numbers typically seen in industry. Instead of using diameter to calculate speed, we did the inverse and used speed to calculate propeller diameter.

By doubling the speed to 200 RPM (20.94 rad/s), the propeller speed meets average industry standards [4] with a propeller diameter of 0.047 meters. While the propeller diameter is small, adding multiple of them throughout the height of the reactor will allow for well circulation of the media as it goes through the bioreactor. With this alteration, the fluid flow of the reactor can be successfully scaled-up.

Putting it Together and Future Steps:

Through scaling up the heat transfer and fluid flow conditions of the lab-scale reactor, we were able to design an industrially-viable bioreactor to license to pharmaceutical companies. A visual representation of the scaled bioreactor, with the parameters necessary to match our model, is shown below.

Figure 37: Industrial-Scale Oncurex Bioreactor Design (Not to Scale)

With the new scale, alterations would have to be made on how to maintain conditions within the reactor. Rather than a sous-vide machine and water bath to maintain internal reactor temperature, a heating and cooling jacket can surround the bioreactor. Liquid such as glycol would flow within the jacket to maintain temperature control. To uphold safety within the bioreactor, a rupture disk would need to be added to prevent the reactor from over-pressurizing as shown on the top right of the reactor schematic.

While an industrial-scale model was created, it is based on the assumptions that were used when calculating the Reynolds and Biot numbers. Similar to all models, reality can look very different from ideal conditions and is affected by other factors, such as uneven flow patterns. As testing continues for the bioreactor, the model can be updated for the lab-scale and in turn will update the industrial-scale numbers as well.