Team:Chalmers-Gothenburg/ModelContruction
Navigation
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).
Elements to model
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 the organisms and cells present 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.
These three models are extensively explained and discussed in following sections
Integrated Modeling
In our project, our modeling outcomes have contributed to our project design in several ways. Firstly, the model serves as a proof of concept to determine whether the project design is feasible or not. If none of the simulations give promising results the latter is likely to be the case. Also, the modeling allows the comparison of the efficacy of the two possible anti-cancer proteins in order to choose the most suited of both and their required effective concentrations. Finally, the modeling will give quantitative indication of the required S. boulardii dosage and treatment time.
How the modeling results influenced our project is further described in the integrated modeling section.
Modeling Overview
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, a dynamic block consisting of a system of differential equations based on exchange metabolite mass balances, and a kinetic block with uptake kinetic expressions for the exchange reactions of the metabolic models. The dFBA iteratively computes and adjusts the boundaries of each exchange of 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 models
Effect of anti-cancer agents
In order to determine the efficiency of p28 and myrosinase on cancer cell proliferation, an ODE model of the cell cycle was implemented. 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) (figure 4) based on the models created by Bagci, Vodovotz, Billiar, Ermentrout, & Bahar (2000) and Aguda (1999), respectively (FIGURES?). Parameter values and initial conditions for these systems were taken from Hamada et al. (2009). An illustration of the cell cycle arrest system is shown in the figure below.
The model by Hamada et al. (2009) was originally created to observe the effect of DNA damage on the cell cycle. For our purposes, the DNA damage part of the system was removed (maybe I should remove it both from code and from figure?). An illustration of the apoptosis system is shown below.
In the cell cycle model, the cell can either go into cell cycle arrest or undergo apoptosis. If the p53 levels are high enough, the concentrations of the species involved in apoptosis are elevated and the cell dies. To simulate this, a threshold was set on the concentration of caspase-3 (abbreviated casp3 in the figure), a species involved in the later stages of the apoptosis system. It was decided that if the caspase-3 concentration becomes higher than the initial concentration of caspase-3, the cell undergoes apoptosis. In conjunction to this, an extra term was added to the system with the motivation that the caspase 3 levels should be able to increase from zero. The term was constructed so that the p53 concentration has a direct effect on the concentration of caspase-8 (casp8), a species that indirectly affects the caspase 3 levels to elevate.
The p53 threshold, p53thresh was set to 0.40 μM (Hamada et al., 2009).
The purpose of the cell cycle model was to find the effectiveness of the anti-cancer agents. To measure this, we decided to simulate the percentage of cancer cells that undergoes apoptosis as a function of environmental anti-cancer agent levels. Stochasticity was introduced into the cell cycle model by two separate means. Firstly, the initial extracellular anti-cancer agent levels are drawn from a normal distribution around the actual environmental concentration. This is done in order to simulate fluctuations in the environment. Secondly, the uptake rate of anti-cancer agent into the cell is drawn from a normal distribution where the mean was found by optimization (see below). For both of the normal distributions, the standard deviation was set to half of the mean.
p53 dynamics
For the p53 dynamics part of the cell cycle model, a new system of ODEs was created. The system describes the degradation of p53 both by Mdm2 and Cop1, as well as by other means not related to Mdm2 and Cop1. An illustration is shown in figure 6.
The differential equations describing the illustrated p53 dynamics follow below.
The parameter values for the p53 dynamics were found with optimization (see below) except for the association, kCop11 and dissociation, kCop12, constants for Cop1. These values were set to kCop11 = 0.1 μM-1s-1 and kCop12 = 1.1 * 10-3 s-1 (Moscetti, Bizzarri, & Cannistraro, 2018).
The initial concentrations of the p53 system were found by the assumption that the system is in steady state prior to addition of anti-cancer agent. The initial concentrations depend on the parameter values as follows.
After the optimization, a stability analysis was performed to make sure that this set of initial concentrations gives a stable steady state.
Models of p28 and myrosinase
To evaluate and compare the effects of p28 and myrosinase on cancer cell proliferation, two different models were made from the cell cycle model. 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.
Model with p28
For the purpose of simulating the effect of p28 on the p53 levels and thereby cell proliferation, the following reactions were added to the cell cycle model.
The association, kp281, and dissociation, kp282, constants for p28 was set to kp281 = 0.86 μM-1 s-1 and kp282 = 5.6 * 10-3 s-1 (Moscetti, Cannistraro, & Bizzarri, 2017). The effect of p28 on the p53 levels are governed by two processes; the uptake of p28 and the binding of p28 to p53. It is assumed that the only reaction affected by p28 is the degradation of p53 via Cop1. This means that it is assumed that p53 can act even when p28 is bound.
The model with p28 was optimized using proliferation data for sulforaphane generated by Yamada et al. (2009) to find values for the parameters involved in the p53 system and the additional parameter added to the dynamics of caspase-8. The proliferation data concern breast cancer cells but for these purposes p28 is assumed to have the same effect on colon cancer cells. The optimization was performed with a genetic algorithm and a set of model parameters were evaluated as follows. The model was run a 100 times using the given parameter set to find the proliferation probability for several concentrations of p28. Preferably, the proliferation would have been estimated from a larger number of runs but in order to keep the running time within a reasonable range, a 100 runs was deemed sufficient. The ODE system was solved for 86,400 time steps, corresponding to the 24 hours at which the proliferation data was recorded. The objective function was set to maximize the inverse of the squared difference between the generated and the experimental data. The parameter values obtained are as follows.
A stability analysis was performed in order to make sure that this set of initial values renders a stable steady state. The stability analysis showed that this steady state is indeed stable. A pdf of the Mathematica code used for the analysis is attached below.
Observing the parameter values and the initial concentrations obtained from the optimization it could be suspected that the initial concentration of Cop1 and the degradation of p53 via Cop1 is overestimated while those for Mdm2 is underestimated. The reason for this is that Mdm2 is supposedly the main regulator of p53 activity (Wang et al, 2011). Another interpretation could be that the regulation of p53 via Cop1 is more pronounced in a cancer cell. However, it is reasonable to assume that the overestimation of the Cop1 activity is a direct result of the fact that the objective function used during optimization favors parameters sets where p28 has a noticeable effect on the cell proliferation. The overestimation of Cop1 activity suggests that the model is faulty or that there are several solutions to the optimization problem. To investigate this further, an identifiability analysis should be performed. Based on the outcome of the analysis it could be decided if the model should be changed or if parameters should be set prior to optimization. Constraints could also be inferred on the Cop1 and Mdm2 activities, to make sure that the Mdm2 in fact becomes the main regulator of p53.
Model with myrosinase
In order to test the effect of sulforaphane on cell proliferation, myrosinase and sulforaphane were added to the cell cycle model as follows.
It was assumed that the myrosinase concentration is constant during the time considered and that Michaelis-Menten kinetics apply (Hochkoeppler & Palmieri, 1992). The kcat value of myrosinase was set to kcat = 1.19 * 102 s-1 and the Km to Kmmyr = 160 nM * s-1 (Hochkoeppler & Palmieri, 1992). The concentration of glucosinolate was approximated based on an average content of 62.3 mg glucosinolate in 100 g of broccoli (Lewis & Fenwick, 1987). Assuming that the molecular weights of the glucosinolates are similar, the molecular weight of sinigrin was used to find the corresponding molar content. The molecular weight of sinigrin is 397.454 g/mol (Pubchem, n.d.), which gave a glucosinolate content of approximately 0.157 mmol per 100 g broccoli. Based on that a person eats a 100 g of broccoli in a day and that 10 liters of liquid passe through the gastrointestinal tract every 24 hours, the glucosinolate content could be approximated to 0.157/(10*24) mM ≈ 650 nM. Note that this is a rough approximation.
To see an effect of sulforaphane on cell proliferation, sulforaphane was assumed to linearly upregulate concentrations of Bax and p21 and to increase the levels of cytosomal cytochrome C (Tortorella, Royce, Licciardi, & Karagiannis, 2015). Bax, p21 and cytochrome C are all involved in apoptosis. The following reactions were added to the respective species in the model.
The parameter values for the uptake and action of sulforaphane was find with a genetic algorithm using data from Frydoonfar, McGrath, & Spigelman (2004) and the same objective function as for the p28 model. However, since the proliferation data is based on a 72 hour exposure time, the ODE system was solved for 259,200 time steps. The following parameter values were obtained.
The mean uptake rate of extracellular sulforaphane was found to be 2.711 * 10-8 s-1.
Outcome of models
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 259,200 time steps, corresponding to 72 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 since it will be limited by the concentration of glucosinolate in the colon. Moreover, the effect of p28 is likely to be underestimated for an exposure time of 72 hours. The reason for this is that the parameters for the p28 model were found using data obtained with an exposure time of 24 hours. To avoid this, the models could have been optimized to data obtained with the same exposure time. However, it would have been even better if the models had been fit to proliferation data obtained at several exposure times. In that way, the output of the models would become more sensitive to exposure time. Disregarding that, it should also be considered that the results from the alpha pheromone model suggest that the cells will produce more p28 than myrosinase at a given cell concentration.
Future work
In the implementation of the cell cycle model, parameter values and initial concentrations were obtained from Hamada et al. (2009). However, these values are based on the dynamics of a normal, healthy cell. For future work, the parameters in the model could be modified in order for the dynamics to better resemble those of a cancer cell.
In the process of optimization, an identifiability analysis should optimally be performed to determine the validity of the outcome. If the work on the cell cycle models would continue, this would be the next step. Based on the identifiability analysis it can be decided if more data is needed in order to fit the model to reality and if there are parameter values that cannot be found with optimization.
To develop the cell cycle models further, the models could be made more time responsive in the sense that cell proliferation is not only a matter of anti-cancer agent concentration but also something that is affected by how long time the cell is exposed to the anti-cancer agent. If the models had been fit to proliferation data available at several time points, they could have been rendered more realistic based on exposure time.
The optimal last step of this modelling project would be to implement the cell cycle models into the genome scale model to get a large-scale model that can be used to evaluate the effects of p28 and myrosinase on cancer cell proliferation.
Kinetic model of the alpha pheromone system
In order to simulate the alpha pheromone system and the production of anti-cancer agent, a model by Kofahl & Klipp (2004) was implemented. The model describes the alpha pheromone pathway from receptor activation to the production of Ste12. The system described by the model is illustrated in figure 8.
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 alpha 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 shown in the equations below.
As can be seen by the equations, it is assumed that the transcription and translation are linear to the Ste12 and the mRNA concentrations, respectively. The purpose of the constant Q is presented below.
The production of alpha pheromone via mRNA was added to the model in the same manner as for the anti-cancer agent production. In this way, a feedback loop was created where the concentration of alpha pheromone is regulated by itself.
Parameter values for the added reactions were found with optimization. The model was fit to the expression of GFP under the FUS1 promoter with the addition of 2000 nM alpha phermone. The data was obtained from iGEM Delft (2012). Prior to optimization, the GFP expression at the first time point was set as the baseline and the data was normalized. In order to fit the data, Bar1 was temporarily added to the model while Far1 was removed from it. In addition, pseudo reactions for mRNA and protein production replaced the separate productions of alpha pheromone and anti-cancer agent in order to fit the total protein production to the data. The initial concentration of alpha pheromone in the model was set to 2000 nM based on the concentration of alpha pheromone used in the experiments. For the optimization, particle swarm optimization was applied, using the squared difference between the data and the model output as the objective function. The following parameter values were obtained.
The figure below shows the model output, total protein concentration, together with the data used for optimization.
It should be noted that the data used for optimization is not quantitative and that the model output as a result is only qualitative. In addition to this, an identifiability analysis should optimally have been performed in order to find out if the solution to the optimization problem is unique. In this case, it is highly suspected that the parameter values that were optimized for are linearly dependent, or correlated, which means that some of them should have been set prior to optimization in order to get a unique solution. Furthermore, it can be questioned wether it is reasonable that the transcription is so much faster than the translation as the parameter values indicate.
In order to find separate parameters for the production of anti-cancer agent and alpha pheromone, translation, transcription and degradation were assumed to be linear with respect to protein length. Q in the equations above was set to $Q = \frac{\nu}{\nu + 1}$ where $\nu = \frac{length(mf\alpha)}{length(anti-cancer agent)}. The following lengths were used: 120 aa for $mf\alpha_2$, (Saccharomyces Genome Database, no date), 28 aa for p28 (Warso et al., 2013) and 527 aa for myrosinase (source?).
The relative p28 production obtained by simulation with 20 nM, 200 nM and 2000 nM alpha pheromone is shown in figure 10.
The production of myrosinase is shown in figure 11.
Observing the figures of p28 and myrosinase production, it can be seen that as the initial alpha pheromone concentration is increased, the maximum concentration of anti-cancer agent produced increases. It can also be noted that the anti-cancer agent concentration becomes zero after some time, illustrating the need for continuous stimulation with alpha pheromone in order to get continuous production of anti-cancer agent.
Comparing figures ?? and ??, it can be seen that, as a direct consequence of the assumption that transcription, translation and degradation are linear to protein length, the production of p28 is higher than that of myrosinase. It can also be observed that the peak of the p28 production is reached within 500 minutes while the production of myrosinase is more prolonged.
To investigate how changes in the alpha pheromone concentration affect the anti-cancer agent production, the model was run for initial concentrations of alpha pheromone between 0 and 1000 nM. Figure 12 shows the maximum concentrations of p28 and myrosinase obtained as a function of initial alpha pheromone concentration.
Figure 12 suggests that the anti-cancer agent production peaks at an initial concentration of approximately 20 nM of alpha pheromone. This would mean that the highest concentration of anti-cancer agent is obtained when the concentration of alpha 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 alpha pheromone, it will diffuse in the colon making the effective alpha pheromone concentration lower...something... Is this reasonable? And what about the threshold that we would like to see in the beginning?? (Can we make a connection to the GEM?). Can we say something about the absolute values of anti-cancer agent concentrations?
Combined Genome scale models
MISSING TEXT
MISSING TEXT
MISSING TEXT
Saccharomyces boulardii Genome Scale Metabolic Model
The first step in our genome scale metabolic modelling approach was to obtain a Saccharomyces boulardii model. Unfortunately, there was no readily available model for this organism in the literature, so we opted to generate our own. Since it can take anywhere from months to years to develop a well curated model, we decided to generate a Saccharomyces boulardii model using the closely related and well curated Saccharomyces cerevisiae GEM. Having access to this GEM allowed us to use a homology based approach to obtain the desired Saccharomyces boulardii GEM with the help of the RAVEN toolbox (Agren, 2013). Please refer to the script reconstruction.m to see our MATLAB code for this GEM reconstruction. The script is divided into the following sections:
1. Install RAVEN & Set Up Working Directory
Download and install the raven toolbox. We used RAVEN 2.0.2. The exact version of RAVEN used for reconstruction, curation, and simulations is also available on our github page
2. Obtain Protein FASTA Files & Perform BLAST
Use following links to obtain protein FASTA files for our organisms of interest:
Run a bi-directional BLAST using the FASTA files for our two organisms to obtain a BLAST structure. This structure contains information about which genes, and therefore reactions, from the template model are present in the target model.3. Obtain High Quality Saccharomyces cerevisiae GEM
Download the yeast GEM here. Check that model is functional.
4. Get Saccharomyces boulardii Homology Model
Use RAVEN function getModelFromHomology, which takes a template model and BLAST structure as inputs. A model is obtained by essentially subsetting reactions from the template model which are also present in the target model, based on the similarity of their protein FASTA files.
5. Saccharomyces boulardii GEM Curation & Gap-Filling
Add non-gene annotated reactions from Saccharomyces cerevisiae GEM. This includes spontaneous reactions, exchange reactions, pseudo-reactions, etc. which were not captured in the BLAST structure. The model is still not able to produce biomass, therefore gap-filling needs to be done. The RAVEN toolbox offers the fillGaps function, which uses a template model as a reference to fill gaps so that biomass can be produced.
6. Add Proteins of Interest to Saccharomyces boulardii Model
After the draft model for Saccharomyces boulardii was reconstructed, we added protein synthesis and exchange reactions for our peptides of interest: MFalpha2, Myrosinase, and P28. The protein synthesis reactions were constructed based on the amino acid content of each peptide, and account for the ATP needed for the addition of each amino acid each (~4 ATP/amino acid (Amthor, 2000)). While the RAVEN toolbox offers functions to add reactions and metabolites to a model, it does not readily construct the lengthy stoichiometric equations needed for each protein. Instead of writing these expressions by hand we wrote a function to do this, called makeProteinEqn, which is under review to be added to the RAVEN toolbox. This function takes amino acid sequences and protein names as inputs, while outputting balanced protein forming equations such as:
'26 Ala-tRNA(Ala)[c] + 20 Arg-tRNA(Arg)[c] + 31 Asn-tRNA(Asn)[c] + 36 Asp-tRNA(Asp)[c] + 11 Cys-tRNA(Cys)[c] + 13 Gln-tRNA(Gln)[c] + 24 Glu-tRNA(Glu)[c] + 39 Gly-tRNA(Gly)[c] + 15 His-tRNA(His)[c] + 28 Ile-tRNA(Ile)[c] + 44 Leu-tRNA(Leu)[c] + 35 Lys-tRNA(Lys)[c] + 10 Met-tRNA(Met)[c] + 30 Phe-tRNA(Phe)[c] + 28 Pro-tRNA(Pro)[c] + 35 Ser-tRNA(Ser)[c] + 31 Thr-tRNA(Thr)[c] + 11 Trp-tRNA(Trp)[c] + 36 Tyr-tRNA(Tyr)[c] + 24 Val-tRNA(Val)[c] + 2107 ATP[c] => 26 tRNA(Ala)[c] + 20 tRNA(Arg)[c] + 31 tRNA(Asn)[c] + 36 tRNA(Asp)[c] + 11 tRNA(Cys)[c] + 13 tRNA(Gln)[c] + 24 tRNA(Glu)[c] + 39 tRNA(Gly)[c] + 15 tRNA(His)[c] + 28 tRNA(Ile)[c] + 44 tRNA(Leu)[c] + 35 tRNA(Lys)[c] + 10 tRNA(Met)[c] + 30 tRNA(Phe)[c] + 28 tRNA(Pro)[c] + 35 tRNA(Ser)[c] + 31 tRNA(Thr)[c] + 11 tRNA(Trp)[c] + 36 tRNA(Tyr)[c] + 24 tRNA(Val)[c] + 2107 ADP[c] + 2107 phosphate[c] + Myrosinase[c]'
Refer to makeProteinEqn.m to see the MATLAB code for this function.
7. Simulations
Check that proteins can be synthesized.
8. Final Clean-up of Saccharomyces boulardii MODEL
Add information about model reconstruction, export model to excel, and save model as a .mat file.
Representative Gut Microbiome Genome Scale Metabolic Models
To model the gut microbiome we decided to use three representative species in terms of composition and biosynthetic capabilities. A similar approach was taken in a gut microbiome study by Shoaie et al. (2013). In this paper, the authors generated genome scale metabolic models for three key gut microbiome member species: Bacteroides thetaiotamicron, Eubacterium rectale and Methanobrevibacter smithii, which represent the main phyla Bacteroidetes, Firmicutes, and Euryarchaeota, respectively. We used these models, which were provided in the supplementary information section. In the subsequent sections, we will discuss the model curation performed on each model. Please refer to sections 1.2-1.4 of the simulationScript.m file to view our MATLAB code.
Bacteriodetes thetaiotaomicron GEM Curation
The objective function was set to growth. Exchange reactions for glutamine, histidine, lysine, phenylalanine, valine, threonine, isoleucine, leucine, methionine, and tryptophan were added to allow for uptake. Additionally, some metabolite names were modified for consistency with other models. For example, the metabolite “CO2” was changed to “carbon dioxide”. This facilitates the identification of relevant fluxes for the automatic formulation of metabolite mass balances in the dFBA simulations.
Eubacterium rectale GEM Curation
The objective function was set to growth. Uptake exchange reactions for the abovementioned amino acids, as well as for succinate, were added. Some metabolite names were modified for consistency, as with the previous model.
Methanobrevibacter smithii GEM Curation
The objective function was set to growth. Curation for this model was more extensive as it appears to be the least well-curated model. Running the gapReport revealed that the model contains many gaps and isolated subnetworks. By using the checkRxn function, it was discovered that a few reactants in the biomass pseudo reaction were not able to be synthesized by the model. To address this, we reformulated the biomass equation to exclude these unreachable metabolites, which allowed the model to grow. Next amino acid and formate uptake exchange reactions were added. Finally, some metabolite names were modified, as discussed previously.
Colorectal Cancer Genome Scale Metabolic Model
As we were unable to find a specific genome scale metabolic model for colorectal cancer in the literature, we consulted with a few researchers in our department. Lucky for us, Avlant Nilsson was working on a Hep G2 model and allowed us to use it. While the model is meant to be specific to liver cancer cells, it was still useful for us since different human cancers similar metabolic hallmarks (Hanahan, & Weinberg, 2011). This liver cancer model was derived from HMR 2.0, a well curated genome scale metabolic model for generic human cells (Mardinoglu, 2014). Curation for our cancer model was very similar to that of the previous ones, uptake exchange reactions were added for some amino acids, as well as the short chain fatty acids propanoate, butyrate, and acetate. The objective function was set to growth. Finally, some metabolite names were modified for consistency with other models. Please refer to section 1.5 of simulationScript.m to see the MATLAB code.
Human Gut Genome Scale Metabolic Model
We obtained a GEM for the human gut from the Human Metabolic Atlas tissue-specific model repository (Uhlen et al., 2015). This model was also derived from HMR 2.0. Since it was not appropriate to set the objective function of this model to growth, we used a modified human maintenance reaction as the objective function for this model. We took the human_proteinPool reaction and metabolite (a stoichiometric sum of amino acids needed for protein turnover) from the HMR 2.0 model, and added it as a necessary reactant for the maintenance reaction: ATP[s] + 6 H2O[s] + proteinPool[s] => ADP[s] + Pi[s] + H+[s]. Note that some additional water was added on the reactant side to make the model uptake water, as this is part of the biological function of gut cells. In the subsequent simulations, this additional water uptake will not have a significant influence as water is not limiting. Uptake exchange reactions for some amino acids, propanoate, butyrate, and acetate were added. Finally, some metabolite names were modified for consistency with other models. Please refer to section 1.6 of simulationScript.m to see the MATLAB code.
Development and use of COM-dFBA: Community Dynamic Flux Balance Analysis
We decided to carry out some form of community dynamic flux balance analysis simulations early on in our modelling approach development. The advantage of dFBA over traditional FBA is that it allows us to simulate dynamic behavior. However it does require knowledge of certain kinetic parameters for the uptake of metabolites from literature. At first we considered compartmentalizing each model to generate one large network, but this would have been tedious and would presented difficulties in determining a community objective function. To avoid this problem, we opted to not merge the models, and instead have them connect through mass balances of some key metabolites. The scheme is presented in the figure below. As can be seen, there are three interacting parts: the FBA block, the dynamic block, and the kinetic block. Based on the desired initial conditions, an initial FBA is performed for each model. Next, the fluxes for the key connecting metabolites are identified and mass balances are generated. An ODE solver then integrates the mass balances to determine the concentrations of the key metabolites, at a given time point. These concentrations are then used to determine what the uptake rates should be for the next FBA simulation, and thus dynamic simulation is obtained.
While there are existing tools for running dFBA simulations, such as the COBRA toolbox, DFBAlab, and DyMMM, we decided to create our own dFBA tool, since the scheme is relatively straightforward to implement in MATLAB. Furthermore, the most recent and cutting edge tool, DFBAlab, does not seem to be easily accessible online. We requested access to the appropriate person, but we have yet to get a response from them.
COM-dFBA, short form COMmunity dynamic Flux Balance Analysis
is a simple framework for carrying out dFBA simulations for complex communities, such as the gut microbiome. It consists of six functions and one script, which can be found in our github repository under the folders COM-dFBA/Scripts. It’s usage will now be described.Initializing Models
Once the models have been curated we move on the initialization step. To view our code, please refer to section 2.0 of simulationScript.m. First, we create a cell structure containing the models to be used: modelSbo, modelBth, modelEre, modelMsi, modelCan, and modelCol representing Saccharomyces boulardii, Bacteroides thetaiotamicron,Eubacterium rectale, Methanobrevibacter smithii, colorectal cancer cells, and colon cells, respectively. Next, we pass this cell structure to the initModels function, which outputs a structure containing the loaded models, called superModel, as well as other relevant information. Refer to initModels.m to see exactly what the function is doing. Next, we pass the superModel structure to the function superModelFluxes, which neatly summarizes the exchange fluxes of each model in a structure called fluxes. Note that these fluxes are calculated with no uptake limitations. Refer to superModelFluxes.m to see what the function is doing. Next, we create a cell structure with the names of the exchange metabolites for which mass balances will be generated. These include biomass for each model, substrates, and products. This structure, as well as the superModel structure, is passed to the checkFluxes function, which in turn outputs the structures metEquations, metIdx, bounds, blocked. This function is mostly meant to serve as a debugging tool to check that the metabolites have been correctly matched in each model, which can be double checked in the metEquations strucure. Refer to checkFluxes.m to see in detail what this function is doing.
Initializing Environment and Kinetic Parameters
Once the models have been initialized we move on the initializing the environment and kinetic parameters. To view our code, please refer to section 3.0 of simulationScript.m. First, we load data from the excel file called paramTables.xlsx, which can be found in the COM-dFBA/Scripts folder. This file contains four sheets, with information regarding Vmax values, Km values, initial and flow concentrations, and other miscellaneous parameters, respectively. The literature search and sources for these parameters are discussed in a further section. The excel data is loaded into a structure called params. Next, we pass the superModel and params structures to the function updateFluxes, which outputs the structures superModelConstrained, fluxesConstrained, boundsConstrained, and lbubConstrained. The function serves as a debugging tool and sanity check. It runs one iteration of the dFBA simulation to check that the bounds and fluxes are being adjusted properly, based on the chosen kinetic parameters and initial concentrations. The most interesting output to look at is perhaps fluxesConstrained, which displays similar information as the abovementioned fluxes structure, the difference being that the fluxes are now cosntrained by initial conditions and kinetic parameters. Refer to updateFluxes.m to see what the function is doing. Once everything is in working order, we finally move onto simulations.
Simulations
Once the environment and models have been initialized we can begin running simulations. Refer to section 4.0 of the file simulationScript.m to see the MATLAB code. To begin a simulation, we simply use the following snippet of code: [t,xa] = ode15s(@(t,x)f(t,x,superModel,params),[t0 tf], initialConditions ,odeoptions). As can be seen, we rely on one of MATLAB’s ODE solver for integration. More specifically, we use the variable step continuous implicit solver ode15s. This solver was determined to be the most efficient through trial and error of different solvers, including ode23 and ode113. Next, we see that the ode15s function calls on another function called f. Additionally, we provide the ODE solver with a time span, initial conditions, and some options This function f takes in the superModel and params structures. Please refer to the file f.m to inspect the code in detail. This function contains three distinct sections: 1. KINETIC BLOCK, 2. FBA BLOCK, and 3. DYNAMIC BLOCK. The kinetic block uses Michaelis Menten kinetics to calculate the maximum possible uptake bound for each metabolite and organism, which depends on how much of a given metabolite is present in the environment. The following expression is calculated in a nested for loop, for each organism i and for each metabolite j.
Next, the function uses these above calculated values to update the appropriate uptake bounds in each corresponding GEM. Finally, we compare the MFalpha2 concentration in the environment with the MFalpha2 threshold for anti-cancer protein production, and update the protein production in the S. boulardii model accordingly. The FBA block then solves the models with updated bounds to obtain exchange fluxes for each organism. If the solution for the S. boulardii model becomes infeasible due to protein production constraints, we attempt to successively diminish the protein production load until the model is solvable or protein production is halted. Finally, in the dynamic block we take the FBA solutions and construct mass balances for each exchange metabolite of interest. This is done in a loop, and the mass balances are formulated as follows:
In other words, the change in the concentration of metabolite X is equal to the sum of production of X, minus the sum of consumption of X, minus the difference between outflow and inflow of X multiplied by the dilution rate. The only exception is the mass balance for cancer biomass, which accounts for the effect of anti cancer proteins instead of having a dilution term. In the following equation the term ACE represents the Anti Cancer Efficiency (1/hr), MW represents the molecular weight of the anti cancer protein (g/mmol), and X(Myr) represents the anti cancer protein (in this case Myrosinase) concentration (mmol/L):
Note that the resulting expressions for biomass have units of g/L*hr, whereas the expressions for metabolites have units of mmol/L*hr.
Determining flows and kinetic parameters
The gastro-intestinal tract can be split up in three main compartments: the stomach, responsible for the acidic and enzymatic break down of the ingested food, the small intestine, which further digest the food and absorbs a large fraction of the nutrients and water, and finally, the large intestine, which includes the colon and hosts microorganisms that ferment non-digestable nutrients (Crane, 1975).
In order to determine the differential equations for all biomasses and metabolites in the dynamic block of the community model (see figure 13), these three compartments have been represented as successive continuously stirred tank reactor-like (CSTR) elements. This implies that they are assumed to be ideally mixed and of constant volume. While these assumptions can seem overly simplifying, they are necessary to be able to model this system.
With this global layout being sketched, volumes and flow rates have to be determined. Later on initial and inflow concentrations of metabolites and biomass will also be discussed.
Volumes and flow rates
A person ingests around 2 liters of food and drinks a day (Wikibooks, n.d.), to digest this, the body produces around 7l of secretions: 2,5l in the stomach, 4,5l in the small intestine and 0,2l in the large intestine (Wikibooks, n.d.). The small intestine in turn absorbs water accounting for 90% of this volume (Azzouz & Sharma, 2018). This sums up to an inflow of:
(4,5 + 2,5 + 2) * 0,1/24 l/day = 0,0375 l/h from the small intestine
0,2 l/day = 0,0083 l/h of secretions
Thus:0,0375 + 0,0083 = 0,0458 l/h round up to: 0,05 l/h
The large intestine absorbs a large part of the remaining fluids but this will not be included in order to have only one outflow to the large intestine; the faeces.
While in reality the volume of the content of the gastrointestinal tract is not constant, we estimate the volume of large intestine based on the transit time of its content. This transit time is highly variable amongst individuals and directly depending on the individual’s diet. The transit time for the large intestine will be approximated to 30h (Chaddock et al., 2014). Assuming that the content passes through the large intestine on a linear time scale we get a volume of: 0,05 l/h x 30h = 1,5l
Initial and inflow concentrations
With the framework of the simulation being set up, both initial and inflow concentrations of metabolites as well as biomass have to be determined. These will allow the model to compute final steady-state concentrations that are of interest.
Initial Biomass Concentrations
Four different types of biomass are included in the model; gut microorganisms, including 3 representative species of the gut microflora, the human colon cells, the human cancer cells and Saccharomyces boulardii. All of these except the last one are assumed to have no inflow but only an initial concentration.
The total dry mass of all microorganisms in the colon is 0,100g (Karlsson et al., 2012). 98% of the microflora consist of Bacteria and the remaining 2% are Archea (Karlsson et al., 2012). Of the Bacteria, 29% belong to the Bacteroidetes phylum, 56% to the Firmicute phylum and the remaining 15% are other species that are not included in the model (Karlsson et al., 2012). To keep the same proportion of bacteria to archea, the percentages of Bacteriodetes and Firmicutes have been proportionally increased to 33% and 65% respectively. With a volume of 1,5l as calculated above, this results in 22,29 g/l of Bacteriodetes represented by B. thetaiotaomicron, 43,04 g/l of Firmicutes represented by E. rectale and 1,33 g/L of Euryarcheota species represented by M. smithii.
Due to a lack of literature data, we have assumed that the number of colonocites is in the same range as the enterocyte cell count in the small intestine, which is 1,67*1010 (Bianconi, 2013).
With a cell volume of 1,4*10-12 liter (Wiśniewski, 2012), a density of 1,05 (Pertoft & Laurent, 1977) and 40% of its mass being dry weight (Kretsinger, Uversky, & Permiakov, 2013) the total dry mass enterocytes, and thus the assumed dry mass of colonocites equals:
1.67*1010 * 1,4*10-12 * 1,05 * 0,40 = 0,98*10-2kg
With a volume of 1,5l this gives a concentration of: 6,55 gdw/L
Metabolite inflow concentrations
The metabolite inflow used in the model is based on daily intake according to a "standard British-type diet" (Stephen & Cummings, 1980). This diet contains 407 grams carbohydrates, including 22 grams of fibers, 85 grams of proteins and 108 grams of fats. Fats are not included in the simulation because they do not play a major role in the fermentation. Additionally, key metabolites such as ammonium and phosphate are included. Respectively 4,018 g/day (WHO, 2003), of which the majority is secreted by the body, and 1,5g/day (EFSA NDA Panel , 2015)
Carbohydrates (glucose)
Dietary carbohydrates are a wide variety of molecules ranging from simple monosaccharides to complex and resistant fibers. The carbohydrates of interest in this model are those that are fermented in the colon, namely fermentable fibers and resistant starch (Chebrut, Barry, Lairon, & Durand, 1995). They are consumed by the gut microorganisms to produce Short Chain Fatty Acids (SCFA). The other types of carbohydrates will either be majoritarily taken up in the small intestine or passing through unaffected. The composition of carbohydrates varies greatly depending on diet type. Here, it is assumed that 20% by mass of dietary carbohydrates are fermentable. To keep the number of exchange metabolites in the model low, all carbohydrates are modeled as glucose.
Amino Acids
In the large intestine, amino acids are both used as building blocks for protein production in all present cells and as substrate for fermentation. 12 grams of proteinaceous material enter the large intestine every day (Hughes, Magee, & Bingham, 2000). It is mainly composed of proteins and peptides and around half of it originates from the diet while the body secretes the rest. All proteinaceous material is modeled as free amino acids and only essential amino acids are included as metabolites. The proportions of amino acids entering the large intestine are assumed to be the same as in food as determined by (Lenders et al., 2009) and scaled down to 12g/day.
Oxygen
While the inside of the colon is anaerobic, blood vessels deliver oxygen to the human cells. This is simulated by adding oxygen in the inflow but only allowing uptake by the human cells. It is assumed that oxygen is in excess thus the concentration is set to a high value.
Secondary Metabolites
Other metabolites included in the simulation are carbon dioxide, acetate, propanoate, butyrate, succinate, formate, ethanol and methane. These are all secondary metabolites, meaning they are not present in the inflow.
Michaelis-Menten Kinetic parameters for substrate uptake
Uptake rates are based on Michaelis-Menten kinetics. Thus, they behave according to the equation:
Vs = (Vmax*S)/(Ks+S)
For each substrate and each organism both a maximal uptake rate Vmax and a Michaelis constant Ks are needed. These are only rarely well researched, thus in some cases rough assumptions had to be made. The values for these constant are summarized in the excel table accessible here as found in literature and converted into the right units for our purpose.
In order to make the simulations output realistic results, parameters have been changed to obtain final biomass and metabolite steady state concentrations that correspond to what would be expected from literature. The changes have been summarized in this excel file. For the third simulation; with the colonocites, gut microbiome and S. boulardii, the final concentrations of the second simulation have been used as initial conditions and an S. boulardii inflow of 1.0 gdw/l is added. The values are summarized in this excel file
Simulation Results
Simulation Design
To investigate the effect of our engineered organism, we set up three distinct in silico experiments:
- S. boulardii and Cancer
- Gut microbiome and Colon
- Gut microbiome, Colon, and S. boulardii
Experiments were run for 200 or 500 hours, depending on what was being investigated. Perhaps the most interesting setup is found in experiment number one, where we can see the interaction dynamics between S.boulardii and colorectal cancer. The rationale behind experiment number two is to obtain a sort of control to compare against experiment number three, where we introduce the engineered S.boulardii to the gut microbiome environment.
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 alpha pheromone. Just before the 50 hour mark, enough alpha 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 alpha pheromone concentration, specifically after the peak around hour 50, is due to the fact that less alpha pheromone can be synthesized once myrosinase production begins. This seems reasonable considering there are amino acid and protein pool constraints. The production levels of myrosinase relative to alpha pheromone concentration is validated by the myrosinase kinetic model.
Results Simulation 2: Gut Microbiome & Colonocytes
Missing text
Results Simulation 3: Gut Microbiome, Colonocytes & S. boulardii
Missing text
Integrated modelling
The purpose of the mathematical models created in this project was to make an impact on the course of the project and to be able to make inferences about the biological systems that we wanted to introduce into yeast. Even though the models we made all have areas for improvement, the modelling still made an impact on 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 and genome scale models 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 alpha pheromone system illustrates how the feedback loop of alpha pheromone works. It shows how the concentration of alpha 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 told 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 genome scale models, 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.
Alpha pheromone threshold
The alpha pheromone model gave additional insight into our anti-cancer agent production system. The results from the simulations indicate that the initiating threshold of alpha pheromone for the production of anti-cancer agent is lower than expected. Based on these findings, we should look into the system more closely and consider to use a weaker version of the FUS1 promoter. This is of importance since we do not want our S. boulardii to produce anti-cancer agent unless yeast cells has accumulated as a consequence of the presence of cancer cells.
Choice of anti-cancer agent
Based on the p28 and myrosinase models, we decided to move on with only myrosinase as the potential anti-cancer agent. The reason for this was that the models indicated that the efficiency of sulforaphane would be better than that of p28. The alpha pheromone model showed that cells would produce more p28 than myrosinase, which would perhaps favor choosing p28. However, myrosinase has the additional benefit that it can keep producing sulforaphane as long as there is any enzyme and substrate (glucosinolate) present, whereas the anti-cancer potential of p28 is limited by its own concentration. Had the results from the cell cycle models been added to the model of S. boulardii and cancer cells, it is expected that even further conclusions could have been drawn regarding anti-cancer agent efficiency.
S. boulardii content in pill
The simulation of
Future work
The models created in this project gave us valuable insights into project design and helped us improve our visualised product. While this is true, the modeling could still be subject to improvement. Below, we bring up model improvements that did not fit the time scope of the project.
For all of the kinetic models, optimization was used to find parameter values. However, in the process of optimization, an identifiability analysis should optimally be performed to determine the validity of the outcome. If the work on the alpha pheromone model and the cell cycle models would continue, this would be the next step. Based on the identifiability analysis it can be decided if more data is needed in order to fit the models to reality and if there are parameter values that cannot be found with optimization.
In the implementation of the cell cycle model, parameter values and initial concentrations were obtained from Hamada et al. (2009). However, these values are based on the dynamics of a normal, healthy cell. For future work, the parameters in the model could be modified in order for the dynamics to better resemble those of a cancer cell. To develop the cell cycle models further, the models could be made more time responsive in the sense that cell proliferation is not only a matter of anti-cancer agent concentration but also something that is affected by exposure time. If the models had been fit to proliferation data available at several time points, they could have been rendered more realistic based on exposure time.
When it comes to the genome scale models, the next step would be to add all of the GEMs implemented in this project, including the gut cell, the cancer cell, the gut microbiota andS. boulardii together in order to get an abstraction of the whole gut system. In addition to this, glucosinolate should be added to the gut inflow to get a better understanding of how the diet can affect the efficiency of myrosinase.
The optimal last step of this modeling project would be to add the results of the alpha pheromone and cell cycle models to the system of GEMs. In this way, we would obtain a large-scale model that can be used to evaluate the effects of our yeast on cancer cell proliferation and to further improve our product.
References
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