Biological systems are complex, yet out of the apparent complexity comes order, and the circadian rhythm is an example of it. In a traditional reductionist point of view, all biological processes ultimately boil down to chemical processes happening inside the body. However, the rate of chemical reactions are generally temperature dependent, yet experiments find that even simple organisms like fruit flies exhibit regular circadian rhythm over a wide temperature range. Thus, the circadian rhythm is indeed a robust system that can adapt to external disturbances. In this project, motivated by the above question, I study a model in a paper on modeling the circadian rhythm of fruit flies. The paper and results therein suggest the possibility of structures involving not only negative feedback of the protein PER on the gene per, but also a positive feedback loop through dimerization of PER. By doing so, it is found that the system can exhibit circadian like behavior over a wide parameter range. In addition, with reference to another paper, I also study a proposal on how a robust circadian rhythm can be obtained through additional mechanism which regulates the rate “constants” of chemical reactions. Lastly, it is also pointed out that further connection with mechanism like inter-cellular communications may be important too.
The anatomical distribution of multiple proteins within complex tissue is predicted to aid in the determination of musculodystrophy. The data from multiplexed proteomic imaging method of array tomography is a very rich data set that enables detection and characterization of each pixel with multiple numbers of markers. The tissue architecture of the samples and thus the inter relationship of multiple cells can be preserved and studied at resolutions that are higher than that of other conventional microscopy methods by a factor of 5-10. This essentially provides us with a rich high-dimensional dataset. Such dataset can be studied in order to characterize Changes in abundance, localization and distribution of single proteins. But also the pairwise relations from such dataset can be studied to characterize changes in the pairwise relationship of proteins in relation to the same complexes. The high dimensionality of this dataset will enable study of higher order and changes in the complex relationships of proteins as part of multi-partner ensembles as well as classification of compartments in terms of their associated heterogeneous protein markers. Having many markers per pixel will help to statistically increase the significance of classification of pairwise relationship, complex associations and classified compartments against sensitivity and specificity of recognition/detection. The underlying hypothesis for my model is that the false positive and false negative errors in the detection will not affect the pairwise and higher-order associations of proteins when analyzing/classifying with large number of markers and large distributions. For this the algorithm for modeling and simulation is as following:
1. Generate rich datasets.
2. Simulate different levels of noise in the dataset.
3. Confirm that with different levels of false positives and false negatives, the protein associations can still be extracted from dataset as long as the number of classifier and dataset is increased.
4. Prediction of the model: the sensitivity and robustness of the statistical model in extracting pairwise and higher order protein associations.
In this project, I attempted to develop a model for the proposed model-based statistical analysis that encompasses the alternate asymmetry of binding of DHPR to RyR in skeletal muscle.
It is commonly thought that every transition in the cell cycle must be bistable since it is dangerous to turn back after you entered the next phase of the cell cycle. In the G1 phase it is thought that the bistability is dependent on the retinoblastoma protein (RB). RB is known to sequester a transcription factor that induces CYE expression. Once a cell is ready for the transition the G1 cell cycle dependent kinase (CYD/CDK4) phosphorylates RB which results in the release of the transcription factor that induces CYE expression. CYE in turn activates the G1/S cell cycle dependent kinase (CDK2) which also phosphorylates RB, creating a positive feedback loop. I used Boolean network analysis and ODE analysis in XPP to show that shows a bistable behavior with hysteresis. Low CYD results in low CYE, once CYD is turned on CYE is high but turning CYD off doesn’t lead to a significant decrease of CYE.
However from experiments we know that there is another pathway wherein CDH1 targets CYE for degradation. Phosphorylation of CDH1 by CYD/CDK4 and CYE/CDK2 leads to inactivation of CDH1.
If the transition of G1 to S-Phase is bistable, the CDH1 branch is supposedly also bistable otherwise there will be degradation of CYE once the CYD signal is turned off. I tested this hypothesis with ODEs in XPP and found that the CDH1 branch indeed has an bistability domain dependent on the levels of CYD.
Next I translated the whole system into a Boolean network. And analysis of this systems shows that the G1-to-S phase transition is indeed bistable and depends on both the RB branch and the FzR branch.
Next I transferred the whole ODE system to COPASY and tried to use time courses from the single branches made in XPP to estimate the parameters of the full system. This system showed with this parameter setting that the bistability is mainly dependent on the pRB branch. This suggests that if you ectopically boost CYE mRNA production, you should be able to artificially induce the G1-to-S phase transition.
I planned to do a parameter sensitivity analysis to see how robust the system is but unfortunately I didn’t have enough time for doing this.
G-Quadruplexes (GQ) are nucleic acid super-secondary structures forming at G rich regions with a (G2+N1-5 )2+ motif. GQs consist of stacked quartets that are in turned formed by Hoogsteen hydrogen bonded Guanines, and are stabilized by potassium (and sodium) ions. GQs were initially discovered in telomeric sequences but later studies showed their enrichment in promoter and other geneic regions (Du et al., 2007; Lam et al., 2013). Their high enrichment in these regions of the genome suggests that they could play a role in gene regulation. However, most studies on quadruplexes have been done in-vitro especially with respect to their biophysical aspects. Even though GQs have been suggested to influence gene regulation, this possibility has not been thoroughly explored.
I explore the potential effects of GQs on transcription using mathematical models. I have used both deterministic model and stochastic model to study the transcription dynamics in the presence or absence of GQ. In both the models it has been assumed that GQs are downstream of the RNA polymerase binding site and four cases have been explored:
- No quadruplex
- Quadruplex in template strand
- Quadruplex in non-template strand
- Quadruplex in both the strands
Most parameters used in the model are adapted from the values reported by biophysical studies, while others have been assumed from guesses.
The models prove that G-quadruplexes can indeed regulate transcription. Future studies would include studying the effect of different parameters on the system and to perform experiments to validate the model.
Modeling Cdc42 GTPase dynamics in the control of polarized growth
The conserved GTPase Cdc42 is a key regulator of cell polarity and morphogenesis that controls processes such as actin-cable polymerization and exocyst localization. In fission yeast, Cdc42 activation exhibits oscillatory dynamics at the two polar growth zones. Anticorrelated oscillations at the cell tips result from positive and delayed negative feedbacks and competition for Cdc42 activators. One of these activators, the guanine nucleotide exchange factor (GEF) Gef1, undergoes a phosphorylation–dephosphorylation cycle that regulates its localization to its site of action at the cell cortex. Cortical localization of Gef1 is hypothesized to be dependent on microtubule-based delivery of a phosphatase to the cell tips.
In this work, I used XPP to model the ordinary differential equations that describe the positive and negative feedbacks and the contribution of a constant level of active Gef1 to Cdc42 activation. Model parameters were tuned to generate oscillations that were dependent on achieving a threshold level of Gef1, which is consistent with experimental observations. At these parameters, phase plane analysis revealed the existence of stable, unstable,
and saddle point nodes and a limit cycle. Bifurcation analysis showed that the system exhibits saddle homoclinic orbit bifurcations that begin at a Hopf bifurcation point and end at the saddle node bifurcation.
To model the role of microtubules in Cdc42 oscillations, I included expressions generating a periodic source of phosphatase to activate Gef1 at the cell tips. The period during which microtubules contact and shrink away from the cell tips was approximated based on published data. The model showed that microtubule contact could maintain Gef1 levels at the cell tips sufficiently for oscillations to occur, and the oscillations were driven by the microtubule period. Changing the microtubule period and amplitude resulted in irregularities in the oscillations. These predictions can be experimentally studied by observing oscillations of active Cdc42 in mutants that have altered microtubule dynamics.
Autophagy is a highly orchestrated degradative process that sequesters and degrades damaged or excess organelles in response to internal and external stress including starvation, oxygen limitation, microbial invasion and hormonal stimulation. The understanding of biochemical pathway of autophagy is important in determining the fate of cells that respond to stress including the challenge of chemotherapeutic agents. This project attempts to build a mathematical model for the activation of autophagy by nutrient deprivation (low concentration of amino acids) with feedback loop when the amino acids generated by protein lysosomal degradation increases the amino acid pool. To achieve this goal, I use MATLAB and XPP learned from the CSHL CCB summer course to build this model and do analysis. The model successfully simulates the change of amino acid level as a response to different nutrient status by autophagy. The system displays oscillation of amino acids level in both nutrient starvation and nutrient enriched status, but does not shows oscillation if starting from the steady state. Total amino acids eventually reach the same level regardless of the initial nutrient status. The bifurcation analysis confirms the simulation result.
Glutamate excitotoxicity is a pathological process implicated in stroke, traumatic brain injury and numerous neurodegenerative diseases. In culture, neurons exposed to the same excitotoxic stimulus can either survive the insult or undergo rapid necrotic or delayed apoptotic cell death. In our lab we are investigating the regulation of this switch-like behaviour. Using fluorescence data obtained in single neurons, we observed that ATP and AMPK rapidly recovered to homeostasis following transient excitotoxic perturbation. In contrast, the recovery of glucose was significantly slower. Interestingly, the extent of this delay correlated with the duration of subsequent survival.
To investigate this behavior from a computational perspective, I previously developed a MATLAB ODE-based model based on prior knowledge of signalling pathways. The model correctly resembled ATP and AMPK kinetics, but failed to resemble the delayed glucose recovery.
In this project, I implemented my model in both Virtual Cell and COPASI, and used the inbuilt parameter estimation functions to fit the glucose response. I confirmed that the model in its current state cannot completely explain the observed glucose kinetics. I next incorporated additional model reactions in an attempt to explain the discrepancy in glucose behaviour. However, extensive parameter estimation for this updated model failed to simultaneously fit the glucose response curve for both the steady-state and perturbation situations. Further analysis is necessary to identify specific reactions that can delay the glucose recovery without simultaneously affecting ATP.
Rac1, a monomeric RhoGTPase plays a critical role in regulating endothelial junctional permeability; however, the causal relationship between Rac1 activity at adherens junctions (AJs) and stability of Vascular Endothelial (VE)-cadherin adhesion, the main protein of AJ complex, is not well understood. My research interest is to determine the underlying molecular and signaling mechanism involved in Rac1-mediated VE-cadherin adhesion. I utilized photo-activatable and photo-switchable proteins, which allowed me to control the spatial and temporal activity of Rac1 and determine the effect of its activity at the level of cell-cell junction on the dynamic of VE-cadherin adhesion. My experimental data suggests that activation of Rac1 at AJs promotes VE-cadherin accumulation by reducing the rate at which VE-cadherin molecules leave the junction (ie. dissociation). Using mathematical modeling, I asked whether the observed decrease in the rate of VE-cadherin dissociation is sufficient to explain the Rac1-induced VE-cadherin accumulation. I generated 11 ODE questions to explain VE-cadherin association and dissociation (including dimerization and complex formation) and by using parameter estimation (COPASI) I showed that VE-cadherin dissociation at the rate of 1.348±0.015 1/sec is sufficient to promote VE-cadherin accumulation as seen in my experimental data. Interestingly, this rate is much faster than what I measured experimentally (0.0018±0.002 1/sec). However, it is important to note that the experimental measurement is the total amount of signal; hence the dimerization and complex formation influence the rate. Another Interesting observation was that the model suggests that the dimerization of VE-cadherin happens at much higher affinity compared to the complex formation, which is both unexpected and interesting since the complex formation has shown to be critical for stability of the junctions. For future studies, I’d like to design and perform experiments that can explain this phenomenon.
In haploid yeast, mating pheromone triggers a fate decision: arrest of the cell cycle in G1 and initiation of mating events. To do that, the MAP kinase of the pheromone response pathway (PRP) Fus3 activates the cyclin dependent kinase inhibitor Far1, which binds to and inhibits all three G1 cyclins complexes, Cdc28/Clns. Cdc28/Cln2 activity is essential to pass the START G1 checkpoint, while Cdc28/Cln3 is required earlier in G1 to drive the expression of Cln2. This has two important consequences: a) inhibition of Cdc28/Clns causes cells to arrest in G1, and b) since Cdc28/Cln2 can block pheromone response, its inhibition by Far1 acts as a positive feedback loop.
However, when cells are exposed to low doses of pheromone a new phenotype appears: some cells start responding to pheromone but after some time, stop the response and divide. We named this behavior “switching cells”
We know this switching behavior is present in cln1cln2 strain but absent in Kar4 strain. Kar4 is a transcription factor induced by pheromone and involved in karyogamy. Kar4 induces, with a time delay, the transcription of Cln3 (cyclin of G1). Therefore, my hypothesis is that the pheromone triggers not only the activation of the pheromone pathway but also, with a delay, the increase of Cln3 (cyclin of G1).
Using Copasi and the idea of a incoherent feedforward loop I built a toy model in order to understand how this two opposite branches can explain dynamically the switching behavior.
The goal is to use the insight that I won with the toy model to build a more realistic model for the pheromone pathway and the activation of Cln3 by pheromone.
The ability of the bacteria to switch from one carbon substrate to another (i.e to reorganize its metabolism) is one of the complex and highly regulated mechanisms that bacteria have developed for responding to their environment changes. The CSR (Carbon Storage Regulator) is a post-transcriptional regulator involved in the regulation of such mechanisms, the protein CsrA (main component of the CSR system) can activate or inhibit the translation of target mRNA’s and has been shown to activate the glycolytic metabolism while inhibiting the gluconeogenetic pathway. We are interested in understanding the role of the CSR system during the switch from glycolysis to gluconeogenesis.
To achieve this goal, we need to integrate the dynamical data in different levels of transcriptomics, proteomics, and metabolomics into one consistent model which could explain the interplay of the CSR system and the transition dynamics of interest. The data at hand includes time courses of interacellular metabolite concentrations measured with mass spectrometry, extracellular concentrations of glucose and acetate measured with HPLC, and enzymatic assays. And finally time series data for mRNA (transcriptomics). So far we only have data for delta-csrA.
We understand that this modeling project is complex and needs refining and going back and forth between experiments and modeling simulations. To put this into a manageable modeling effort in the context of this course, two different approaches were taken. First, we started at the metabolic level, and only used the affinity constants (Km) from the literature and estimated the missing Km’s and all Vmax’s by fitting to the metabolomics data. We got a relatively good fit. The next step would be to use the enzymatic assays to update the Vmax’s since we assume they scale linearly with the amount of enzymes and check to see if the predictions agree with the concentration data for the knock out strains. We expect to see disagreements which will help us learn more about the dynamics and update the structure of the model.
As a second strategy we looked at each enzymatic reaction one at a time, and incorporated the three levels of gene to mRNA to protein/enzyme which regulates that reaction in a sub-model. This model is supposed to help us locate the point of csrA regulation and confirm and add to what is in the literature about where csrA regulates the glycolytic and gluconeogenic pathways.
The simulations and optimizations were done in MATLAB for the first part of the project. We used COPASI for the modeling of the second part in order to explore that software and its capabilities.
Once Manon obtains the complete data sets, we could proceed with integrating all the sub models with the metabolic model and learn about the complex interactions of the different moving pieces of this multi-level integrative model.