Model

Growth Models

---

Why Modeling Growth is Important for Our Project

Accurate prediction of the growth of our engineered commensal bacteria is essential for the safety and efficacy of our therapy. It is crucial that our bacteria do not overgrow within the gut, as it is important to maintain the natural human microbiome ratios. At the same time, our bacterium strain should not get outcompeted by other bacteria in order to be able to secrete clinically relevant concentrations of the therapeutic drug. Getting accurate information on the growth rate is therefore valuable for our project.

Additionally, we intended from an early point in the project to include biocontainment strategies like kill switches. If our bacteria were to escape from the host into the environment, it could have unknown but disruptive effects on ecosystems. That is why biocontainment strategies are increasingly important in synthetic biology applications.^1,2 However, it is becoming an increasing concern that kill switches produce the toxin even when not induced, which is called leaky expression. This low-level toxin slows growth and leads to evolutionary pressure to lose the kill switch gene.¹ Therefore, it was important for us to measure the potential effect of leaky expression.

Our Approach: Exponential, Logistic and Gompertz Growth Models

To accurately measure growth for the different plasmids, the lab team measured the OD600 & GFP/mCherry production of transformed bacterial plasmids (see Wet Lab: Engineering - cloning, cycle 5&6). To better analyse the experimental results, we developed a MATLAB code which fits 3 of the most common bacteria growth models used: the Exponential growth model, the Logistic growth model, and the Zwietering Modified Gompertz (Gompertz) growth model.

The code can fit the three models to the obtained experimental data and calculate the best-fitting model of the 3. This fitting is performed using MATLAB's built-in solver, lsqcurvefit (nonlinear least squares), which estimates the fitting by minimising the sum of squared residuals between the model predictions and the experimental data. Once the fitting has been completed, it calculates the corrected Akaike Information Criterion (AICc), a metric that evaluates the goodness-of-fit and penalises unnecessary parameters of the models. The combination of these two metrics was especially important because the exponential, logistic and Gompertz models have 1, 2 and 3 free parameters, respectively. Apart from AICc, root mean squared error (RMSE) is another metric used in our post-run analysis to compare how closely our plots compare to the experimental data. Models with a higher number of free parameters should almost always be able to provide a better fit, causing false positive interpretations if only using RMSE to quantify which is best.

All models contain a growth rate, which is a measure of how quickly the bacteria multiply. However, bacteria cannot grow forever in a dish. Eventually, they will run out of space or out of nutrients, and therefore stop growing. So the logistic and Gompertz models include a parameter for max capacity, which approximates the maximum amount of bacteria you can get. The Gompertz model also has another parameter, the lag time. When bacteria are switched from one environment to another, they have to adjust for a while before they grow. This value is modeled by the lag time. These concepts are also shown in Video 1, which uses our self-made explanation software of growth models. The explanation tool can be found in the resources.

Modeling is describing what happens with mathematics. All three of the models we use can be modeled mathematically with a direct formula, as well as a derivative. In our code, we used the direct formula, but we decided to provide both of them on our wiki. The technical and mathematical details can be seen in Table 1. In addition, the table also highlights the limitations and assumptions of the three models.

Table 1: Technical Details of the Exponential, Logistic and Gompertz Growth Models

Build: The Process of Coding the Growth Models

We started by coding the mathematical formulas for the three models, and fitting them to data we made up, but showed growth in a somewhat asymptotic manner. To see if the fitting worked, we plotted the data and models in graphs.

We then started simulating data with the Gompertz formula and a random numerator, and fitted the models against this to have more accurate representations of experimental data.

We then started to code tools that allow us to compare which performs best. The programme calculates RMSE and AIC. To save these results and the fitted parameters, we exported them to an Excel file.

By then, we started to prepare to get data from the lab team, so we built functions that could import data from Excel files. This can be done by selecting the Excel file with the data. Once we got the data, however, realized that we needed to take the plate layout into account. We had multiple repetitions of every condition per plate, so we added a function to the model that would group data based on the condition. For this, you need to provide an Excel file that describes what is seeded per column and per row, and the tool combines this information automatically and groups conditions together. For an example see our raw data files in resources

We then also started to incorporate statistical comparisons between conditions with analysis of variance (ANOVA). Because empty wells were muddling our statistics, we added a section in the code that excludes any wells from rows or columns that are marked as 'Empty'.

At the end of this design cycle, we have a model that can:

• Load Excel files with repeated measurements
• Group wells by condition based on a plate layout
• Fit the growth data from each non-empty well to three growth models
• Record their fitted parameters
• Compare them based on AIC
• Export a table with all relevant data to Excel
• Produce graphs showing the growth and modeled fit for every single well
• Produce graphs for all grouped conditions together, also plotting the mean fit for the condition
• Statistically compare conditions with ANOVA and produce boxplots to visually support the analysis
• statistically compares the mean RMSE, the mean AIC and the mean R-squared of the 3 models.

So now we are ready to apply it to a real experiment for the first time.

Experimental Cycle #1

The goals we had for this experimental cycle were to compare the three models and find out which one was best, so that we could use that for further analysis. Additionally, we were wondering if we could use fluorescent reporters as a proxy for the production of the drug. For the therapy it is important that the drug is produced at a steady rate. However, during the treatment the bacteria will first colonize the gut and grow in number, but then reach max capacity and maintain a steady level. If the number of bacteria does not increase, does the protein production still happen? Therefore, we hoped to see the number of bacteria flatline, and GFP/mCherry fluorescence continue to increase.

So we can summarize the goals as follows:

• Goal 1: Determining the best model
• Goal 2: Determining if bacterial growth is different from the increase in fluorescence intensity
• Goal 3: Determining if fluorescence intensity keeps increasing even when bacterial growth reaches max capacity
• Goal 4: Compare the growth rates for different level 1 plasmids
• Goal 5: Identify a good starting OD

Testing the models and exploration with initial data

The lab team acquired overnight growth data, with 192 measurements of OD600, GFP fluorescence intensity and mCherry fluorescence intensity for a total of 576 measurements. A 96-wells plate was seeded with E. coli that was transformed with geneblocks 1-8 and an untransformed E. coli (see Wet Lab - Engineering Cloning cycle 5&6). For analysis, we excluded geneblock 1-4 however, as sequencing showed that the constructs we ordered were incorrectly synthesized by IDT.

When we ran the code, we could see that the growth experiment was successful and that the models worked as intended (Fig. 1A, C, E). The results showed that the Gompertz model had a better fit in terms of AIC than both the exponential and the logistic growth model for the majority of the experimental data set obtained (Fig. 1B, D, F-H). This was what we expected.

An observation we did make though is that GFP and mCherry fluorescence intensity would reach a plateau long before the end of the experiment. This is because the CLARIOstar has a maximum value of 260000. The gain of the fluorescence intensity was set so that the highest value at the start of the experiment would be 10%. This meant however that the mCherry and GFP signal could only increase 10-fold before reaching the max, meaning they hit it too fast. For GFP this effect was a bit less, because the medium is slightly autofluorescent, meaning that not all of the green signal was GFP. If for example GFP was 50% of the green signal at the start, it would mean it was only 5% of the machine's max value, allowing a 20-fold increase.

Knowing that the Gompertz model worked best we started to use it for all the analysis from now on. So we used it first to see if seeding at different starting OD's changed the growth rate. We expected it would be the same, but we needed to confirm before pooling the data together. The different starting ODs were determined to have no differences with ANOVA (Fig. 2A).

So we pooled the data together and analysed if E. coli with different geneblocks showed different growth. We used geneblock 8 as a control for all analysis, unless mentioned otherwise. We chose geneblock 8 as it was expected to only produce cytosolic mCherry and was part of every single level 2 plasmid we would test. Compared to geneblock 8 none of the transformed bacteria showed different growth rates (Fig. 2B).

However, a post hoc analysis did show that the GFP-producing bacteria (GB 5 and 6) grew slower than the L-DOPA producing bacteria (GB7), for which we have no real explanation. Perhaps the metabolic cost of the GB7 is lower than for GFP, or it is less toxic to the cell when overexpressed but we don't know why this would be the case.

This made us wonder if the increase in GFP or mCherry fluorescent signal would be faster than the growth rate of the actual bacteria. However, it seems that the increase in fluorescent signal corresponds to the growth rate of the bacteria (Fig. 2C). It does show a statistical significance for the GFP to bacteria growth rates, but we think this interpretation is not meaningful due to noise by the green autofluorescent medium that we talked about before, and the extremely low variance in the fitted growth rate for bacteria.

What We Learned from Experimental Cycle #1

• Goal 1: Gompertz model performs the best
• Goal 2: Increase in fluorescent signal and number of bacteria happen at similar rates, or perhaps fluorescence increases slightly more quickly for GFP-producing bacteria
• Goal 3: We were unable to answer. But we learned that we should adjust the gain for fluorescent measurements such that the highest starting value is 2% of the saturation limit of the CLARIOstar instead of 10%
• Goal 4: The L-DOPA producing plasmid seemed to show a higher growth rate than the GFP containing plasmids. However we do not have an adequate explanation for this, other than the hypothesis that the L-DOPA producing plasmid imposes a lower metabolic cost on the bacteria
• Goal 5: Starting OD does not significantly change the data. So we should utilize the lowest OD that can be reliably measured. We decided to use a starting OD of 0.10 from now on.

Additional learning: If your experiment reaches max capacity too quickly it will muddle the results. In our case mCherry reached max capacity way before the end of the experiment, and before we could reliably extract biological parameters. This could explain why the logistic model, which is punished less by AIC due to having less parameters, performs better for mCherry increase, while Gompertz performed much better in the bacterial growth dataset.

What We Improved in the Model Before Doing Experimental Cycle #2

Setting Limits of the Y-axis

We found out that the limits of the y-axis should be set consistently by ourselves. So we introduced variables that stored the limits. This was retroactively added to the figures shown above.

OD to Bacteria Conversion

So far we approximated the OD to bacteria conversion as 800 million per million per OD value of 1.0. But the lab team by this time had done an experiment to determine the colony-forming units (CFUs) per OD per mL. This we turned into a formula that converts measured OD to the number of CFUs:

The last factor of 2×10⁵ was due to the dilution that allowed us to accurately count the CFUs. The protocol we used to determine CFU can be found in Wet lab - experiments.

Experimental Cycle #2

For this experimental cycle we wanted to test and validate our most important level 2 plasmids. We now also had a correct plasmid with the arabinose kill switch (Wet Lab - engineering cloning). So we were very curious if these constructs showed reduced growth even when not induced because of leaky expression, or if the kill switch could be activated to halt growth. In addition, in the first cycle we only looked at E. coli. But now we wanted to include the commensal bacteria: Pseudomonas alcaligenes.

Goals:

• Goal 1: Investigate if the kill switch stops growth
• Goal 2: Investigate if leaky expression of the kill switch slows growth
• Goal 3: Investigate if Pseudomonas alcaligenes can also be accurately modeled with the Gompertz model
• Goal 4: Investigate if protein production continues even when bacterial growth flatlines, using fluorescent reporters

So the experiment was set up in such a way that for 2 wells per construct had no arabinose added during the experiment, 2 wells had arabinose added at the start of the experiment, and 2 wells had arabinose added after 4 hours.

Unfortunately, we could not see a complete halting of growth for the induced arabinose kill switch (Fig. 3). However, in Pseudomonas we found evidence for a small, but significant (3-4%) reduction in growth speed when arabinose was added to the kill switch. After we confirmed the sequence again, we looked in literature for explanations. We found that the arabinose promoter can be suppressed by glucose and other catabolites.³ Our medium contains yeast extract, and is therefore likely a repressor of the arabinose promoter. The lab team did additional experiments with minimal medium supplemented with glycerol to check. The data seemed to suggest that the kill switch only worked for the Pseudomonas, but it was not very conclusive (see Wet Lab - Engineered Function). Together these results indicate that Pseudomonas is somehow more vulnerable to leaky expression of the kill switch, but also that the kill switch does not work as intended.

Observation: Lower Growth Rates Across the Board Compared to Experimental Cycle #1

As can be seen in the representative images, for some reason, the bacteria grew at a lower pace in this experiment than the first. We checked if this was due to the formula we used to calculate the number of bacteria per mL from the OD. But that was not the case. The OD-values themselves also increased at a way slower pace. The most likely explanation is that the temperature from the CLARIOstar was different between the two runs, which impacted the growth rate.

What We Learned from Experimental Cycle #2

• Goal 1: The kill switch did not work. This is likely due to glucose and catabolite repression of the arabinose promoter.
• Goal 2: There is some evidence for the effect of leaky expression of the kill switch toxin on Pseudomonas growth. This is not conclusive, as the observed effect was very small (3-4% reduction of growth rate) and was not decidedly reproduced in other experiments by the lab team.
• Goal 3: The Gompertz model can accurately model growth of Pseudomonas.
• Goal 4: We were unable to see if protein production continues after growth flatlines. This was because the fluorescent measurements got corrupted after adding the arabinose for reasons we don't understand. This meant that for both GFP and mCherry fluorescence intensity we could not use the data in any meaningful way.

Additional learning: Small unintended changes might have drastic effects on the growth rate of bacteria. For some reason our bacteria grew way slower than in the first cycle. We expect that the temperature settings of the CLARIOstar were different, but can not prove that with the data currently available.

The code and all generated figures are available

The final code we used is available on the resources page to download and use yourself. We added a "Read me" file with instructions on how to use it. We added both the raw data and plate lay out files for both cycles so you can try the code for yourself if you would like!

We fit many datasets, and showing all of those graphs (>100) would make this page a pain to read. However, if you are interested in checking out the data you can find it on the resources page. We have divided the results from the first experimental cycle and the second experimental cycle.

Flux Based Analysis

---

Flux Balance Analysis (FBA) is a constraint-based modeling approach that allows us to predict metabolic flux distributions in our engineered bacteria. This is crucial for optimizing L-DOPA production and understanding the metabolic constraints that may limit therapeutic efficacy.

Why we chose FBA

Our initial goal was to model Pseudomonas alcaligenes; we wanted to find a way to predict L-DOPA output to adjust our practical work and to know what to expect from our in-vitro experiments. FBA modelling was chosen for said task.

FBA - flux based analysis, is a common technique that is used to look at how an organism in a steady state distributes its resources or fluxes across its metabolic network. Since such models are not available for our bacteria we took a metabolic model of Pseudomonas putida (iJN1463), which is a close relative with a similar metabolic network. Both species of Pseudomonas also show similar ability to process aromatic compounds, which is important since both catechol and L-dopa are aromatic. As our experimental approach included adding a gene for tyrosine-phenyl lyase (TPL) which is the enzyme that produces L-DOPA we added this enzyme and its related reactions to the model.

Setting the constraints and limitations of the model

We first had to set the major limitations of the FBA: it being state-state and its constraints. FBA is known as a constraint based model since it assumes optimal usage of resources. So we researched and set constraints that we found in literature for uptake of catechol, ammonia and pyruvate; precursors for the TPL L-DOPA reaction. Leonie van Steijn told us in an interview that the uptake rates are the most important constraints of the model, so we spend additional effort to make sure we could argue for reasonable rates.

Balance between biomass growth and L-DOPA production

We were wondering what the effect of our L-DOPA production could be on biomass growth. So first we let our FBA find the maximum L-DOPA production possible, which was determined to be 0.821 minimol per gram dry weight per hour (mmol/gDW/h). We then wanted to calculate the cost of producing a certain amount of L-DOPA. So we set the limits of L-DOPA production at 0 to the max, in intervals of 0.1 mmol/gDW/h, and then have the model find the maximum biomass growth rate for each value. From this we found that every 0.1mmol/gDW/h reduces the biomass growth rate by 0.0131 (Fig 1.).

Then we wanted to better understand what this meant on the relation between L-DOPA yield and biomass growth. Therefore we set L-DOPA production at 0-100% of the maximum in intervals of 10%, and made the model calculate what the effect was on growth for each level of L-DOPA production. Through this we found out that the effect of L-DOPA production on biomass growth follows a linear relationship, with 100% L-DOPA production reducing biomass growth rate by 13.3% (Fig. 2). The biomass growth is in mmol/gDW/h, but is calculated for 1gDW. This means it can be seen as the fraction of its own mass it gains per hour.

Now we understand the cost of L-DOPA production and optimizing yield on biomass growth rate. However, Leonie van Steijn told us that the model should optimize biomass first instead of L-DOPA production if we wanted it to be more biologically relevant. So we optimized the biomass growth rate, and then manually lowered it in steps of 10% to see what level of growth rate would allow L-DOPA production. Consistent with our earlier finding we found that L-DOPA production was maximed at a biomass growth rate at 86.7% (Fig. 3).

Discovering the limiting metabolites: catechol

Given the chemical pathway we expected that pyruvate and catechol could potentially be limiting factors in L-DOPA production. As pyruvate uptake is around 5 times that of catechol, we started by investigating the impact of catechol on L-DOPA production. However, given the earlier finding that biomass growth rate is partially sacrificed to synthesize L-DOPA we also took this into account. So we changed the catechol uptake rates from 0 to 1.2, which is about 0-1.5 times the normal uptake rate. Then we set biomass growth at different levels. For each catechol uptake and biomass growth rate we then extracted the L-DOPA production (Fig. 4).

Our analysis shows that at maximum biomass growth rate L-DOPA production starts only if the uptake rate of catechol is higher 0.8 mmol/gDW/h. However, we are surprised there is any L-DOPA production at all. This means that catechol is not the only limiting factor, but that an initial proportion of catechol is essential for maintaining optimal growth. Additionally, we found that at a growth rate of 50%, virtually all catechol is used for L-DOPA production, until values higher than 0.8 mmol/gDW/h are reached. At this point they diverge, indicating that above this value catechol again becomes necessary to maintain growth. This shows a delicate metabolic balance of catechol, where an initial portion of catechol is necessary for optimal growth, but also that an over-supply will eventually not directly translate to increased L-DOPA production.

Because of the interesting divergence at high catechol uptake rates we decided we wanted to study this further, so we made the model go from an uptake rate of 0 to 3 mmol/gDW/h. This gave us even further insight. Eventually all growth rates achieve the same increase in L-DOPA production per increase in catechol uptake, no matter what growth rate it is. But high growth rates initially increase less in terms of L-DOPA production than the low growth rate ones. This means that at low growth rates initially all catechol is available to be spent on L-DOPA, but as this requires other resources it eventually converts into L-DOPA only at a certain conversion rate of around 3 units of catechol per 2 units of L-DOPA. For high growth rates, there is initially no catechol available, but as the uptake increases it will eventually get converted into L-DOPA at the same rate.

Unraveling the effect of pyruvate as a carbon source on biomass growth and its interaction with catechol

Now that we knew that catechol is vital for a balance between optimal growth and L-DOPA production we started to look at the impact of pyruvate. Pyruvate can be used immediately as a primary energy source in the citric acid cycle (PubChem), and is therefore a critical component of cell growth. Therefore we started to wonder if changing the ratio of pyruvate to catechol would change the carbon yield. The carbon yield is what the amount of biomass growth you get per unit of substrate used.

Normally the ratio of pyruvate to catechol uptake is around 5, so we decided to model pyruvate uptake to catechol uptake rates between 0-10, to include both lower and higher than usual ratios. As can be seen in figure 6, increasing the pyruvate to catechol ratio decreases the efficiency with which carbon is incorporated into biomass. We had honestly expected the reverse effect, with increasing pyruvate uptake leading to increased growth. However, this result might indicate that increasing the pyruvate uptake means that other non-carbon components of the cell are undersupplied, reducing carbon-specific yield.

Initial modelling of the kill switch through ATP used for maintenance

In addition to that we attempted to simulate a killswitch in our FBA; but the design was very limiting. We tried to mimic the effect of MazF: fast degradation of mRNAs leading to increased metabolic cost of mRNA production and lower levels of proteins and enzymes. We approximated this by manually increasing the maintenance cost of the cell through increasing ATPM by 50%, and approximating the reduction in translation of enzymes by reducing the activity of all enzymatic reactions to 10%. From this we learned that biomass growth rate is reduced to 0.041gDW/h, but L-DOPA synthesis is possible even when the kill switch is activated (Fig. 6).

Finding that L-DOPA production happens even when the kill switch was activated was quite surprising, as we did not expect the bacteria to have resources left to produce L-DOPA. So we wanted to understand what part of this was due to the increase in maintenance cost compared to the reduction in translation rates. Therefore we decided to keep the 50% increase in ATP maintenance cost, but vary the ratio of translation (Fig. 7). Interestingly, we found that the reduction in translation rates did not impact biomass growth rate, and that we could actually get a maximum level of L-DOPA production that was equal to the level we could achieve without the increased maintenance costs. This suggests that L-DOPA synthesis does not require energetic input, but is also limited by available metabolites.

Our results made us wonder if the FBA model is limited in its capacity to model processes like transcription and translation effectively. As mentioned before, we modeled reduced translation by cutting enzyme reaction speed as transcription and translation are not processes that are part of the FBA model. So this attempt demonstrated to us some limitations of this simple FBA model. Not only is it in steady state and assuming optimal convergence, it also does not contain a full genome and important cellular processes like transcription and translation. So the accuracy of the model was not ideal.Therefore we decided that it would be better to increase mRNA degradation rates itself, as that more closely resembles the mechanism of action. We decided to attempt to incorporate gene-product reactions in a new model, which we will describe on the next page.

What we learned with the FBA

Our estimates for the L-DOPA production allowed us to attempt to estimate our yields in the wet lab. More importantly, it also allowed us to more intricately understand our model organism and the biosynthesis pathway of L-DOPA. We identified crucial limiting factors like catechol and pyruvate, and showed how the competition for these resources affect the growth of our bacteria.

We learned that L-DOPA production siphons away resources which hampers optimal growth. Supplying additional catechol can remedy this at a specific rate of around 3 units catechol to 2 units L-DOPA. The growth rate of the cell strongly influences how many resources can be shifted towards L-DOPA production and therefore dictates L-DOPA production at low and basal levels of catechol uptake.

Additionally, we found that the ratio of pyruvate to catechol impacts the biomass carbon yield. Although we had expected that increasing pyruvate to catechol ratios would lead to increased usage of pyruvate for cell growth, we found that it actually decreased the yield. This likely means that other resources in our model are relatively less available at the high pyruvate to catechol ratios, leading to reduced efficiency.

Lastly, we learned that biomass growth is almost zero at an increased ATP maintenance of 50%, and that L-DOPA production is not impacted by this. This means that L-DOPA production is not very energetically costly. However, we also learned that this version of the FBA is not suited to model the mazF toxin mechanistically. So we decided to incorporate gene product reactions and expand the capabilities of the model. This would also allow us to screen knock-out genes that would increase yield, and allow us to model the biological effect of our kill switch toxin on transcription and translation.

GPR-Extended Model

---

The final phase of our modeling work involved creating a complete genome-scale metabolic model with gene-protein-reaction (GPR) associations. This allows us to link genetic modifications directly to metabolic outcomes and perform in silico gene knockout studies.

Why we chose GPR/ME-model

Gene-protein-reaction / Metabolic and expression model was chosen due to its more realistic representation. As mentioned previously a basic FBA model is missing a lot; because in reality TPL reaction is not as simple, TPL enzyme consumes cellular resources, requires transcription and translation, consumes amino acids and ATP etc. We wanted to mimic those biological processes, so we added the GPR component, which incorporated a full genome, following proteins and reactions for our organism. And components of ME-model, so for each gene we had corresponding transcription, translation and corresponding ATP costs to facilitate all of those processes. This will give us even more insight into L-DOPA synthesis and metabolic balancing of the costs.

This would allow us to calculate and add energy and resource burdens caused by MazF and L-dopa production.

Model base

The Model was still adapted from Pseudomonas putida (iJN1463) genome scale metabolic model. However since that model lacked a full genome we utilized Pseudomonas alcaligenes ASM159728v1, available on genbank and ran BLAST sequences to map homologous genes between the two. Additionally, we utilized fasta sequence from the same source to account for protein folding achieving nearly complete mapping Pseudomonas alcaligenes.

We also incorporated:

• Standard TPL reaction: linked it to a TPL gene for modelling purposes.
• MazF pathway: modelled in a more biologically relevant manner. We modeled it at both transcriptional and translational levels and the MazF toxin was accurately represented to cleave mRNA and assigned an energy cost in its inactive state.
• Quorum sensing (QS): added to better demonstrate the build up and gene population dependent effects on production.

Model Environment

Uptake constraints from different environments were incorporated in the model. For our project we modeled the following two environments:

• In vitro: to model in vitro growth in our medium, for proof of concept.
• In vivo: to model a more gut-like environment, with limited nutrient and oxygen resources

Dynamic simulations

We attempted to perform

• single cell - metabolic expression in 1 bacterium under varying conditions and changes (ex. Catechol knockouts and MazF activation).
• multi cell - similar conditions to single cell but used to understand population effect; enhanced via addition of quorum sensing previously.

Technical challenges

Initially CobraME-CobraPy framework was utilized seeing as it was quite popular. Sadly in our experience it proved to be extremely unreliable and incompatible with modern dependencies/addons. Even when forced into a legacy version and all of the dependencies were downgraded it sadly refused to work.

We sadly only recently discovered that there is a much better modern alternative based on ETFL - docker setup; connecting to a jupyter notebook . That is also the model that will be available in the supporting resources.

All of the dependencies worked and the model ran all the way to the end. Sadly it seems there is no internal Ldopa flux as seen in Figure 1. Under further investigation we found that the model lacked sources for all amino acids, which of course disables translation and therefore enzyme production.

Conclusion

Unfortunately we lacked time to fully utilize, validate and analyse the extent that this model can bring. However, we did get a full run of the model that did track ATP usage. L-DOPA was not produced, which might be due to how our model optimized for goals. Under further investigation we found that for some reason there was no source of amino acids in the cytosol, which of course abolishes translation. This means that adding viable amino acid sources might resolve the issues with the model, but we did not have enough time to implement this.

We did however validate that ATP was assigned properly to transcription, which we consider a huge step towards a functioning GPR/ME model. Furthermore, we could confirm a metabolic cost for the introduced plasmid, showing that transformation indeed puts a burden on cells. Concerns about this impacting the genomic stability have been raised in literature, showing it is a truly important consideration for future synthetic biologists (Stirling et al., 2017).

We aimed high, and got quite far. We are proud of what we tried to accomplish, even if we have to admit some level of defeat.