Model | Stony-Brook

What is Mathematical Biology?

Mathematical biology integrates theoretical frameworks, mathematical models, and conceptual representations to explore the fundamental principles driving biological systems' structure, function, and behavior. By employing tools like MATLAB, we can develop precise mathematical descriptions of complex biological processes. These models enable the simulation and prediction of biological behavior, offering a powerful means of conducting quantitative analysis. They play a crucial role in validating experimental data, testing theoretical predictions, and revealing insights that may not be directly observable in experimental settings. Through modeling, we can construct systems of equations that mimic cellular activities and interactions.

Deterministic and Stochastic Modeling

Developing an effective miRNA detection kit requires knowledge of how the gene system works, and how sensitive the system is to be applicable. The concentration of miRNA in the blood varies per type of miRNA and the method of extraction, and the magnitude of concentration in blood can get as low as an attomolar. With a low working concentration in a noisy biological environment, it’s necessary to understand how the system behaves with these conditions.

Figure 1. Concentrations of individual miRNAs in the blood (Biovendor, 2018)

Modeling biological gene circuits is a fundamental approach to understanding how complex biological systems function. Modeling techniques can provide unique insights into the behavior of a gene circuit, helping to elucidate how genes and proteins are regulated and expressed. To model our system, we created deterministic and stochastic models, each of which helps us answer different questions.

Our deterministic model uses ordinary differential equations (ODEs) to describe the average behavior of a genetic circuit. These models assume a predictable path, given initial conditions of reactant concentrations and reaction rates. A stochastic model, or probabilistic model, represents the dynamics of the system using a series of random, discrete events that occur over time. Stochastic models are best for when there are a small number of molecules in the system and random fluctuations greatly impact the behavior. We used Gillepsie’s algorithm to simulate individual events such as the binding of miRNA to the Lac operator or the translation of mRNA.

Why use both deterministic and stochastic models?

Deterministic models give us an average or "expected" behavior of the system by assuming large numbers of molecules and ignoring noise. This is useful for establishing a baseline or a reference point for how the system should behave under ideal or average conditions. It helps in understanding the underlying dynamics, identifying key parameters (like binding rates or degradation rates), and validating initial assumptions. Furthermore, deterministic models allow for quick exploration of parameter space without the need for extensive simulations, identifying which parameters are most critical to system behavior.

Finally, by comparing the results from deterministic and stochastic models, we can validate our understanding of the system. For instance, if the stochastic model shows significant deviations from the deterministic predictions, it highlights the role of noise and stochastic effects, confirming that these effects are indeed important in our system.

Deterministic Model

Our miRNA detection system involves the inhibition of the dual regulation system of LacI and L7AE in the presence of miRNA. With the GFP gene further downstream of the plasmid, the dual regulation system will inhibit the expression of GFP under normal conditions. However, in the presence of miRNA-326, Ago2 will bind to miRNA to form a RISC complex. When the RISC complex binds to the complementary sequence of miRNA after the start codon of the mRNA, it slices the dual regulation system, resulting in the expression of GFP.

In MATLAB, we looked to model the concentration of GFP to test our various genetic constructs.

RISC Equation

The constants for these equations were mostly found through literature review. However, some had to be calculated through other means.

We used the J23100 Anderson promoter whose transcription rate was documented as 0.042663 PoPS (Gaston Day School, 2020). We multiplied this number by the copy number of 15 for our Ago2 plasmid and the number of plasmids present in the system. We know the molecular weight of our plasmid as each base pair is roughly 0.65 kDa (Eurofins Genomics, n.d.), and we can find the total number of plasmids through a known concentration of plasmid inserted into our CFS. We can then convert this to M/s by dividing by Avogadro’s number and using 1.5 mL as the volume as that was the size of our CFS. We repeated this for the dual repression plasmid, which has a copy number of 10. This is the new transcription rate. Our “old” transcription rate before we changed to CFS was calculated with a E. coli volume of 6e-16 L.

We’re using the NEB cell free system, and the exact translation rate of amino acids is not documented. However, it is derived from cell extracts of Escherichia coli S30, so we used the typical translation rates of E. coli instead. We divided the 12.1 amino acids per second translation rate by the number of total amino acids in our proteins to determine the translation rate for each one.

The protease activity in the NEB CFS is also not defined. We are still in the process of characterizing the protein degradation constants, so we will be using an E.coli dilution rate constant as a placeholder for the degradation rate. This dilution rate is modeled after the half life of the E. coli.

Since we don't have accurate degradation rates for our cell-free system at the moment, we don't think it would make sense to change our transcription rates as we believe this would lead to a more inaccurate model.

Constant	Description	Value	Source
kTRA(old)	Transcription rate of Ago2 plasmid (M/s)	1.77e-9	(Gaston Day School, 2020)
kTRA(new)	Transcription rate of Ago2 plasmid (M/s)	1.81e-11	(Gaston Day School, 2020)
kTLA	Translation rate of Ago2 (1/s)	0.0194	(Harvard University BioNumbers, 2024)
kdeg	general protein dilution in E. coli 1/sec	0.000577	(McKernan, 2015)
miRNA	Upper bound magnitude of miRNA in the body, depends on the miRNA	1e-12	(Khashayar et al., 2022)
k_on_ago	RISC formation constant in 1/Ms	2e7	(Wee et al. 2012)
k_off_ago	RISC dissociation constant in 1/s	7.7e-4	(Wee et al. 2012)

Table 1: Constants for RISC ODE equations

The change in Ago mRNA, Ago protein, and RISC complex, are written using standard synthesis - degradation.

As mentioned before, the translation rate of LacI and L7ae were calculated using a basal translation rate in E. coli. The hill's coefficient is assumed to be 1 as there is only one binding site for miRNA in Ago2, resulting in no cooperativity.

Constant	Description	Value	Source
kTLLI	translation rate of LacI (1/m)	0.0472	(Harvard University BioNumbers, 2024)
kTLLA	translation rate of L7ae (1/s)	0.1417	(Harvard University BioNumbers, 2024)
q	Hill's coefficient of miRNA and Ago2	1	-

Table 2: Constant for LacI and L7AE ODE equations

The change in Ago mRNA, Ago protein, and RISC complex, are written using standard synthesis - degradation.

The translation rate of GFP was calculated using the 12.1 amino acids/seconds basal rate found in E. coli. The LacI-LacO dissociation constant was found to be 592 nM in our literature search, however we’re using an enhanced LacI with a binding strength of greater than 100 compared to normal LacI (Semsey et al.). Therefore, we divided the dissociation constant by 100 to get 592e-11 M.

Constant	Description	Value	Source
kTR(old)	Transcription rate of dual repression plasmid	1.182e-9	(Gaston Day School, 2020)
kTR(new)	Transcription rate of dual repression plasmid	1.51e-11	(Gaston Day School, 2020)
kTLGFP	Translation rate of GFP (1/s)	0.0714	(Harvard University BioNumbers, 2024)
n	Hill's coefficient for LacI - LacO	2	(Semsey et al., 2013)
m	Hill's coefficient for L7Ae-k-turn	1	(Lilley, 2014)
k_on_LE	bimolecular rate of association (1/M*s)	8.4 × 10*6	(Wang et al, 2012)
k_off_LE	dissociation of L7Ae from Kt-7 (1/s)	.002	(Wang et al. 2012)
k_dLE	L7AE-kturn dissociation rate in M	2.38e-10	k_off_LE/k_on_LE
k_dLI	LacI - LacO dissociation rate in M	592e-11	(Du et al. 2019)

Table 3: Constants for mRNA and GFP ODE equations

The change in LacI, L7AE, mRNA, and GFP concentrations, are written using Hill's repression terms. LacI and L7AE are inhibited by the RISC complex. mRNA is inhibited by the concentration of LacI. GFP is inhibited by the concentration of L7AE.

Hill functions are mathematical expressions used to model cooperative binding or regulation in biological systems, especially where multiple molecules interact to produce a non-linear response, such as in gene regulation, ligand-receptor binding, or enzyme activity. They help capture how the binding of one molecule influences the binding of others, leading to a sigmoidal (S-shaped) response.

Hill's Coefficient:
n=1: The system shows no cooperativity, and the response is hyperbolic (similar to Michaelis-Menten kinetics). This means that each binding event is independent.
n > 1: The system shows positive cooperativity. The binding of one molecule makes it easier for subsequent molecules to bind. This is often seen in cases where multiple molecules (e.g., proteins) bind to the same substrate or regulatory site, leading to a sigmoidal curve.
n < 1: The system shows negative cooperativity, where the binding of one molecule makes it harder for others to bind.

Figure 2: Concentration of GFP expression over time without any repression, schematic from Shu, et al.

It is important to note that the basal concentration of GFP falls at ~2.5M and reaches equilibrium at around 1e4 seconds.

Figure 3: Concentration of GFP expression over time with dual repression, schematic from Shu, et al.

Figure 4: Concentration of GFP expression over time with LacI repression, schematic from Shu, et al.

Figure 5: Concentration of GFP expression over time with L7ae repression, schematic from Shu, et al.

We made two separate models for LacI transcriptional repression and L7ae translational repression to prove the necessity of a dual regulation system. Although individual repression from each do make an impact on GFP concentration, complete inhibition shows to be most effective.

Figure 6: Concentration of GFP expression over time with miRNA and Ago2 repressing the dual repression system, schematic from Shu, et al.

In this model, we can note that the concentration of GFP actually returns to its basal level of 2.5e-4 M with RISC inhibiting the dual regulation system. It takes roughly 5x more time to equilibrate to this concentration.

Figure 7: GFP expression with varying concentrations of miRNA

Figure 8: Using fluorescence intensity to find the concentration of miRNA

From Khashayar et. al. 20222022 and Špringer et al. 2021, we know that the standard magnitude of miRNA in the blood ranges from attomolar to picomolar. Analyzing these graphs, it is clear the sensitivity of our system is in the correct range of magnitude for miRNA detection. Our system seems to be the most sensitive from 10^-16 to 10^-14. There is a notable dose-dependent relationship, with increasing miRNA concentrations leading to higher GFP expression. The sensitivity of the system to miRNA concentrations can be inferred from the slope of the curve. If the system is highly sensitive, small changes in miRNA levels should lead to noticeable changes in GFP output, which is crucial for detection in low-concentration scenarios like attomolar ranges. Since it’s difficult to quantify the exact magnitude of miRNA-326 in humans, we want to make sure our system is sensitive to a wide range of miRNA concentrations.

Stochastic Model

Systems within cells do not occur as uniformly as described by our ODE modeling. Small changes can occur in the concentration of species due to shifts in equilibrium caused by a multitude of factors– like temperature or the concentration of other species.

Biological systems are inherently noisy due to this discrete randomness, aka stochasticity.

Incorporating an element of randomness into models allows for a closer look into a system’s resilience to change and consistency.

The Gillespie algorithm, also known as the Stochastic Simulation Algorithm (SSA), is designed to simulate the time evolution of chemically reacting systems with inherent randomness.

In our miRNA detection system, we model the behavior of miRNA, Ago2, the RISC complex, LacI, L7AE, and GFP. Given the low concentrations of miRNA (in the attomolar to picomolar range) and the potential for significant stochastic fluctuations, applying the Gillespie algorithm allows us to capture the probabilistic nature of these interactions.

Our Program

Our Stochastic Model was designed using an SSA. Just as in our ODE which simulated the behavior of miRNA, Ago2, the RISC complex, LacI, L7AE, and GFP, the initial concentration of miRNA was on the order of 1E-12 as is common in the blood of B-Cell Lymphoma patients. The initial concentration was 9E-12. This, because the stochastic model was run on the molecular level with each species in the system starting at 10 molecules. miRNA, however, had an initial molecule count of 15,000, to represent an initial amount of miRNA that is crucial to our system. After the stochastic model was run, molecule counts were converted to molar concentration through the formula

For our SSA to work properly, an initial concentration of 0 was impossible. The number of molecules was multiplied by 1E13 to scale up the system to an appropriate molecular count.

The same set of constants were defined for the Stochastic Model, some of which were changed by +/- 10-100 /s to tune the model for a more probable distribution.

1) kTR = 1.1820; Transcription rate for mRNA
2) kTLGFP = 0.714; Translation rate of GFP
3) kTLLI = 0.00472; Translation rate of LacI
4) kTLLA = 0.001417; Translation rate of L7AE
5) kdeg = 0.0000577; Degradation rate (protein dilution in E. coli)
6) kTRA = 1.77e-4; Transcription rate of Ago2
7) kTLA = 0.194; Translation rate of Ago2
8) k_dRI = 0.526e-7; Dissociation rate of RISC

ALL CONSTANTS ARE DEFINED IN M/S EXCEPT FOR kTR WHICH IS M/S2

A second random variable between 0 and 1 was defined at each time step. If this random number was between the values of 0 and the value of propensity 1, event 1 would occur. If the value of this variable was between propensity 1 and propensity 2, event 2 would occur. This continued for all 8 events. The interaction of all of the species were approximated with +/- 1 to show repression, formation, degradation, and other processes that ‘create or destroy’ molecules in the system. This shown here:

The final graphs were obtained :

Figure 9: Stochastic representation of our gene construct

Considering the approximate nature of this Stochastic Model, the prediction of GFP concentration reaching a steady state at around 2.25E-5 for most runs, compared to the 2.50E-4 predicted by our deterministic modeling, GFP expression in a randomized system is still sufficient for the detection of miRNA-326 given an miRNA concentration on the order of 1E-12. It is important to note that this model is particularly sensitive to changes of the number used to scale the system in this equation:

It is also very sensitive to change of the variable cell.free.system.volume. It is also important to highlight that the rate constants kTR, kTLLGFP, kTLLI, kTLLA, kdeg, kTRA, kTLA, k_dRI were all subject to change of between 10-100 if it helped tune the model to the expected final concentrations. This was accepted as the SSA makes a gross generalization that at each time step, concentrations change by 1 M. Additional scaling past that provided through propensities was necessary.

Parameter Sensitivity Analysis

To ensure our miRNA detection system is robust and sensitive across a range of miRNA concentrations, we try to understand how changes in model parameters affect the system's output. In biological systems, parameters such as reaction rates, binding affinities, and degradation rates can vary significantly due to experimental conditions, individual biological variability, and measurement errors. Using deterministic or stochastic models alone may not capture the full extent of this variability.

Therefore, we used Monte Carlo simulations to observe how sensitive the system's behavior is to changes in specific parameters. We specifically modified the following:

Parameter	Mean Value	Gaussian Standard Deviation
kTR	1.1820e-9	2.0e-10
kTLGFP	0.0714	.01
KTLLI	0.0472	0.005
KTLLA	0.1417	0.02
kdeg	0.000577	0.00005
k_dLI	592e-11	60e-11
k_dLE	2.38e-10	2e-11
kTRA	1.77e-9	2e-10
kTLA	0.0194	0.002
n	2	.5
m	1	.5
q	1	.5
k_on_ago	2e7	2e6
k_off_ago	7.7e-5	8e-6

Table 4: Standard deviation numbers for monte carlo simulations

Variations for these were determined based on assumptions made in literature, with most standard deviations ~10%. The hill’s repression coefficient for LacI and DNA had a deviation of 0.5 in literature, so this was extrapolated to the hill’s coefficient for miRNA and Ago2.

Through Monte Carlo simulations, we can observe a strong deviation in the miRNA sensitivity range compared to that of the original. This is likely due to the nonlinear effects of our genetic system. To quantify how variability in our system and which parameters are responsible for this, we ran monte carlo simulations to calculate the sobol indices of our parameters.

In global sensitivity analysis, the goal is to understand how much of the uncertainty (variance) in the output Y of a model can be attributed to the uncertainty in its input parameters X1,X2,…,Xd. The total variance of the output V(Y) can be decomposed into contributions from individual inputs and their combinations:

The first-order Sobol index S_{i} measures the proportion of the total variance in the output that is caused by variations in a single input parameter X_i, while all other input parameters are fixed. Mathematically, it is defined as:

E(Y∣Xi) is the expected value of Y when the input Xi is fixed, and the remaining inputs are allowed to vary. V(E(Y∣Xi)) is the variance of this expected value.

Si=1 means that the variance in the output is entirely due to Xi, and no other parameters contribute.
Si=0 indicates that Xi does not affect the output, and the variance is entirely due to other parameters or their interactions.
Typically, Si lies between 0 and 1, showing the fraction of output variance due to the specific input.

Assumptions:
The first-order Sobol index assumes that the output variance attributed to each input is additive and separable. This means that interactions between inputs are not considered in the first-order effect. If interactions between variables exist, they are captured by higher-order Sobol indices (e.g., second-order, third-order indices), but these are ignored in the first-order index.

Simulating Ago2-miRNA Binding with GROMACS

Our iGEM project focuses on developing a bacteria-detector system that glows green when it detects high levels of miRNA-326, a microRNA associated with B-cell lymphoma. As part of the dry lab component, we are using GROMACS to optimize the binding of human argonaute protein (AGO2) and miRNA-326 with potential applications to enhance its binding affinity. This work supports our main project by providing insights into the molecular interactions critical for the detection system.

Modeling Approach

We utilized GROMACS, with the leap-frog integrator and Particle Mesh Ewald (PME) for long-range electrostatics, to model and analyze the structural dynamics of Ago2 and miRNA-326. By simulating the binding of these two structures under accurate physical conditions, we can determine the most optimal binding position of Ago2 and miRNA-326. Further techniques such as umbrella sampling could allow us to identify mutations that could improve the binding affinity, thereby enhancing the sensitivity of our bacteria-detector system.

During this process, we also noticed that there is an incomplete crystallized protein file on RCSB for miRNA-122 bound to human argonaute. As a result, there is currently no accurate PDB representation of the Ago-miRNA complex. We decided to also generate a complete version of Ago2 bound with miRNA-326 using molecular docking and MD simulations.

Progress and Achievements

Tutorial Completion: We completed the GROMACS tutorial for lysozyme, which provided essential knowledge for setting up and executing molecular dynamics simulations. (http://www.mdtutorials.com/gmx/lysozyme/06_equil.html).

System Stability: Successfully stabilized AGO2, miRNA-326, and the docked version, in an aqueous environment, ensuring realistic simulation conditions that mimic physiological settings.

Simulation Visualization: Obtained a visual representation of Ago2-miRNA trajectory in free floating space, where we can observe binding affinity and molecular motion.

Further simulations and analysis are still in progress

Molecular Docking

We used Autodock Vina to find energetically optimized binding configurations of the two structures. We used an exhaustiveness of 8 and generated 9 docking poses.

Limitations: Autodock Vina only allows for 32 active torsions, and miRNA-326 has 168 active torsions. We were therefore greatly restricted in accuracy and flexibility, but the way we made it work was by keeping the C3’--O3’ bonds, O3’--P bonds, and C1’--N bonds, the most active as they were the most crucial for intermolecular interactions and folding in miRNA.

Steps Performed in GROMACS

(1). Generating the Topology File:
We created a topology file for both AGO2 and miRNA-326 using GROMACS tools like pdb2gmx. The topology file includes detailed information about the molecular structure, atom types, bonds, angles, dihedrals, and force field parameters needed for the simulation. This file is crucial for setting up the system, defining how atoms will interact during the molecular dynamics simulation. Additionally, it ensures compatibility with the chosen force field, allowing accurate simulation of protein-nucleic acid interactions.

(2). Defining the Simulation Box and Solvation:
After generating the topology, we placed the AGO2 and miRNA complex in a defined simulation box using the editconf command. We used a rhombic dodecahedron geometry as it allowed us to maximize space efficiency. Once the box was set, we solvated the system by adding water molecules (genbox or solvate tool). This step is important to mimic the aqueous environment in which biological molecules naturally exist, ensuring realistic simulation conditions.

(3). Adding Ions to Neutralize the System:
To neutralize the overall charge of the system, we added counterions with the gmx genion tool, using a Verlet cutoff scheme (rcoulomb = 1.8 nm, rvdw = 1.8 nm) to maintain electrostatic stability during the molecular dynamics run. This step prevents artificial charge buildup, which could otherwise lead to unrealistic interactions in the simulation. Physiological salt concentrations can also be included to better replicate in vivo conditions.

(4). Energy Minimization:
Before running any molecular dynamics, energy minimization was performed using the mdrun command with a cutoff of 1.8 nm for short-range electrostatics and van der Waals interactions, ensuring that the system’s potential energy was lowered and no steric clashes or bad contacts remained. Energy minimization helps ensure that the system is in a stable state before equilibration and production MD.

(5). Equilibration:
Equilibration is a critical step to allow the system to relax and reach a stable state under simulated physiological conditions. It was performed in two phases: NVT (constant Number of particles, Volume, and Temperature) and NPT (constant Number of particles, Pressure, and Temperature) ensembles. During the NVT phase, temperature was maintained at 300 K using the V-rescale thermostat with a time constant of 0.1 ps. In the NPT phase, pressure control was handled by the C-rescale algorithm at 1 bar with a compressibility of 4.5e-5 bar⁻¹ to ensure the system mimicked physiological conditions.

(6). Production MD Simulation:
After equilibration, the system was ready for the production run, where the actual molecular dynamics simulations were carried out over 1000 ps with an integration step size of 1 fs. Data from the trajectory were output every 10 ps for further analysis. This step allows for observing the dynamic behavior of the AGO2-miRNA complex, capturing important interactions, conformational changes, and potential binding events under physiological conditions. The trajectory file generated during this phase is analyzed to extract meaningful insights, such as binding affinity or molecular motions.

(7). Results:
These steps were performed for Ago2, miRNA-326, and the docked version, allowing us to successfully stabilize the systems, setting the stage for further mutation analysis.

AGO2 in Water:

miRNA-326 in Water:

AGO2 - miRNA326 in Water:

(8). Technical Challenges and Learning:
Force Field Integration: Most GROMACS tutorials, including the one we followed, utilize the OPLS force field, which lacks parameterization for protein-nucleic acid interactions. Therefore, we had to implement the CHARMM36 force field, a newer version not prepackaged with GROMACS, enhancing the accuracy of our simulations.

Topology Building: We encountered challenges with GROMACS not detecting certain atoms, which necessitated manual editing of the PDB files to remove these atoms, ensuring accurate simulations. The topology generators such as CGenFF were not built for ligands like microRNAs. We had to use Automated Topology Builder, which wasn’t as compatible with the CHARMM forcefield. This caused misalignment in atoms and errors with structure.

GROMACS Compilation: Compiling GROMACS on Stony Brook University's High-Performance Computing (HPC) cluster was a significant challenge. This process required learning to use HPC resources and Slurm, essential for managing computational jobs efficiently.

HPC Utilization: The experience with HPC has facilitated our GROMACS simulations and enabled us to explore other computational tools, such as Monte Carlo simulations, for different aspects of our project.

Extending from the Main Project

By understanding the optimal binding between Ago2 and miRNA-326, we can understand how the sensor might detect the presence of miRNA-326 based on its binding to Ago2, allowing you to fine-tune sensor design to maximize sensitivity and specificity. We can use Pymol or Chimera to create Ago2 mutants, then run through MD simulations to observe stability and conformational changes, and finally use umbrella sampling to investigate the detailed binding mechanism and free energy profiles.

There are also therapeutic implications such as aiding in designing small molecules or therapeutic interventions that modulate this interaction. This could be valuable in cases where miRNA-based therapies are explored, such as miRNA mimics or anti-miRs.

Conclusion:

Our use of GROMACS for modeling AGO2 mutations is a strategic approach to enhancing miRNA binding, contributing to the broader field of synthetic biology. This work not only propels our iGEM project forward but also offers a robust framework for exploring protein-RNA interactions at the molecular level, supporting the optimization of our bacteria-detection system.