This is the mathematical and computational modeling page! Here, you will learn about the engineering process behind each model we built for Helo. This includes what questions we wanted to answer, why we chose to use models we did, how our experimental work relates to modeling, and where you can go to download our software.
"Computational modeling is the use of computers to simulate and study complex systems using mathematics, physics and computer science" [1]. Models are composed of several variables, formulas, and schematics that help researchers explain what they're studying. There are three ways of applying modeling to a project:
The power of modeling comes from its efficiency. Hundreds of in silico experiments can be run in seconds simply by adjusting variables and equations. Then, the data collected from these simulations can inform our project in any of the three ways listed above.
We decided to utilize all three methods of applying modeling to our project. Testing a drug's efficacy in vitro is integral to proving its viability, especially for Exendin-4 which competes with the endogenous hormone GLP-1. Efficacy is tied to drug-receptor binding; "the pharmacological effect is produced by the drug binding to the receptor to either mimic or antagonize the receptor" [2]. The competitive binding model (CBM) predicts how well a given concentration of Ex-4 will bind to its receptor, taking into account how native GLP-1 interferes with binding. This model also compares dissociation constants for Ex-4 found throughout scientific literature we've researched, and how it affects binding potential.
Once Ex-4 is bound, several downstream events occur that eventually lead to insulin secretion, the hormone responsible for lowering blood sugar levels. Modeling these downstream effects details Ex-4's mechanism of action, pairing well with the binding kinetics described in the CBM. The dosage-effect model (DEM) clarifies this complex process and provides insight on what effect a given dose of our drug would have.
Another crucial aspect to our project is collaborating with others. One team we've helped is the Stony Brook University. Their project also involves producing a recombinant protein, so we developed a protein production model (PPM) to show the validity of S. cerevisiae as a producer of their protein. We were well-equipped for this task because we used S. cerevisiae as a host and had already written a protein production model for Ex-4.
The CBM follows traditional fractional occupancy conventions, as per the model outlined by Salahudeen and Nishtala [2]. We chose this formula for its accurate depiction of our system, where Ex-4 and GLP-1 behave like two drugs with the same effect, but different potency. The DEM adopts the Emax model, one of the most fundamental depictions of the dosage to effect relationship in pharmacodynamics. What made this model so attractive was its simplicity, and how accessible its derivation was. Both the CBM and the DEM were coded on Jupyter Notebook in Python 3. The PPM was built from the ground up using SimBiology, a MATLAB package for modeling biological systems. All parameter values were found in the literature across multiple databases, such as PubMed, NCBI, and various course textbooks.
A great way to further legitimize a mathematical model is to bring it outside of the literature and code by connecting it directly to results from the wet lab. Our models achieve this by influencing the experiments themselves, and taking inputs from experimental data. For example, the binding affinity assay we plan to do will demonstrate how well our product binds to the receptor of interest, GLP-1 R, via dissociation constant values. The values obtained from this experiment can be included in the CBM. This will allow us to evaluate our Ex-4's ability to compete compared to literature values from previous Ex-4 binding experiments. Another example comes from our partners in the Ku lab at the University of California, San Francisco, who will perform an in vitro pancreatic beta cell assay using our recombinant protein. Since the DEM estimates the relationship between Ex-4 and insulin production, the model gives a range of Ex-4 concentrations to test for the expected effect.
This model is based on fractional occupancy, the ratio of the number of receptors bound to Ex-4 to the total number of receptors, whether they are bound to Ex-4 or not [3]. The basic fractional occupancy equation (1), however, is not robust enough to account for the Ex-4 and GLP-1 relationship. As previously mentioned, both proteins compete for the same receptor on pancreatic beta cells, and produce the same downstream effect. This means that both proteins must contribute to the number of receptors bound.
Equation (2) is a more accurate representation of our system. Ex-4 and GLP-1 increase the number of receptors bound in relation to their dissociation constants, Kd, as well as the total number of receptors (see Table 1 for parameter descriptions). Ex-4 has been shown to bind to GLP-1 R with greater affinity than GLP-1 [4], which translates to a lower Kd for Ex-4. Dividing each ligand concentration, [L], by its respective Kd, reveals how this difference in binding affinity changes fractional occupancy.
We applied this formula over a range of Ex-4 concentrations, keeping GLP-1 concentration and GLP-1 Kd as constants. We also graphed multiple Ex-4 Kd values from the literature, so we could easily examine our experimental Kd in comparison. These literature values were initially collected using their own unique methods, therefore classifying fractional occupancy by Kd has the added benefit of exhibiting how these different experiments changed their reported values.
The DEM falls under the classification of pharmacodynamic modeling. "Pharmacodynamics is the study of a drug's molecular, biochemical, and physiologic effects or actions" [5]. The specific equation we selected is the Emax model, which assumes that all biological mechanisms have a maximum output [2]. Here, that correlates to a maximum effect, Emax, of insulin secretion from a pancreatic beta cell (see Table 2 for parameter descriptions).
Unlike the CBM, equation (3) does not involve native GLP-1, making this a more appropriate model because of its parallel to our current wet lab capabilities. The in vitro pancreatic beta cell test we intend to connect to the DEM consists of adding varying concentrations of Ex-4 to a culture of human beta cells, then measuring the amount of insulin secreted. This system does not contain any GLP-1, so we felt comfortable forgoing GLP-1 entirely. Future research could perform experiments on mammalian hosts, where basal activity caused by GLP-1 will be present alongside Ex-4's effect. Then, a more holistic model can be developed to account for that.
Ease of use is a key characteristic when choosing a model. Searching for each parameter's literature value can be extremely difficult and time consuming. Due to the relatively short duration of our project, maximizing the work not done was imperative. We found values for equation (3) in most of the papers we read describing Ex-4's effect, compared to other equations' parameters which required more work to ascertain. Therefore, picking the Emax model saved time by expediting our background research process.
Unlike the prior two models, the PPM finds its basis in MATLAB's biology-focused app, SimBiology. SimBiology is a program used to draw visual schematics of our system, and automatically write ordinary differential equations (ODEs) describing each step of said schematic. Biological mechanisms contain a virtually innumerable amount of factors, but establishing reasonable assumptions can help manage these factors.
Equations (4) through (10) demonstrate how every element associated with protein production changes with concentration over time (see Table 3 for parameter descriptions). Elements such as transcription factors, although known to affect protein synthesis, remain omitted to help streamline the model. Much like the DEM, searching for the literature values of these more obscure factors can also be too time consuming and inconsistent. It's important to note that these parameters represent initial values, i.e. the very beginning of protein synthesis. Over time, they will fluctuate based on the equations listed above. All formulas obey the law of mass action, which states a reaction's velocity stems from the concentrations of its reactants [6]. Take equation (7), for example, describing the change in mRNA concentration over time. Rate constants, k, dictate how significantly each of the constituents affects this change. Positive terms such as (k3 · DP) and (k5 · MR) promote change because their reactions generate more mRNA. The opposite effect occurs for the two negative terms, -(k4 · M · R) and -(k6 · M), because their reactions consume mRNA.
Figure 1. The competitive binding model, showing binding affinity of Ex-4 in the presence of GLP-1. Each line represents a different Kd value for Ex-4, acquired from the literature. Concentrations of Ex-4 that produce half maximal occupancy are plotted.
Figure 2. The dosage-effect model, showing insulin secretion from human beta cells after exposure to a range of Ex-4 concentrations. The EC50 value is plotted: the point where additional Ex-4 causes a slowly diminishing effect.
Figure 3. The protein production model, showing the steps contributing to protein S expression in a S. cerevisiae cell. Blue ovals represent biological reactants, and yellow circles represent biochemical reactions. This model is designed for Stony Brook University. See Table 3 for variable descriptions.
Results from the CBM have a sigmoidal shape (Figure 1). This can be attributed to the assumptions made in the fractional occupancy approach (1), where a maximum number of receptors exists. Once all of the receptors are saturated, excess Ex-4 can no longer bind, causing the curve to plateau. The location of each curve on the axes is exactly as expected. Higher Kd values graph themselves further up the x-axis, meaning more Ex-4 is required to reach half maximal occupancy. This corresponds to horizontal translation between the graphs. Half maximal occupancy is the point where half of all available receptors are bound by a ligand. These results give us confidence in our model to predict binding kinetics, allowing for the implementation of experimental Kd values. In other words, the CBM is another step to confirming our recombinant protein's viability, supported by the shape of the graphs matching the binding relationship we're studying. Now, we can attempt a binding affinity assay to reveal the dissociation constant of our Ex-4, and instantly determine its effectiveness by adding it into the model.
The curve resulting from the DEM has a sigmoidal shape (Figure 2), much like the CBM. In this system, fractional occupancy directly correlates with biological effect, so both figures resembling each other is reasonable. Sigmoid graphs are common among pharmacodynamic models due to the mechanism therapeutics undergo within the body. Early addition of Ex-4 induces little change, followed by a small range of significant effect, which then plateaus regardless of more Ex-4 supplement. This graph enables us to pursue an in vitro mice pancreatic beta cell test by providing concentrations of Ex-4 to evaluate. The DEM shows what concentrations will produce a measurable effect, and what concentrations to avoid testing. Then, experimental values can be cross-referenced with the predicted values to analyze our protein's potency. The range centers around the EC50 value, calculated to produce 6.65 pg of insulin per islet per minute.
For the PPM, the design of the schematic itself dictates how its equations are written (Figure 3). All factors mentioned in Table 3 are pictured within one S. cerevisiae cell, along with their respective reactions. The ribosome binding reaction, M + R → MR, is non-reversible despite evidence showing ribosomal dissociation can occur [7]. This assumption is justified because of the sheer difficulty in searching for literature values describing ribosome kinetics. To circumvent this, the PPM relies upon ribosome abundance rather than kinetics, clarifying the binding interaction into one step. Designing the PPM for Stony Brook University not only promotes collaboration, but also gives our team valuable experience working with tools like MATLAB and SimBiology. In the future, these skills can be applied to model protein production of Ex-4 in various conditions, giving us insight into future implementations of Helo. To read more about implementations, check out the Implementation page!
While we were searching for a suitable model for our system, we realized that computational modeling is a rigorous field of study, requiring an ample background in biology and programming. Our modeling team’s biggest challenge was finding an appropriate model with the resources we had in this niche subject. Specifically, how both Ex-4 and GLP-1 cause the same downstream effect, but still compete against each other. We decided to make our models adaptive to assist future iGEM teams studying similarly unique binding interactions, letting them model their protein of interest against endogenous ligands inside their system.
Our models are adaptive because they are dynamic and adjustable. That is, we avoided hard coding our parameters into the models. For example, in the CBM, our graph is able to calculate and plot fractional occupancy for the Kd values that users input. Our graph will also automatically find where fractional occupancy is exactly half.
All documentation can be found on our GitLab repository. Download information and instructions for usage can be found in the documentation. These notebooks can also be executed via Google Colab, so anyone with a web browser may access them.
Competitive binding of two-ligand and one binding site is expressed as
(1) Another way to describe this is as foccupancy because the efficacy of a drug/ligand in question is measured by the amount of receptor-ligand complexes that are produced. This means the fractional occupancy can be represented by the efficacy of the drug in comparison to the total maximum efficacy. For the rest of this derivation, we will be using foccupancy.
(2) The total receptor concentration is dependent on the total free receptors and the total receptors bound to either the labeled ligand or the unlabelled inhibitor or competitor.
Meanwhile, we can also associate from equilibrium that
(3) Rewrite the equation based on what we found in (2).
(4) By manipulating this statement
We can replace the original equation’s numerator, such that
(5) From this equation, we can multiply by (1/[R]).
(6) Cancel out the denominator from the numerator by multiplying each half by Kd1.
(7) Factor out the Kd1.
(8) If we want to, we can consider more specificity-the two ligands have different potency.
(9) Once again we can cancel out the [R] to simplify the equation, which gives us a final equation of
Ligand, L, binds reversibly to its receptor, R, which undergoes a process resulting in the effect, E. Kon and Koff are the rate constants for ligand-receptor binding and dissociation respectively [8].
(1) Kd can be described using the same variables and constants.
(2) The rate of binding can be described using a differential equation depending on change of concentration of ligand-receptor complex.
(3) Concentration of free receptor, [R], can be rewritten as the total concentration of receptor, [Rtot], minus the concentration of bound ligand, [LR].
(4) Assume the system is in equilibrium, i.e. [LR] won't change.
(5) Koff / Kon can be substituted with Kd, and Kon can be canceled out.
(6) Assume that the effect, E, has a direct relationship to [LR], and there exists a maximal effect, Emax, when the total amount of receptors are bound.
(7) Substitute the equation from step 6 in, taking the place of [LR].
(8) Because only one ligand can bind to GLP-1 R, we can assume that receptor occupancy has a direct relationship to effect. Therefore, Kd = EC50.
Each reaction abides by the law of mass action, relating the concentrations of the reactants and products via their respective rate constants, k.
D + P ↔ DP translates to k1 * DP - k2 * DP
DP → M + D + P translates to k3 * DP
M + R → MR translates to k4 * MR
MR → S + R + M translates to k5 * MR
M → null translates to k6 * M
S → null translates to k7 * S
| Parameters | Values | Notes | Sources |
|---|---|---|---|
| L1 | 0.01 to 10,000 nM | Conc. of Ex-4 in the system | Experimental range |
| L2 | 1 nM | Conc. of GLP-1 in the system | Experimental value |
| Kd1 | 6 nM, 1.38 nM, 0.19 nM | Dissociation constants of Ex-4 to GLP-1 R | Yap and Misuan [9], Gao and Jusko [10], Gao and Jusko [11] |
| Kd2 | 25.1 nM | Dissociation constant of GLP-1 to GLP-1 R | Zhao et al. [12] |
| Parameters | Values | Notes | Sources |
|---|---|---|---|
| D | 1 molecule | DNA abundance in the system | Experimental value |
| P | 30,000 molecules | RNA polymerase abundance in the system | Borggrefe et al. [14] |
| DP | 0 molecules | DNA · RNA polymerase complex abundance in the system | Relies upon previous reactions: DNA · RNA polymerase binding |
| M | 0 molecules | mRNA abundance in the system | Relies upon previous reactions: transcription |
| R | 187,000 molecules | Ribosome abundance in the system | von der Haar [15] |
| MR | 0 molecules | mRNA · Ribosome complex abundance in the system | Relies upon previous reactions: mRNA · Ribosome binding |
| S | 0 molecules | Protein S abundance in the system | Relies upon previous reactions: translation |
| k1 | 0.0216 per second | RNA polymerase binding rate | Darzacq et al. [16] |
| k2 | 0.1450 per second | RNA polymerase dissociation rate | Darzacq et al. [16] |
| k3 | 0.0197 per second | Transcription rate | Provided by SBU: 40 bp/second for 2031 bp |
| k4 | 1 per second | Ribosomal binding rate | Assumption |
| k5 | 3.024e-7 per second | Translation rate | Provided by SBU: 8 aa/second for 677 aa |
| k6 | 5.0217e-4 per second | mRNA degradation rate | Provided by SBU: avg. half-life of 23 minutes |
| k7 | 2.188e-5 per second | Protein S degradation rate | Provided by SBU: avg. half-life of 8.8 hours |