Modeling Overview

Modeling Objective

Before elaborating on the modeling part of this iGEM project, a good question to ask is why would we model our project? In general, creating a mathematical model for synthetic biology can give precious insight into a complex system, which ideally leads to predictive power over the outcome of this system and allows the improvement of the experimental design.

As described by Ledley (1960), modeling is a cyclical and dynamic process where the formulation of the model is repetitively changed to account for new literature data, experimental data and previous modeling outcomes. This new formulation yields a new outcome, which in turn can serve as a basis for a new modeling cycle or an improvement of the experimental design (see figure 1).

Modeled Elements

In this project, several aspects can be modeled. Because the yeast is modified to produce anti-cancer proteins as a response of an engineered pheromone sensing feedback loop, one model of interest is the pheromone production and signaling cascade and its corresponding anti-cancer protein production (see blue dotted line in figure 2). Another model consists of the anti-cancer protein effect on the cancer cells, specifically the effect on the cell cycle of the cancer cells (see green dotted line in figure 2). While these models are based purely on reaction kinetics, the third element to model is more complex: the growth of, and interaction between all the types of involved cells in the gut (see orange dotted line in figure 2). Thus including the probiotic S. boulardii and the cancer cells but also the key microbiome species and the healthy gut cells. All these types of cells interact by competing for metabolites and through the anti-cancer protein. The community interaction model thus includes metabolic models of all cells and the two first models. To read all about the construction of these models, please consult the modeling design page.

To give an overview of the modeling done in this project, the core consists a community interaction model supported by two kinetic models. The community interaction model is a dynamic Flux Balance Analysis (dFBA) framework, consisting of three blocks: an FBA block with genome scale metabolic models (GEMs), a dynamic block consisting of a system of differential equations based on exchange metabolite mass balances, and a kinetic block with kinetic uptake expressions for the exchange reactions of the metabolic models. The dFBA iteratively computes and adjusts the boundaries of each exchange metabolite of every metabolic model, based on Michaelis-Menten kinetics. The kinetic models additionally provide kinetic input to the dFBA for the production of pheromones and anti-cancer proteins and the efficacy of the anti-cancer protein. Both the kinetic and interaction models are formulated using values and parameters from literature. The modeling output enables us to improve our experimental design, which can yield different experimental results that can in turn be used to improve the model. This cyclic flow of information is illustrated as a flowchart in figure 3.

Kinetic Model - Alpha Pheromone Feedback Loop

In order to simulate the α pheromone system and the production of anti-cancer agent, a model by Kofahl & Klipp (2004) is implemented. The model describes the α pheromone pathway from receptor activation to the production of the transcription factor Ste12. The system described by the model is illustrated in figure 4.

In this project, Bar1 is removed from the genome of the yeast. Because of this, the concentrations and rate of changes of Bar1 and all of its deriviatives included in the model by Kofahl & Klipp (2004) was set to zero. To get the anti-cancer agent production as a function of the α pheromone concentration in the environment, additional reactions where added to the model. The additional reactions include production of anti-cancer agent mRNA and protein as well as production of α pheromone. For more details about the construction of this model, please consult the modeling design page.

To investigate how changes in the α pheromone concentration affect the anti-cancer agent production, the model was run for initial concentrations of α pheromone between 0 and 1000 nM. Figure 5 shows the maximum concentrations of p28 and myrosinase obtained as a function of initial α pheromone concentration.

Figure 5 suggests that the anti-cancer agent production peaks at an initial concentration of approximately 20 nM of α pheromone. This would mean that the highest concentration of anti-cancer agent is obtained when the concentration of α pheromone in the environment is around this value. A value of 20 nM can seem low considering that we do not want production of anti-cancer agent unless several yeast cells have accumulated. However, it should be taken into account that as the yeast cells produce α pheromone, it will disperse in the colon making the effective α pheromone concentration lower than the one produced.

Kinetic Model - Cell Cycle

In order to determine and compare the efficiency of p28 and myrosinase on cancer cell proliferation, one ODE model for each anti-cancer agent is built. The first model includes p28 and its effect on the p53 levels in the cell. The other model includes the enzymatic activity of myrosinase and the effect of sulphoraphane on cell survival. The cell cycle model is based on the work of Hamada et al. (2009) and consists of three parts; p53 dynamics, cell cycle arrest and apoptosis. The cell cycle arrest and apoptosis parts of the model were implemented according to Hamada et al. (2009) (simplified schematic shown in figure 6) based on the models created by Bagci, Vodovotz, Billiar, Ermentrout, & Bahar (2000) and Aguda (1999). The construction of these two models is described in details on the modeling design page.

To test the effect of p28 and myrosinase on cell proliferation, the models were run for anti-cancer agent concentrations between 0 and 200 μM. For simplification, the concentration of sulforaphane was set directly and the conversion of glucosinolate by myrosinase was left out from the simulations. For each concentration of anti-cancer agent, cell proliferation was calculated based on a 1000 realizations evaluated for 86,400 time steps, corresponding to 24 hours. The results are shown in the figure below.

Observing figure 7, it can be seen that the effect of sulforaphane seems to be stronger than that of p28. It should be noted, however, that the effect of myrosinase might be less strong than that of sulforaphane since it will be limited by the concentration of glucosinolate in the colon. Moreover, it should also be considered that the results from the α pheromone model suggest that the cells will produce more p28 than myrosinase at a given cell concentration.

Community dynamic Flux Balance Analysis

The first step in designing a community dFBA framework is to determine which cell types are to be included. Firstly, our community dFBA framework must contain the GEMs for colorectal cancer cells, human gut cells, and S. boulardii, since they are key aspects of our project. The question of how to model the complex system that is the gut microbiome is a more complicated problem. We decided to use GEMs of three representative species, in terms of composition and biosynthetic capabilities, which were taken from a gut microbiome study by Shoaie et al. (2013). In this paper, the authors generated GEMs for three key gut microbiome member species: Bacteroidetes thetaiotaomicron, Eubacterium rectale and Methanobrevibacter smithii, which represent the main phyla Bacteroidetes, Firmicutes, and Euryarchaeota, respectively. Below is a schematic depicting the system modeled by our community dynamic Flux Balance Analysis framework. Note that the schematic below corresponds to a scenario of a healthy person without colorectal cancer.

However, we are also interested in modeling the specific metabolic interactions between colorectal cancer and our engineered yeast. Below is a schematic showing this interaction.

Once we established the system to be modeled, we move onto the implementation of our dFBA scheme. We developed COM-dFBA, short for COMmunity dynamic Flux Balance Analysis, a simple framework for carrying out dFBA simulations for complex communities, such as the one shown above. COM-dFBA makes use of the RAVEN toolbox functionalities, and consists of six original functions and one script, which can be found in our GitHub repository under the folders COM-dFBA/Scripts. To understand the workings of COM-dFBA it is important to know that the chore of this simulation is an iterative cycle consisting of three parts: the kinetic block, the FBA block, and the dynamic block. Based on user-defined initial metabolite concentrations, the kinetic block first determines the upper uptake rate of each metabolite j by each organism i, based on Michaelis-Menten kinetics. These rates are then used to update the allowed uptake bound of each exchange metabolite j by each organism i. Next, an FBA simulation is performed for each model in the FBA block. In the dynamic block, the fluxes for each exchange metabolite j is identified in the FBA solution of each organism i, and mass balances are generated for each exchange metabolite. An ODE solver then integrates these mass balances to determine the concentrations of each exchange metabolite j at a given time point t. These new metabolite concentrations are then used by the kinetic block to recalculate the allowed uptake bounds of metabolite j by organism i. This continues in a cyclic fashion until a steady state is reached, and thus dynamic simulation is obtained.

To investigate the effect of our engineered organism, we set up two distinct in silico experiments:

1. S. boulardii and Cancer
2. Gut Community Interaction
   • Gut Microbiome and Colon before S. boulardii Inflow
   • Gut Microbiome and Colon with S. boulardii Inflow
   • Gut Microbiome and Colon after interrupting S. boulardii Inflow

The first simulation was run for 500 hours. In this in silico experiment we can see the interaction dynamics between S.boulardii and colorectal cancer.

In experiment number two, which is actually three consecutive simulations, we aim to investigate the effect of adding S. boulardii to the gut microbiome. In the first sub-simulation, we determine the steady state composition of the gut microbiome without S. boulardii. In the second sub-simulation, we then introduce an inflow of S. boulardii to the steady state obtained in the first sub-simulation. Finally, in the third sub-simulation, we stop the inflow of S. boulardii. The rationale behind this is to see the extent to which S. boulardii disrupts the gut microbiome composition, to see if S. boulardii will become out competed by the gut microbiome once the dosage/inflow is halted, and to see if the gut microbiome can return to its pre-S. boulardii composition once the dosage/inflow is halted.

Results Simulation 1: S. boulardii & Cancer Cells

In this simulation we aim to investigate the interaction between a growing colorectal cancer and our engineered yeast. As can be seen on the top subplot, both the colon cancer and S. boulardii biomass are near zero at the starting time point. The cancer begins to grow in its characteristic exponential fashion, while the yeast accumulates and begins producing α pheromone. Just before the 50 hour mark, enough α pheromone accumulates to trigger the production of the anti cancer protein Myrosinase. This causes the growth rate of the cancer to slow down, resulting in a global maximum of cancer biomass between 200 and 250 hours. After this phase there is enough myrosinase in the environment to “overtake” the growth rate of cancer, and we see cancer biomass begin to drop until it reaches zero just after 400 hours. S. boulardii approaches a steady state of around 12 g/L at the end of the simulation. The observed behavior of the α pheromone concentration, specifically after the peak around hour 50, is due to the fact that less α pheromone can be synthesized once myrosinase production begins. This seems reasonable considering there are amino acid and protein pool constraints. The production level of myrosinase relative to α pheromone concentration is validated by the myrosinase kinetic model. The system of ODEs defined by the exchange metabolite mass balances was solved for 1342 time steps, and the simulation took approximately 4.92 hours to complete. This simulation was run using MATLAB 2017b on the Hebbe computer cluster, which is part of the Chalmers Center for Computational Science and Engineering. The Hebbe cluster is built on Intel 2650v3 CPU's, and the system consists of a total 315 compute nodes (total of 6300 cores) with 26 TiB of RAM and 6 GPUs.

Results Simulation 2: Gut Community Interaction

As discussed previously, this simulation is actually composed of three successive simulations steps with the final concentrations of each step determining the initial concentrations of the successive one. The first simulation step serves to initialize the model, the initial biomass and metabolite concentrations found in literature and summarized on the Model Design page HERE PUT LINK are used to find the final steady-state reached. In the next step, S.boulardii is added in the inflow at a concentration of 1.00g/l, which corresponds to a continuous daily intake of in total 1.2 g of yeast. THis is done to observe how it will establish itself and how this influences the rest of the gut microbiome. Then, in the third and final step, the inflow of S.boulardii is put to 0g/l again to simulate the end of the treatment and to check whether S.boulardii will remain in the gut or be washed-out.

The simulation result can be seen in the figure below, where the biomass and metabolite concentrations are plotted over time. Vertical delimit the simulations parts.

Firstly, it is noteworthy that the final S. boulardii concentration is 1.05g/l while the inflow is 1.00g/l. While this is is not an enormous difference, it does mean that the S. boulardii is not only reaching the concentration in the inflow but is also able to grow.

S. boulardii induced change of gut microbiome concentration, this can be quantified as the final divided by initial concentration for each cell type, as presented in Table 2 below.

The change of Firmicute biomass between after introduction of S.boulardii is negligible. However, the Bacteroidetes and Methanogen growth is being influenced by the S. boulardii presence. It is though important to note that we know from literature that the exact composition of the gut microbiome is naturally variable (Karlsson et al., 2012). Especially the proportion of methanogen is variable, Karlsson et al. (2012) measured it to be 2.0%±4.3. These results are still reasonable within the boundaries from literature and thus we can conclude that the introduction of S. boulardii as simulated here does not disrupt the healthy gut microbiome.

Additionally, we can conclude that, because the S. boulardii concentration returns to values very close to 0g/l in around 50 hours after stopping the inflow of S. boulardii, it is unable to survive on its own, without continuous dietary intake. Also, the other gut microbiome species return to concentrations virtually equal to those before the treatment. This is a positive outcome because patients would want the engineered yeast to leave their gastro-intestinal tract shortly after the treatment, without having their gut microbiome altered.

Integrated modeling

The purpose of the mathematical models created in this project is to make an impact on the course of the project and to be able to make inferences about the biological systems that we introduce into yeast. Although the models all have areas for improvement, the modeling made an impact on the overall project and the work in the wet lab. In this section we explain in what ways it did so.

Proof of concept

First and foremost, the kinetic models and GEMs implemented in this project could be used as proof of concept and for illustrations of how our product will work in practice.

When it comes to the kinetic models, the kinetic model of the α pheromone system illustrates how the feedback loop of α pheromone works. It shows how the concentration of α pheromone in the environment induces the MAPK cascade that in turn leads to the production of anti-cancer agent. The cell cycle models with p28 and myrosinase tell us how the anti-cancer agents potentially affect the cell cycle. They also illustrate how p28 and sulforaphane can lead to the induction of apoptosis and thereby effectively kill cancer cells.

For the GEMs, the results from the simulations of the gut microbiota with and without S. boulardii indicate that the yeast has no dramatic effects on the composition of the gut microbiota. It demonstrates how S. boulardii can survive alongside gut microbial species while also not harming the patient. Moreover, simulations with S. boulardii , suggest that S. boulardii can produce myrosinase without depleting its amino acid resources. This means that the yeast can kill cancer cells while also growing.

Aguda, B. D. (1999). A quantitative analysis of the kinetics of the G2 DNA damage checkpoint system. Proceedings of the National Academy of Sciences, 96(20), 11352–11357. https://doi.org/10.1073/pnas.96.20.11352

Wiśniewski, J. R., Ostasiewicz, P., Duś, K., Zielińska, D. F., Gnad, F., & Mann, M. (2012). Extensive quantitative remodeling of the proteome between normal colon tissue and adenocarcinoma. Molecular Systems Biology, 8, 611. https://doi.org/10.1038/msb.2012.44

Crane, R. K. (1975). The Physiology of the Intestinal Absorption of Sugars This review has to do with the physiology of the intestinal absorption of sugars and should properly begin with a brief discussion of the components of the physiological system which carries out this i, 2–19. https://doi.org/10.1021/bk-1975-0015.ch001

EFSA NDA Panel. (2015). Scientific Opinion on Dietary Reference Values for phosphorus. EFSA Journal, 13(7), 4185. https://doi.org/10.2903/j.efsa.2015.4185

Wikibooks. (n.d.). Medical Physiology/Gastrointestinal Physiology/Secretions - Wikibooks, open books for an open world. Retrieved October 7, 2018, from https://en.wikibooks.org/wiki/Medical_Physiology/Gastrointestinal_Physiology/Secretions

Stephen, A. M., & Cummings, J. H. (1980). THE MICROBIAL CONTRIBUTION TO HUMAN FAECAL MASS. Journal of Medical Microbiology, 13(1), 45–56. https://doi.org/10.1099/00222615-13-1-45

Kofahl, B., & Klipp, E. (2004). Modelling the dynamics of the yeast pheromone pathway. Yeast, 21(10), 831–850. https://doi.org/10.1002/yea.1122

Stephen, A. M., & Cummings, J. H. (1980). THE MICROBIAL CONTRIBUTION TO HUMAN FAECAL MASS. Journal of Medical Microbiology, 13(1), 45–56. https://doi.org/10.1099/00222615-13-1-45

Wang, L., He, G., Zhang, P., Wang, X., Jiang, M., & Yu, L. (2011). Interplay between MDM2, MDMX, Pirh2 and COP1: the negative regulators of p53. Molecular Biology Reports, 38(1), 229–236. https://doi.org/10.1007/s11033-010-0099-x

Tortorella, S. M., Royce, S. G., Licciardi, P. V., & Karagiannis, T. C. (2015). Dietary Sulforaphane in Cancer Chemoprevention: The Role of Epigenetic Regulation and HDAC Inhibition. Antioxidants & Redox Signaling, 22(16), 1382–1424. https://doi.org/10.1089/ars.2014.6097

Lev Bar-Or, R., Maya, R., Segel, L. A., Alon, U., Levine, A. J., & Oren, M. (2000). Generation of oscillations by the p53-Mdm2 feedback loop: A theoretical and experimental study. Proceedings of the National Academy of Sciences, 97(21), 11250–11255. https://doi.org/10.1073/pnas.210171597

Kretsinger, R. H., Uversky, V. N., & Permi︠a︡kov, E. A. (Evgeniĭ A. (2013). Encyclopedia of metalloproteins.

Bianconi, E., Piovesan, A., Facchin, F., Beraudi, A., Casadei, R., Frabetti, F., … Canaider, S. (2013). An estimation of the number of cells in the human body. Annals of Human Biology, 40(6), 463–471. https://doi.org/10.3109/03014460.2013.807878

Bagci, E. Z., Vodovotz, Y., Billiar, T. R., Ermentrout, G. B., & Bahar, I. (2006). Bistability in Apoptosis: Roles of Bax, Bcl-2, and Mitochondrial Permeability Transition Pores. Biophysical Journal, 90(5), 1546–1559. https://doi.org/10.1529/biophysj.105.068122

Azzouz, L. L., & Sharma, S. (2018). Physiology, Large Intestine. StatPearls. StatPearls Publishing. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/29939634

Pertoft, H., & Laurent, T. C. (1977). Isopycnic Separation of Cells and Cell Organelles by Centrifugation in Modified Colloidal Silica Gradients. In Methods of Cell Separation (pp. 25–65). Boston, MA: Springer US. https://doi.org/10.1007/978-1-4684-0820-1_2

Bianconi, E., Piovesan, A., Facchin, F., Beraudi, A., Casadei, R., Frabetti, F., … Canaider, S. (2013). An estimation of the number of cells in the human body. Annals of Human Biology, 40(6), 463–471. https://doi.org/10.3109/03014460.2013.807878

Karlsson, F. H., Fåk, F., Nookaew, I., Tremaroli, V., Fagerberg, B., Petranovic, D., … Nielsen, J. (2012). Symptomatic atherosclerosis is associated with an altered gut metagenome. Nature Communications, 3(1), 1245. https://doi.org/10.1038/ncomms2266

Frydoonfar, H. R., McGrath, D. R., & Spigelman, A. D. (2004). Sulforaphane inhibits growth of a colon cancer cell line. Colorectal Disease, 6(1), 28–31. https://doi.org/10.1111/j.1463-1318.2004.00488.x

Pubchem. (n.d.). Compound Summary for CID 23682211 - Sinigrin. Retrieved October 7, 2018, from https://pubchem.ncbi.nlm.nih.gov/compound/Sinigrin#section=Top

Hughes, R., Magee, E. A., & Bingham, S. (2000). Protein degradation in the large intestine: relevance to colorectal cancer. Current Issues in Intestinal Microbiology, 1(2), 51–58. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/11709869

Moscetti, I., Bizzarri, A. R., & Cannistraro, S. (2018). Imaging and kinetics of the bimolecular complex formed by the tumor suppressor p53 with ubiquitin ligase COP1 as studied by atomic force microscopy and surface plasmon resonance. International Journal of Nanomedicine, Volume 13, 251–259. https://doi.org/10.2147/IJN.S152214

ledley, R. S. (1960). Report on the Use of Computers in Biology and Medicine - Robert Steven Ledley - Google Boeken. Retrieved from https://books.google.se/books?id=J5grAAAAYAAJ&printsec=frontcover&hl=nl&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false

WHO. (2003). Ammonia in Drinking-water - Background document for development of WHO Guidelines for Drinking-water Quality. Guidelines for drinking-water quality (Vol. 2). World Health Organization. Retrieved from http://www.who.int/water_sanitation_health/dwq/ammonia.pdf

Warso, M. A., Richards, J. M., Mehta, D., Christov, K., Schaeffer, C., Rae Bressler, L., … Das Gupta, T. K. (2013). A first-in-class, first-in-human, phase I trial of p28, a non-HDM2-mediated peptide inhibitor of p53 ubiquitination in patients with advanced solid tumours. British Journal of Cancer, 108(5), 1061–1070. https://doi.org/10.1038/bjc.2013.74

Moscetti, I., Cannistraro, S., & Bizzarri, A. R. (2017, November 20). Surface plasmon resonance sensing of biorecognition interactions within the tumor suppressor P53 network. Sensors (Switzerland). https://doi.org/10.3390/s17112680

Chaddock, G., Lam, C., Hoad, C. L., Costigan, C., Cox, E. F., Placidi, E., … Spiller, R. C. (2014). Novel MRI tests of orocecal transit time and whole gut transit time: studies in normal subjects. Neurogastroenterology and Motility : The Official Journal of the European Gastrointestinal Motility Society, 26(2), 205–214. https://doi.org/10.1111/nmo.12249

Hamada, H., Tashima, Y., Kisaka, Y., Iwamoto, K., Hanai, T., Eguchi, Y., & Okamoto, M. (2009). Sophisticated Framework between Cell Cycle Arrest and Apoptosis Induction Based on p53 Dynamics. PLoS ONE, 4(3), e4795. https://doi.org/10.1371/journal.pone.0004795

Lewis, J., & Fenwick, G. R. (1987). Glucosinolate content of brassica vegetables: Analysis of twenty-four cultivars of calabrese (green sprouting broccoli, Brassica oleracea L. var. botrytis subvar. cymosa Lam.). Food Chemistry, 25(4), 259–268. https://doi.org/10.1016/0308-8146(87)90012-4

Hochkoeppler, A., & Palmieri, S. (1992). Kinetic properties of myrosinase in hydrated reverse micelles. Biotechnology Progress, 8(2), 91–96. https://doi.org/10.1021/bp00014a001

Yamada, T., Mehta, R. R., Lekmine, F., Christov, K., King, M. L., Majumdar, D., … Das Gupta, T. K. (2009). A peptide fragment of azurin induces a p53-mediated cell cycle arrest in human breast cancer cells. Molecular Cancer Therapeutics, 8(10), 2947–2958. https://doi.org/10.1158/1535-7163.MCT-09-0444