Team:Virginia/Model

Modeling

MANIFOLD

Index:

Modeling

Overview

Our scaffold-incorporated bacterial microcompartment (BMC) system, called Manifold, incorporates targeted enzyme fusions to the inner BMC shell proteins to spatially orient select enzymes to minimize metabolic flux leakage. For the Summer of 2020, the team thoroughly developed a proof-of-concept for the production of model molecule trans-resveratrol for its compact size and compatibility.

With the onset of the global pandemic and need for quick project adaptations to a virtual environment, the modeling team sought to put a positive twist on the unprecedented circumstances and took the opportunity to flesh out a thorough multi-scale model to support experimental design in wet lab and explore various other facets of resveratrol metabolism optimization. We utilized a multitude of mathematical models in both Matlab and Python to understand the viability of our BMC system. We chose to embrace the complexities that modeling could envelop and built a multi-scale model balancing assumptions with the realistic conclusions we were to make.

From our endeavors, we determined the optimized expression of our 14 different parts utilizing cataloged parts from the iGEM Parts Registry and the expected increase of biosynthesized trans-resveratrol when utilizing the Manifold system. Our silico chronology began with a considerable amount of literature review to familiarize ourselves with recent advances in BMC research. From this, we were able to construct an optimization model for ideal part expression used in future wet lab practices. We then delved into the reaction kinetics of our system to quantify the improvement of Manifold on current methods of trans-resveratrol biosynthesis.

The modeling pipeline for our project was driven by two main questions:

In what relative concentrations should the twelve parts of our system be expressed?
How much can a fully assembled Manifold system be expected to increase resveratrol production when compared to a control system?

I. Research-Driven Analysis

The in silico investigation of Manifold was guided by these questions relating to the improved assembly and efficiency of our system. In order for our design to present a significant breakthrough in metabolic engineering, we needed quantitative answers. To first construct these questions, the modeling team conducted research driven analysis to determine what parts of our system may already be characterized and which may require further modeling efforts.

The bacterial microcompartment shell became the first topic of conversation. The [21-gene] Pdu operon (named for its natural 1,2-propanediol metabolism pathway) derived from Citrobacter freundii has been only recently elucidated in literature. The shell is known to be composed of certain Pdu proteins with Pdu enzymes naturally localized to the lumen (or interior) of the compartment. The proteinaceous shell ideally forms an icosahedral shell with a mean face-to-face diameter of 117.99 nanometers (Figure 1)

[1] Kennedy et al., “Apparent Size and Morphology of Bacterial Microcompartments Varies with Technique.” 2020.

. Our study utilized seven Pdu proteins (PduABJKNUT) which are necessary for the formation of this nanocarrier and therefore lead to BMCs of the same size.

Fig 1. 3D rendering of the icosahedral shape of an assembled BMC. Dimensions formulated from findings in Kennedy et. al. to be 117.99 nanometers with a side length of 68.12 nanometers

[1] Kennedy et al., “Apparent Size and Morphology of Bacterial Microcompartments Varies with Technique.” 2020.

. Note: Images represent the idealized BMC structure based on hypothesized methods of Pdu protein assembly.

Our selection of this synthetic Pdu operon was bolstered by a recent study elucidating the stoichiometric composition of these proteinaceous enclosures

[2] Parsons et al., “Synthesis of Empty Bacterial Microcompartments, Directed Organelle Protein Incorporation, and Evidence of Filament-Associated Organelle Movement.” 2010.

,

[3]Yang et al., “Decoding the Stoichiometric Composition and Organisation of Bacterial Metabolosomes.” 2020.

. This resource quantified the ratio between Pdu proteins apparent in BMC enclosures (Figure 2) which provided a solid quantitative basis for our desired parts analysis. With twelve parts required for full assembly of the Manifold system, the challenge of correct, relative expression became a pressing issue. With the realization that wet lab data was not forthcoming to assist in testing protein copy number through fluorescence or other microscopy methods, our attention shifted towards expanding upon methods used in literature to create an optimization model in silico (Part II).

Fig 2. Pdu BMC shell schematic representing all constitutive shell proteins

[3] Yang et al., “Decoding the Stoichiometric Composition and Organisation of Bacterial Metabolosomes.” 2020.

. Note: Model based on absolute quantification data from Liu et. al. (2020) but protein arrangement is still only speculative.

II. Relative Part Expression

BMC Background

The PduJ protein is the most abundant protein (5212 monomers) of the BMC shell

[3]Yang et al., “Decoding the Stoichiometric Composition and Organisation of Bacterial Metabolosomes.” 2020.

. Therefore, any Manifold system with correct construction of the BMC would mandate the expression of at least 5212 PduJ proteins. We therefore used the required copies of PduJ as a relative lower bound for computing idealized expression of all other parts of the Manifold system. These predicted stoichiometries in part expression were deduced from theoretical relationships within the bacterial microcompartment.

One other key stoichiometric relationship is mediated through the PduA protein. The N-terminal of the PduD protein is hypothesized to bind with the C-terminal (luminal-facing) residues of PduA

[4]Jorda et al., “Exploring Bacterial Organelle Interactomes.” 2015.

. Manifold utilizes the first eighteen amino acids of PduD as a tag for lumen localization to the BMC as has been performed before

[5]Fan and Bobik, “The N-Terminal Region of the Medium Subunit (PduD) Packages Adenosylcobalamin-Dependent Diol Dehydratase (PduCDE) into the Pdu Microcompartment.”2011.

. This peptide sequence was combined with the zinc finger domain to enable DNA scaffold localization to the lumen of the BMC for subsequent enzyme attachment. This connectivity is significant in establishing the quantified relationships between Manifold components which are further explored in the part expression model below.

Model Activation

Correct stoichiometry between Manifold parts required a method of quantification for both transcriptional and translational parameters. Due to the large number of required components, part expression would have to be carefully modulated to be in high enough quantity for a high-output system without killing the E. coli host. Transcriptional rate relativity has already been quantified through the Anderson Promoter Library. The Anderson RBS Library has introduced several RBS combinations but the translational rates of these sequences has yet to be quantitatively compared. We then looked at the RBS Calculator by the Salis Lab. This software calculates the translation initiation rate (TIR) based on a thermodynamic model of free energy on and around the start codon

[6]Salis, Mirsky, and Voigt, “Automated Design of Synthetic Ribosome Binding Sites to Control Protein Expression.” 2009.

. Using 21 Anderson RBS sequences and two others found naturally in acquired vectors, we found the calculated TIR values for our eight parts requiring an RBS. These values were used in conjunction with the relative transcriptional values of 19 Anderson Promoters to find optimized sequence combinations to achieve correct relative part expression. Due to the sheer number of combinations, the process yielded computational ratios with minimal error compared to ideal, theoretical ratios (Figure 3). Important assumptions for this model are listed below.

Assumptions

Assume steady state (allows for transcriptional and translational rates to be multiplied)
Assume ACC and HIV RT protein copies are limited by subunit with lowest translation initiation rate
Assume ribosome binding only occurs at primary start codon less than 20 base pairs downstream of selected RBS sequence

Part	# of molecules	Ideal Ratio	RBS	Promoter	Computational Ratio	Error
PduD fusion	1841	0.353	BBa_J61107	BBa_J23101	0.357	0.00335
GFP fusion	3682	0.706	BBa_J61111	BBa_J23106	0.711	0.00500
PduJAKBNUT	5212	1.000	Warren RBS	BBa_J23114	1.000	0

Fig 3. Results from the relative expression model which show ideal and computed ratios for individual part expression relative to PduJ. Top: The bar graph represents the differences between ideal and computationally-derived ratios for each part. Ideal ratios are derived from stoichiometric comparisons between parts in the idealized system derived from literature. Computational ratios are derived from the model output using the product of Anderson promoter transcriptional rates and translation initiation rates of Anderson RBS sequences. Bottom: The table presents sample calculations for finding both ideal and computational ratios along with the corresponding error. Note: All r_oligo parts are uniformly represented as one value and PduA’BJKNUT is omitted for simplicity.

Exploring Model Validity

The relative expression model detailed the success in customizing Anderson Promoter and RBS sequences for individual part expression relative to PduJ. This success is characterized by the minor deviations between the ideal and computed ratios. A supplementary model was created in complement to the relative expression model to verify its findings. Specifically, the central dogma theory and the creation of our DNA Scaffolds was investigated. The relative ratios of the human immunodeficiency virus reverse transcriptase compared to PduJ is observably low (Figure 3). Efficient functioning of the Manifold system requires a sufficient amount of DNA Scaffolds to be produced and encapsulated by the BMC before the compartment becomes fully enclosed. This led us to further analysis to determine the applicabilty of this computational-derived ratio for HIVRT expression in establishing correct assembly of our system.

Model Constraints

Research in BMC assembly remains primarily speculative. Therefore, we used information on viral capsid assembly (a similar proteinaceous compartment) for timeline comparison. Analysis of MS2 bacteriophage viral capsid assembly details the protein shell formation required around 280 seconds for nucleation, growth, and lag time

[7] Garmann, Goldfain, and Manoharan, “Measurements of the Self-Assembly Kinetics of Individual Viral Capsids around Their RNA Genome.” 2009.

. These assembled particles are around thirty nanometers in diameter upon completion. Assuming assembly time only increases with particle size, we used this timestamp as a lower bound for our BMC formation. In addition to the expected assembly time, it is important to know the amount of DNA Scaffold we may expect to be encapsulated within the BMC upon formation. The PduD localization peptide has a measured Incorporation Efficiency of 24.7% to the lumen of the BMC

[8] Juodeikis et al., “Effect of Metabolosome Encapsulation Peptides on Enzyme Activity, Coaggregation, Incorporation, and Bacterial Microcompartment Formation.” 2020.

. Using these parameters, our system would require the presence of 37,267 (or 4.571 on the log10 scale) molecules of DNA Scaffolds before t = 280 seconds.

We defined our relationships according to the schematic below (Figure 4). This process develops the connections between the input transcriptional/translational controlling factors (Promoter and RBS) and the output number of DNA Scaffold molecules. Our model used a numerical solver of the differential equations below to test this input-output relationship (Table I). Our assumptions are also listed below.

Fig 4. Central Dogma Schematic. Where ORI represents the plasmid copy number; kr represents the transcription rate; kd represents the degradation rate; kl represents the translation rate; krHIVRT represents the HIV-RT transcription rate; and ka represents the annealing rate.

Model Equations

Table I: Differential Equations representing components of the central dogma theory.
Component	Differential Equation
HIV-RT mRNA	d(HIV RT\|mRNA) / dt = k_{rHIV RT}[HIV gene] - k_{dHIV RT\|mRNA}[HIV RT\|mRNA]
HIV-RT	d(HIV RT) / dt = k_l[HIV RT\|mRNA] - k_{dHIV RT}[HIV RT]
r oligo	d(r_oligo ) / dt = k_{rHIV RT\|HIV RT}[r_oligogene] - k_{dr_oligo} [r_oligo] - ( k_a [r_oligo] / AV_cell ) ²
DNA Scaffold	d( DNA_scaffold ) / dt = ( k_a [r_oligo] / AV_cell ) ² - k_dDNA [ DNA_scaffold ]

Assumptions

Assume MLRT expedites process and is therefore not significant in calculating lower bound
Assume optimized IPTG concentration for part expression and E. coli survivability

Fig 5. The kinetics involved in DNA Scaffold formation in the Manifold system. Molecules of HIV RT mRNA represented by light blue curve; HIV RT protein represented by dark blue curve; r oligo strand represented by red curve; DNA Scaffold represented by magenta curve. The black circle indicates the point on the DNA Scaffold curve at t = 104 seconds which corresponds to 38,244 molecules.

Model Conclusions

The results of the model effectively predict the presence of 38,244 molecules of DNA Scaffolds after 104 seconds in the E. coli host (visualized by black circle in Figure 5). There are sufficient copies of DNA Scaffold available before the projected time of full encapsulation. This leads us to assume the promoter/RBS combinations for parts required for formation of DNA Scaffolds demonstrate ideal part expression for proper DNA Scaffold inner-luminal binding given our parameters.

III. Resveratrol Production

How much can a fully assembled Manifold system be expected to increase resveratrol production when compared to a control system?

This question became the root of our modeling endeavours. Controlling a limited reaction space inferred improvement in reaction kinetics but drawbacks in molecular availability within the microcompartment. We generated a mass action kinetics model to understand the true viability of compartmentalizing the reactions responsible for resveratrol synthesis.

Fig 6.Trans-resveratrol metabolic pathway, a two-step process, specifying substrate, product, and key enzymes.

Michaelis-Menten Reaction Kinetics

Enzyme-catalyzed reaction equations are commonly derived through Michaelis-Menten kinetics. Using these equations relies on steady-state assumptions that substrate concentration is much greater than enzyme concentration and enzyme-substrate complexes are reached rapidly. Resveratrol biosynthesis involves reactions with multiple substrates and competitive inhibition (Figure 6). The single-substrate Michaelis-Menten equation along with its derivatives are shown below.

Michaelis-Menten Single Substrate

V =
k_cat [E]_T [S] / k_m + [S]

Michaelis-Menten Two-Substrate

V =
k_cat [E]_T [A] [B] /
1 / k_A k_B
+ ^[A]⁄_{k_B} + ^[B]⁄_{k_A} + [A] [B]

Michaelis-Menten Two-Substrate Competitive Inhibition

V =
k_cat [E]_T k_A k_B [A] [B] [I] / k_I + k_I k_A [A] + k_I k_B [B] + [I] + k_A k_B [A] [B] [I]

Assumptions

Assume enzyme-substrate complex reached rapidly
Assume ATP and bicarbonate in excess
Assume reactions are all irreversible
Assume binding is independent for all substrates in all reactions (n=1)
Assume competitive inhibition not significant in STS reaction (Ki >> Km)

Table II: Michaelis-Menten derived Differential Equations for Resveratrol Synthesis
Enzyme	Equation	Reaction Type
ACS	V_ACS = k_cACS [ACS]k_m\|CoA k_m\|ac [CoA][ac][a-CoA] / k_I\|ACS +k_I\|ACS k_m\|CoA [CoA] +k_I\|ACS k_m\|ac [ac]+[a-CoA] +k_m\|CoA k_m\|ac [CoA][ac][a-CoA]	Competitive Inhibition Multi-Substrate
ACC	V_ACC = k_cACC [ACC][acetyl-CoA] / k_m\|acCoA + [acetyl-CoA]	Single Substrate
4CL	V_4CL = k_c4CL [4CL] [p-coumaric acid] [CoA] / 1 / ( k_m\|pco ) (k_m\|CoA4CL) + [p-coumaric acid] / k_m\|CoA4CL + [CoA] / k_m\|pco + [p-coumaric acid] [CoA]	Multi-Substrate
STS	V_STS = k_cSTS [STS] [mal] [4-co] [a-CoA] / 1 / ( k_m\|4-co ) (k_m\|mal) + [mal] / k_m\|4-co + [4-co] / k_m\|mal + [mal] [4-co]	Multi-Substrate

Control Model Conception

We used the single substrate Michaelis Menten equation and derivative formulas in our system reactions. The Acetyl-CoA Synthetase (ACS) reaction requires a competitive inhibition multi-substrate equation due to the competitive binding of Acetyl-CoA in the active site. Acetyl-CoA Carboxylase (ACC) proceeds as a single substrate reaction due to assumed excess of bicarbonate ion and ATP. 4-Coumarate-CoA ligase (4CL) and Stilbene Synthase (STS) enzyme reactions are both characterized by multi-substrate binding. These reaction kinetics equations are described below (Table I) along with the necessary assumptions for characterization.

We used a numerical solver in MATLAB to display the change in concentration of metabolites responsible for resveratrol production in a free enzyme E. coli culture (Figure 7). Concentrations of acetate were assumed to not change significantly due to glucose supplied by the medium. Likewise, the concentration of Coenzyme A was assumed to not change due to 4:4 stoichiometric recycling within the reactions. Our first model demonstrates a control reaction. This is accomplished through examining overall metabolite levels per cell without introducing the BMC. As seen in Figure 7, we looked at the overall metabolite concentrations in an individual cell when engineered to express 18,410 molecules of each enzyme. We then measured the concentration of resveratrol based on an established cell density of E. coli cells per liter

[9]Glazyrina J, Materne EM, Dreher T, Storm D, Junne S, Adams T, Greller G, Neubauer P. High cell density cultivation and recombinant protein production with Escherichia coli in a rocking-motion-type bioreactor. Microb Cell Fact. 2010 May 30 9: 42. doi: 10.1186/1475-2859-9-42. p.4 left column 5th paragraph

. This system utilizes the equations from Table II and initial concentrations based on the E. coli cell volume which may be found here.

Fig 7. A: Resveratrol production over 12 hours in a free enzyme E. coli system. The overall metabolite concentrations in an E. coli cell are observed over 12 hours. B:The yield of resveratrol in bulk is viewed over 12 hours based on a total cell density of 4.8 ⋅ 10¹¹cells/L. We realize the presence of metabolites at concentrations between 10^-20 and 10^-6 mM are not biologically possible (i.e. suggests the productions of fractions of molecules) as seen in the top two figures. It is important to understand the model approximation of these values is theoretical in nature because it best describes the kinetics in a single cell. The average cell resveratrol production is described over this 12 hour period which means some individual cells may be producing more resveratrol or none at all. The estimation of overall cell culture resveratrol production is displayed in image B which uses these individual cell averages to calculate a bulk resveratrol titer in a greater cell population. More mathematically, the curve represents the resveratrol concentration (in mM) multiplied by the cell density (4.8 ⋅ 10¹¹cells) and the molar mass (228.25 g/mol) to examine the resveratrol titer in the bulk sample of cells (in mg/L) over 12 hours.

Model Results

The free enzyme metabolite concentrations display the impacts of starting concentrations on metabolite curves. Acetyl-CoA, malonyl-CoA, and p-coumaric acid are present at much greater concentrations in the initial phases of the reaction and therefore reach a steady-state immediately. The 4-coumaroyl-CoA and resveratrol curves have a pronounced burst phase of approximately 3 minutes before reaching steady state increases.

The major takeaway from this model is the resveratrol titer after 12 hours. Our system concludes a titer of 636.1 mg/L after 12 hours. This result is consistent with previous in vitro studies using the same derived enzymes as our project (4CL enzyme from Arabidopsis thaliana and STS from Arachis hypogea) at similar conditions

[9]Lim et al., “High-Yield Resveratrol Production in Engineered Escherichia Coli.” 2011.

. This consistency with literature demonstrates more confidence in our conclusions when developing comparisons of the control system to the BMC model.

The BMC Pore

The successful modeling of the Manifold system requires two important changes from the control reaction kinetics model. First, the adjustment of enzyme concentrations. Though we propose consistency in the number of overall enzyme molecules present for reaction, the enzyme concentrations through the BMC model will be increased due to a decrease in the reaction space. The other change results from the activity of the BMC pore and the transport of reactants into our system. This requires a deeper look into the pore to understand if, and by how much, the pool of reactants may change.

The BMC Pore is a complex system. Previous works have proposed the mechanism of pore transport is mediated through hydrogen bonding with the PduA Serine-40 residue. This proposition relies on findings in a natural propanediol BMC where molecules such as glycerol and 1,2-propanediol are observed to freely diffuse through the pore

[10]Crowley et al., “Structural Insight into the Mechanisms of Transport across the Salmonella Enterica Pdu Microcompartment Shell.” 2010.

. Both molecules present functional groups which facilitate hydrogen bonding. Further, propionaldehyde has an aldehyde functional group with minimal ability for hydrogen bonding leading to limited pore diffusion.

Manifold Pore Dynamics

For the scope of our project, we were interested in the proposed pore dynamics of our reactant species. Originally, our system was designed to utilize 4CL and STS inside a BMC requiring Malonyl-CoA as the primary reactant. Upon inspection of the pore profile of Malonyl-CoA (Figure 8), we were able to visualize the interactions between Malonyl-CoA and the PduA serine residue. However, additional factors including size and bulky side groups introduce a higher energy profile which discourage pore diffusion. We amended our design and proposed the use of two additional enzymes for better availability of our precursor pool for BMC reactions. Through introduction of ACS and ACC enzymes, our reaction relies on acetate and p-coumaric acid diffusion through the BMC pore in place of the bulky Malonyl-CoA.

Fig 8. Visual of malonyl-coa with observed hydrogen bond interactions in the PduA pore. The protein structure (green) of the PduA residues surround the main chain of a malonyl-CoA molecule (predominantly light blue). Yellow dashed lines represent hydrogen bonds between malonyl-CoA and the serine residue of the PduA shell protein. Note this image presents one possible conformation of the malonyl-CoA molecule in the BMC pore. The image is purposely chaotic to prove the predicted minimal diffusion of malonyl-CoA due to molecule size and other steric hindrances which increase the energy barrier for particle entry.

Acetate molecules are on the length-scale of 1,2-propanediol and propionaldehyde molecules which leads to better conclusions about its pore characteristics. Most importantly, acetate has a carboxyl functional group which has ability for hydrogen bonding. We examined acetate in the BMC pore alongside propionaldehyde and 1,2-propanediol to develop conclusions about its pore dynamics (Figure 9). The obvious difference in hydrogen bonding is easily seen between propionaldehyde (1 hydrogen bond acceptor) and 1,2-propanediol (2 sites for both hydrogen bond donation and acceptance). Acetate also has two sites for hydrogen bonding but relies on hydrogen bond acceptance due to the deprotonation of the acetate ion at physiological pH. P-coumaric acid presents two favorable hydrogen bonding functional groups as well (not shown) and is therefore assumed to present similar pore energies to acetate. Developing quantitative conclusions between these images may be misleading due to the ability of several more molecular arrangements within the pore. However, the visual comparison of acetate’s functional groups with those of 1,2-propanediol and propionaldehyde may lead to claims of relative pore energies. Due to the number of functional groups in the three molecules and their interactions with PduA, we determined a qualitative relationship between the BMC pore permeability (P) of the three molecules as follows:

P_{Propionaldehyde} ≤ P_Acetate & P_{p-coumaric acid} ≤ P_{1,2-Propanediol}

and

P_{Propionaldehyde} ≠ P_{1,2-Propanediol}

Due to the availability of pore parameters for the more commonly studied propionaldehyde and 1,2-propanediol, we used this relationship to infer the acetate and p-coumaric acid energy profile in the BMC pore.

Fig 9. Visual of molecules with observed hydrogen bond interactions in the PduA pore. Molecular representations of propionaldehyde (left), acetate (middle), and 1,2-propanediol (right) with serine residues on the PduA protein. Hydrogen bonds less than 4 angstroms are expressed as yellow dashed lines. Note all molecules were constructed and inserted into the pore along the same plane and at the same point in space. These molecular conformations in the pore represent a sample of possible arrangements.

Pore Transport Model Conception

The main goal of the pore dynamic model is to investigate the mass exchange of natural BMC metabolites between the cytoplasm, pore, and BMC. A steady state concentration of these molecules inside the BMC would inform us of the theoretical reactant pool available in our Manifold system. We investigated the pore dynamics through a three-compartment mass action model. The flux between these compartments was derived from Fick’s First Law:

J = -D ⋅
dC / dz
= A ⋅
dM / dt

Where J is mass flux (kg / s⋅m²
), D is the diffusion coefficient ( m² / s
), A = area (m² ), and M = mass (kg). We also used a diffusion equation for inclusion of an activation energy corresponding to the energy required for pore diffusion:

D = D_maxexp(
E_a / RT
)

Where D is the diffusion coefficient ( m² / s
), D_max = maximum diffusion ( m² / s
), E_a is the energy difference between the pore and cytoplasm/bmc ( J / mol
), R is the gas constant ( J / mol⋅K
)and T is the Temperature. Combination of these two equations led us to our mass transfer equations:

[
dM_CP / dt
= -D ⋅ A ⋅
( ^P_i⁄_{V_P} - ^C_i⁄_{V_C} ) / ^d⁄₂
] and [
dM_PB / dt
= -D ⋅ A ⋅
( ^B_i⁄_{V_B} - ^P_i⁄_{V_P} ) / ^d⁄₂
]

Where P is the pore mass (kg), C is the cytoplasm mass (kg), B is the BMC mass (kg), V_P is the pore volume (m³ ), V_C is the cytoplasm volume (m³ ), and V_B is the BMC volume (m³ ).

We used a numerical solver to evaluate the differential equations with a dt of 50 picoseconds to observe the logarithmic shape describing the concentrations of metabolites in each compartment based on previously observed pore profile energies of propionaldehyde and 1,2-propanediol (Figure 10). Assumptions for this model are listed below:

Assumptions

Assume the concentration at certain point represents homogeneity to the entirety of the compartment
Assume symmetry in the energy profiles on either side of the pore
Assume entry into pore is strictly along z axis with no radial considerations for change in the energy barrier

Fig 10. The change in concentration of propionaldehyde and 1,2-propanediol in the three-compartment system. Curves represent the BMC concentrations of 1,2-propanediol and propionaldehyde based on predetermined pore energies after 7 milliseconds. All constants including energies and other important parameters are found here. Note that due to the computational burden by using a mandatory dt of 50 picoseconds, 7 milliseconds becomes the assumed quasi-steady-state of the system.

Using the initial concentration of 15.542 mM in the cytoplasmic space, 1,2-propanediol and propionaldehyde have observed concentrations of 14.775 and 7.020 mM, respectively after 7 milliseconds. 1-2-propanediol is examined to reach a quasi-steady-state concentration in the BMC with the cytoplasmic concentration after 7 milliseconds. Due to the higher energy profile of propionaldehyde in the pore, we see a well-defined amount of propionaldehyde molecules which lack sufficient energy to diffuse through the pore after this same time period. We used the expected propionaldehyde and 1,2-propanediol concentrations in the BMC (computed as the BMC concentrations after 7 milliseconds) as the lower and upper bounds, respectively, for the normalized initial concentration of compartmentalized acetate and p-coumaric acid.

BMC Model Conception

In order to determine the comparative increase of resveratrol yield in the Manifold system, we had to describe the BMC reaction kinetics as they compared to the free enzyme control. This model would utilize the same equations and constants from Table II with a few important changes. As previously mentioned, the enzyme concentrations were heightened due to the decrease in reaction space between the E. coli cell volume and the BMC volume (approximately a 1000 fold increase). However, the total enzyme molecular count per cell persists between these examples (18410 molecules in free enzyme E. coli cell versus 3682 molecules per BMC with an assumed 5 BMCs in each E. coli cell). Also, as previously explored, the initial concentration of acetate has been normalized as a range between 7.020 and 14.775 mM compared to the free enzyme initial concentration of 15.542 mM. Similarly, p-coumaric initial BMC concentration was normalized to be between 0.1129 and 0.2377 mM compared to the free enzyme initial concentration of 0.25 mM. We used a numerical solver to model the metabolite concentrations over time in E. coli expressing BMCs (Figure 11). Important assumptions are addressed below:

Assumptions

Assume quasi-steady state reached after 7 milliseconds
Assume 5 idealized BMCs created per E. coli cell
Assumptions made in the free enzyme reaction model
Assume prior, in silico derived pore energies for propionaldehyde and 1,2-propanediol are well-representative of true pore energies

Fig 11. A: Resveratrol concentrations in a free enzyme and BMC system over 12 hours. All metabolite concentrations measured per BMC reactor. Note the same computational logic explained in Figure 7 applies to this model.B: The comparative resveratrol yields displayed over 12 hours in a bulk cell population. Blue lines represent the lower and upper bound estimates for resveratrol production based on pore calculated concentrations of reactants p-coumaric acid and acetate at t=0. The green line (B) represents the free enzyme production of resveratrol after 12 hours. The resveratrol yield was computed as the resveratrol concentration over time (in mM) multiplied by the cell density (1.0 ⋅ 10¹¹cells/L) and molecular weight of resveratrol (228.25 g/L) as well as the assumed presence of 5 BMCs per E. coli cell to achieve the resveratrol titer (in mg/L) in the BMC system. Note the cell density was adjusted in accordance to prior work displaying OD600 changes upon BMC protein expression when compared to a control (1.0 ⋅ 10¹¹cells/L compared to 4.8 ⋅ 10¹¹cells/L). Find more information here.

Model Conclusions

The improvement of resveratrol production is easily visible in the comparison between the free enzyme and BMC system. The final resveratrol titer in the BMC kinetic model was 131.5 g/L assuming 1.0 ⋅ 10¹¹cells/L (cell density adjusted due to plasmid expression). This would infer the fold change of total resveratrol titer is between 206.5 and 434.2 times higher when comparing the Manifold system to a free enzyme system.

It is important to note the importance of assumptions when creating statements of this magnitude. Our models are only as conclusive as our assumptions are reliable. Of highest importance are the assumptions of this perfect Manifold system. These include ideal and complete binding of our enzyme molecules in pre-specified ratios, correct absolute initial conditions of precursors in the BMC reaction space, and the homogeneity of our system across tera-magnitudes of cells. Corrections for biological inconsistencies following the successful implementation of our project in wet lab practices (such as observed changes in cell densities) may better inform both reaction kinetic models through parameter validation. The assumption of Michaelis-Menten kinetics should also be reviewed in future work to establish the validity of these equations and the theory of enzyme saturation when enzyme concentrations are noticeably increased. Wet lab inspired corrections may also be applied to the relative part expression model after preliminary wet lab testing for validation of promoter/rbs combos in producing results similar to those found in silico. We realize the estimation for the increase of our compartment yield may be an overestimate due to a flaw in the pore energy found specifically for propionaldehyde. This may lead to an overestimate of approximated propionaldehyde found within the BMC over time seen in Figure 10, leading to an overestimate of the progression of resveratrol yield seen in Figure 11. Unfortunately, much of our final findings depended on this value because we were not provided any other values for comparison. This disparity in available information is something we hope to address in future work. The field behind BMCs is still continuing to grow and we hope the methodologies of our silico practices can aid others in its advancement.

We also want to reiterate our intention in utilizing multiple models. Our modeling pipeline was purposeful in answering our two questions. In order to arrive at realistic conclusions to these questions, we relied on different models to pinpoint certain areas of our complex system. This allowed us to resolve issues (Malonyl-CoA availability) and perform quality control checks on prior models (central dogma validation of the part expression model) along the way. The combination of models therefore minimizes potential for error and unrealistic assumptions for improved confidence in end results. A more detailed diagram explaining the input-output relationships between these models is found below (Figure 12).

Fig 12. Visualization of the Modeling Pipeline and the corresponding input-output relationships between models.

Acknowledgments

We want to thank Ryan Taylor for his commitment to our derivation of the Pore Calculator and continued moral support throughout the project. We also would like to thank Dr. Papin and Dr. Kozminski for project advice and for overall iGEM team advocacy. Other noteworthy contributors to our project included Jainam Modh and Vignesh Valaboju.

Files

All MATLAB files may be found at our github repository. Our parameter table may be found here.

Sources

[1] Kennedy et al., “Apparent Size and Morphology of Bacterial Microcompartments Varies with Technique.” 2020.
[2] Parsons et al., “Synthesis of Empty Bacterial Microcompartments, Directed Organelle Protein Incorporation, and Evidence of Filament-Associated Organelle Movement.” 2010.
[3] Yang et al., “Decoding the Stoichiometric Composition and Organisation of Bacterial Metabolosomes.” 2020.
[4]Jorda et al., “Exploring Bacterial Organelle Interactomes.” 2015.
[5] Fan and Bobik, “The N-Terminal Region of the Medium Subunit (PduD) Packages Adenosylcobalamin-Dependent Diol Dehydratase (PduCDE) into the Pdu Microcompartment.” 2011.
[6] Salis, Mirsky, and Voigt, “Automated Design of Synthetic Ribosome Binding Sites to Control Protein Expression.” 2009.
[7] Garmann, Goldfain, and Manoharan, “Measurements of the Self-Assembly Kinetics of Individual Viral Capsids around Their RNA Genome.” 2009.
[8] Juodeikis et al., “Effect of Metabolosome Encapsulation Peptides on Enzyme Activity, Coaggregation, Incorporation, and Bacterial Microcompartment Formation.” 2020.
[9] Lim et al., “High-Yield Resveratrol Production in Engineered Escherichia Coli.” 2011.
[10] Crowley et al., “Structural Insight into the Mechanisms of Transport across the Salmonella Enterica Pdu Microcompartment Shell.” 2010.