COMPUTATIONAL MODELING

INTRODUCTION

Computational modeling applies theoretical and mathematical frameworks to describe and predict the behavior of biological processes. In our project, these models are used to support the design, refinement, and prediction of our engineered CRISPR-Cas13a system. With the help of computational tools such as MATLAB, we can develop complex systems of reactions and derive systems of ordinary differential equations (ODEs) to reflect such biological processes. This approach allows us to simulate system dynamics, generate quantitative predictions, and help optimize experimental design ahead of wet lab testing.

Both deterministic and stochastic models were used to capture the dynamics of our system. Deterministic models use ODEs to describe the average behavior of the system under ideal conditions, while stochastic simulations, such as Gillespie's algorithm, capture random fluctuations that become significant at lower molecule counts. Together, both models provide a robust view of system dynamics where deterministic models identify key trends and stochastic simulations show potential variability and random noise that may be of importance.

General Model Assumptions:

  • All reactions are assumed to occur in a homogeneous, well-mixed environment.
  • All species start at a known, fixed initial concentration.
  • Concentrations are modeled on a continuous scale (ex. nM), assuming sufficient molecule numbers for deterministic behavior.
  • Binding interactions are assumed to occur with 1:1 stoichiometry unless otherwise specified.

DETERMINISTIC MODELING

Deterministic models give us the average or “expected” behavior of the system by assuming large numbers of molecules and ignoring random noise. Given initial concentrations and reactions, these models result in predictable trajectories of system dynamics, providing a baseline for how the system should behave under ideal conditions. This is useful for identifying key parameters, validating initial assumptions, and guiding experimental design.

IN-VITRO

Model Assumptions

  • Cas13a and gRNA binding is treated as irreversible. Literature values for the dissociation rate are very small, making unbinding negligible within the experimental timeframe (Sinan et al.).
  • RNP Complex degradation is excluded since this experiment is run over a relatively short duration (less than a day) and the effects of complex degradation are minimal and unlikely to alter system dynamics.
  • To further simplify the model, multi-step processes are put together into one single reaction for rapid equilibrium.
    • RNP Binding and Activation is represented as one step
    • RNP & HIV Binding and Cleavage is represented as one step

Cas13a binds to gRNA to form a ribonucleoprotein (RNP) complex:
Equation 11
The RNP complex binds to HIV RNA and cleaves it, activating the trans subunit of the protein and turning it into a hyperactive RNP. This step represents the first cis-cleavage of the target sequence that triggers collateral activity.
Equation 11
The hyperactive complex (RNP*) cleaves both HIV RNA and nonspecific RNA sequences nearby at an accelerated rate:
Equation 11
Degradation rate of gRNA, HIV, and non-specific RNAs:
Equation 11
From this system of reactions, the following ODEs were derived to describe how these concentrations change over time:
Equation
From literature, we obtained key rate constants such as degradation and synthesis rates. For degradation, we used reported half-lives of RNA or protein and converted them into rate constants using the formula:
Equation 11
This allows us to represent biological decay as a continuous rate in the model. Similarly, other parameters such as synthesis rates and binding affinities were derived or approximated from published experimental data and biological databases.

Parameter Description Value Unit Reference
kRNP Cas13a–gRNA RNP complex formation rate 2.9 × 10−3 nM−1·sec−1 (Sinan et al.)
kdeg,gRNA, kdeg,NS, kdeg,HIV Degradation rate of RNA in-vitro 0.0009 sec−1 (Shin and Noireaux)
kcis Cas13a–gRNA RNP complex activation and cleavage rate of HIV 8.3 × 10−4 nM−1·sec−1 (Feng et al.)
ktransNS Hyperactivation state of Cas13a–gRNA complex (RNP*) cleavage of nonspecific RNAs 2×, 3×, 4×, 5× of kcis nM−1·sec−1 (Feng et al.), (Nalefski et al.)
ktransHIV Hyperactivation state of Cas13a–gRNA complex (RNP*) cleavage of HIV RNA 3× of kcis nM−1·sec−1 (Feng et al.), (Nalefski et al.)
Table 1: Table defining rate constants used for our in-vitro Cas13a cleavage model.
Concentration of HIV RNA over time

Figure 1: Concentration of HIV RNA over time.

HIV concentration over time for different values

Figure 2: HIV concentration over time for different values of ktransHIV, scaled from 2x - 10x relative to cis-cleavage rates (kcis).

HIV RNA levels decline rapidly within ~500 seconds at baseline rates (ktransHIV = 3 x kcis), demonstrating the efficiency of Cas13a-mediated cleavage once activated. (Figure 1) By varying the hyperactivation rate (ktransHIV) and observing HIV RNA concentrations over time as ktransHIV increases from 2× to 10× relative to kcis, we find that HIV RNA clearance accelerates, confirming that higher trans-cleavage activity substantially enhances overall knockdown efficiency. (Figure 2)
Deterministic simulation showing HIV RNA concentration over time

Figure 3: Deterministic simulation showing HIV RNA concentration over time for different Cas13a:gRNA ratios (1:1, 1:2, 1:3, 1:5, 1:10) across three total Cas13a + gRNA concentrations (1.0 nM, 10.0 nM, 100.0 nM). Each subplot represents a different total concentration, with higher concentrations showing faster knockdown.

We observe that knockdown is faster when gRNA is in excess relative to Cas13a and at higher total concentrations. Both stoichiometry and molecule abundance influence cleavage efficiency in vitro. These results not only help identify optimal conditions for in-vitro Cas13a experiments but also provide the concentrations derived from our graphs, which were used to guide wet-lab testing.

IN-VIVO

Model Assumptions

In addition to the assumptions applied to our In-Vitro model,

  • We assume CD4 receptor expression on our HEK293 cells is already saturated and remains constant throughout the simulation. This reflects our experimental design where creating a stable cell line with our CD4 receptors will take ~ 1 month.
  • Cas13a delivery is modeled as two compartments, where Cas13ex is extracellular and Cas13in is intracellular. Transport from Cas13ex to Cas13in occurs via receptor-mediated uptake by CD4 receptors and is modeled as a first order reaction.

While the in-vitro model captures the core cleavage dynamics of Cas13a and HIV cleavage, it does not account for the additional complex factors present inside a cell. Building on the In-vitro framework, the in-vivo model introduces additional steps:

Extracellular Cas13a binds to the CD4 receptor (Y) and enters the cell. This process is modeled with Michaelis-Menten kinetics, where the binding saturates at high Cas13a concentrations. The half-saturation constant (Km = 100 nM) represents the extracellular Cas13a concentration at which the uptake rate is half of its maximum, given the receptor density.
Equation 11
Inside the cell, Cas13a binds to gRNA to form a ribonucleoprotein (RNP) complex:
Equation 11
Since we are now working with an in-vivo system, we consider that proteins and RNAs naturally degrade inside the cell over time, unlike the in-vitro system where degradation was negligible.
Equation 11
We also consider the synthesis of new gRNA, HIV RNA, and nonspecific RNAS:
Equation 11
These reactions remain the same from our in-vitro model, where the cis-cleavage of HIV activates Cas13a into the hyperactive RNP complex, which then cleaves both nonspecific RNA and additional HIV sequences at a faster rate.
Equation 11
From this system of reactions, the following ODEs were derived to describe how these concentrations change over time:
Equation 11

Parameter Description Value Unit Reference
Y Concentration/Density of Receptors ~100–200 nM (Wang et al.), (“Size of HEK (Human Embryonic Kidney) 293 Cell - Human Homo Sapiens - BNID 108918”)
Km Half-saturation constant 100 nM (Hoare)
kdeg,HIV Degradation rate of HIV in CD4 cells 0.000004 sec−1 (Zhou et al.)
kdeg,gRNA Degradation rate of gRNA 0.00003 sec−1 (Saifullah et al.)
kdeg,NS Degradation rate of NS 0.00002 sec−1 (“Typical MRNA Half Life in Human Cells - Human Homo Sapiens - BNID 104747”)
kbind Binding affinity of CD4 to general antibody 4 × 10−4 nM−1·sec−1 (Kim et al.)
kRNP Cas13a:crRNA complex formation & activation 0.25 sec−1 (“Team:Newcastle/Model - 2019.Igem.org”)
kdeg,cas13in Degradation rate of internal Cas13 0.00003 sec−1 (Zhou et al.)
kdeg,cas13ex Degradation rate of external Cas13 0.000001 sec−1 (Wan et al.)
ksyn,gRNA, ksyn,HIV, ksyn,NS Synthesis of gRNA and HIV under CMV promoter 0.002 mRNA/sec (“Transcription Termination Rates - Mammalian Tissue Culture Cell - BNID 106381”)
kcis Cas13a-gRNA RNP complex activation and cleavage rate of HIV 8.3 × 10−4 nM−1·sec−1 (Feng et al.)
ktransNS Hyperactivation state of Cas13a-gRNA complex (RNP*) cleavage of nonspecific RNAs 5x of kcis nM−1·sec−1 (Feng et al.), (Nalefski et al.)
ktransHIV Hyperactivation state of Cas13a-gRNA complex (RNP*) cleavage of HIV RNA. 3x of kcis nM−1·sec−1 (Feng et al.), (Nalefski et al.)
Table 2: Table defining rate constants used for our in-vivo Cas13a uptake model.
Extracellular Cas13a (blue) binding to CD4 receptors

Figure 4: Extracellular Cas13a (blue) binding to CD4 receptors expressed on HEK293 cells and intracellular Cas13a accumulation (green) in the cell.

After reaching a peak, intracellular Cas13a gradually decreases due to degradation and being consumed by the Cas13a-gRNA RNP complex, while extracellular Cas13a approaches zero.

STOCHASTIC MODELING

While deterministic models work well for describing the average behavior of a system, they ignore the random noises and fluctuations bound to be present in a realistic biological environment. In living cells, molecule numbers can be very small, and intrinsic noise can strongly influence system behavior. To account for this, we utilized the Gillespie algorithm, a stochastic simulation algorithm (SSA), to simulate our system over time as a series of discrete, probabilistic molecular events.

Our stochastic model was developed in MATLAB using SSA. Unlike deterministic models that treat concentrations as smooth and continuous, this program simulates individual molecular reactions.

RATE SCALING:

To convert our deterministic model into a stochastic one, all concentrations and reaction rates were scaled to reflect molecular behavior. HEK293 cells have a total volume of approximately 700 fL to 1.8 pL (BNID 108918). Using the full cell volume in this simulation would result in high molecule counts, not reflecting the effects of intrinsic noise. To better represent random fluctuations, our model uses a 1 fL subvolume to represent an intracellular region.

Each species’ initial concentration (nM) was converted to the number of molecules present in a 1 fL volume using:

Equation 11
Where:
  • C = concentration (M)
  • V = system volume (L)
  • NA = Avogadro's number (6.022 * 1023 molecules/mol)
Bimolecular reactions depend on the number of possible molecule to molecule collisions. This was scaled from concentration-based ODE rates (nM/s) to molecule-based propensities using:
Equation 11
Unimolecular rates such as decay rates have the units of s-1 and do not need to be scaled (Lecca).

To make sure simulations ran properly, we avoided true zero values since SSA requires at least one molecule for a reaction to occur. After simulations, outputs were scaled back into concentrations for direct comparisons with our deterministic results.
50 stochastic simulation runs of HIV RNA cleavage

Figure 5: 50 stochastic simulation runs of HIV RNA cleavage by Cas13a-gRNA complexes in-vitro.

We find that HIV RNA levels steadily decrease over time as the Cas13a-gRNA complex works to cleave their target. Each colored line represents one possible trajectory, and the variability between them highlights the random nature of molecular interactions. Despite this, we see an overall trend, confirming consistent and effective HIV RNA knockdown (Figure 5).

STEADY STATE ANALYSIS

Steady state analysis in mathematical biology is a crucial part of analyzing system behavior when there is the constant production of one of more variables. Steady states occur when the rate of production of a variable is balanced with its rate of degradation, resulting in a constant level of expression (Steady State Kinetics). Knowing when steady states occur in a system and what their level of stability is allows for greater insight into the conditions which allow for homeostasis and what variables can cause disruption. It is of particular interest for us to know what steady states, if any, exist in our in-vivo system, as this will allow us to model realistic homeostatic conditions and understand which variables are responsible for potential instability.

To find potential steady states, each differential equation in our in-vivo system was set equal to zero (indicating no change in concentration) and solved simultaneously. One steady state was found:

Equation 11
Steady states can be stable, unstable, or partially stable. In stable steady states, the biological system returns to steady state after experiencing small perturbations, while in unstable steady states, the system does not (Tyson and Novak). In partially stable steady states, the system is stable with respect to some states but not others, resulting in varied behavior (W. M. Haddad and Somers). To determine the stability of our steady state, we performed stability analysis.

Due to the complex nonlinear nature of our in-vivo system, it was necessary to create a linear approximation at steady state in order to analyze stability. To do so, we computed a Jacobian matrix and evaluated it at our system’s steady state. A Jacobian matrix is a matrix consisting of the gradient of each function in a system, allowing for the behavior of that system to be examined at a specific point. The gradient of a function f is the vector-valued function Δf(x) which contains all first-order partial derivatives of that function with respect to inputs x1 … xn :
Equation 11
Each partial derivative represents the change in f that can be attributed to the input xi, where i = 1, …, n. In our system, there are seven different inputs, so each gradient takes the following form:
Equation 11
Each row of the Jacobian matrix contains the gradient of the i-th function in the system. Our functions are each ODE in the in-vivo system:
Equation 11
Therefore, the Jacobian matrix of our in-vivo system is the transpose of the gradient of each function f1, ..., f7:
Equation 11
Equation 11
where
Equation 11
The numerical values of each variable at steady state are substituted into the Jacobian matrix to compute its eigenvalues. If all eigenvalues have negative real parts, the steady state is stable, but if any eigenvalue has a positive real part, it is unstable. A system with no positive eigenvalues but one or more eigenvalues equal to zero is not fully stable and requires further analysis (Roussel).

The Jacobian matrix of our system computed at steady state has 6 negative eigenvalues and one equal to zero, indicating the steady state is not fully stable and requires further analysis. Examining the steady state, it occurs only when there is no Cas13a present in the system, resulting in values of zero for Cas13ex, Cas13in, RNP, and RNP*. To further examine the stability of this steady state, we modeled it with differing initial concentrations of HIV, NS, and gRNA in proximity to their respective numerical steady state values.
Concentration of HIV RNA over time

Figure 6: Varying initial concentrations of HIV, non-specific RNA, and gRNA all converge to steady state values in the absence of Cas13a in-vivo.

Since the concentrations of HIV, NS, and gRNA all converge to their steady state values regardless of their initial concentration, the steady state is stable in the absence of Cas13a. It cannot be considered completely stable in the context of the entire system, however, since the introduction of Cas13a causes the system to be unable to return to this steady state due to RNP complex formation, cis cleavage, and trans cleavage. This is what the eigenvalue equal to zero indicates.

To further explore the effect of Cas13a on our in-vivo system, we modeled the effect of various Cas13a initial concentrations on HIV knockdown. We utilized the numerical steady state values as initial concentrations of HIV, NS, and gRNA to simulate Cas13a entry into a HEK293 cell at steady state.
Concentration of HIV RNA over time

Figure 7: The effect of varying initial concentrations of Cas13a on the knockdown efficiency of HIV RNA in-vivo.

These results indicate that even small amounts of Cas13a are highly effective in reducing HIV concentration, with reduced efficiency enhancement beyond roughly 50 nM. This is due to the production of gRNA not being fast enough to quickly form more RNP complexes. These results show high efficacy of our therapeutic and will be crucial in guiding wet lab’s future in-vivo HIV cleavage experiments.

REFERENCES

  1. Feng, Wei, et al. “A Sensitive Technique Unravels the Kinetics of Activation and Trans‐Cleavage of CRISPR‐Cas Systems.” Angewandte Chemie International Edition, vol. 63, no. 22, 25 Mar. 2024, https://doi.org/10.1002/anie.202404069.
  2. Haddad, Wassim M, and Luke Somers. “Partial Stability of Nonlinear Dissipative Feedback Systems.” 2022 American Control Conference (ACC), 31 May 2023, ieeexplore.ieee.org/document/10156296, https://doi.org/10.23919/acc55779.2023.10156296.
  3. Hoare, Sam R. J. “The Problems of Applying Classical Pharmacology Analysis to Modern in Vitro Drug Discovery Assays: Slow Binding Kinetics and High Target Concentration.” SLAS Discovery: Advancing the Science of Drug Discovery, vol. 26, no. 7, 11 June 2021, pp. 835–850, https://doi.org/10.1177/24725552211019653.
  4. Kim, Youn H, et al. “Clinical Efficacy of Zanolimumab (HuMax-CD4): Two Phase 2 Studies in Refractory Cutaneous T-Cell Lymphoma.” Blood, vol. 109, no. 11, 1 June 2007, pp. 4655–4662, https://doi.org/10.1182/blood-2006-12-062877.
  5. Lecca, Paola. “Stochastic Chemical Kinetics.” Biophysical Reviews, vol. 5, no. 4, 30 July 2013, pp. 323–345, www.ncbi.nlm.nih.gov/pmc/articles/PMC5425731/, https://doi.org/10.1007/s12551-013-0122-2.
  6. Nalefski, Eric A., et al. “Kinetic Analysis of Cas12a and Cas13a RNA-Guided Nucleases for Development of Improved CRISPR-Based Diagnostics.” iScience, vol. 24, no. 9, Sept. 2021, p. 102996, https://doi.org/10.1016/j.isci.2021.102996.
  7. Roussel, Marc. Stability Analysis for ODEs. 2005.
  8. Saifullah, et al. “Effective RNA Knockdown Using CRISPR-Cas13a and Molecular Targeting of the EML4-ALK Transcript in H3122 Lung Cancer Cells.” International Journal of Molecular Sciences, vol. 21, no. 23, 24 Nov. 2020, www.ncbi.nlm.nih.gov/pmc/articles/PMC7727695/, https://doi.org/10.3390/ijms21238904.
  9. Shin, Jonghyeon, and Vincent Noireaux. “Study of Messenger RNA Inactivation and Protein Degradation in an Escherichia Coli Cell-Free Expression System.” Journal of Biological Engineering, vol. 4, no. 1, 2010, p. 9, https://doi.org/10.1186/1754-1611-4-9.
  10. Sinan, Selma, et al. “Kinetic Dissection of Pre-CrRNA Binding and Processing by CRISPR-Cas12a.” 25 July 2023, pmc.ncbi.nlm.nih.gov/articles/PMC10402064/, https://doi.org/10.1101/2023.07.25.550589.
  11. “Size of HEK (Human Embryonic Kidney) 293 Cell - Human Homo Sapiens - BNID 108918.” Harvard.edu, 2025, bionumbers.hms.harvard.edu/bionumber.aspx?s=n&v=1&id=108918.
  12. Steady State Kinetics in Regulating Homeostasis of Cells. 16 Mar. 2023, https://chem.libretexts.org/@go/page/418917.
  13. “Team:Newcastle/Model - 2019.Igem.org.” Igem.org, 2019, 2019.igem.org/Team:Newcastle/Model#ref-2.
  14. “Transcription Termination Rates - Mammalian Tissue Culture Cell - BNID 106381.” Harvard.edu, 2025, bionumbers.hms.harvard.edu/bionumber.aspx?id=106381&ver=4&trm=CMV+promoter+transcription+rate&org=.
  15. “Typical MRNA Half Life in Human Cells - Human Homo Sapiens - BNID 104747.” Harvard.edu, 2025, bionumbers.hms.harvard.edu/bionumber.aspx?id=104747&ver=12&trm=mRNA+decay&org=.
  16. Tyson, John, and Bela Novak. Cell Cycle Vignettes Bistability and Oscillations. 2011.
  17. Wan, Yiming, et al. “Optimizing a CRISPR-Cas13d Gene Circuit for Tunable Target RNA Downregulation with Minimal Collateral RNA Cutting.” ACS Synthetic Biology, vol. 13, no. 10, 2024, pp. 3212–3230, pubmed.ncbi.nlm.nih.gov/39377757/, https://doi.org/10.1021/acssynbio.4c00271.
  18. Wang, Meiyao, et al. “Quantifying CD4 Receptor Protein in Two Human CD4+ Lymphocyte Preparations for Quantitative Flow Cytometry.” Clinical Proteomics, vol. 11, no. 1, Dec. 2014, https://doi.org/10.1186/1559-0275-11-43.
  19. Zhou, Yan, et al. “Kinetics of Human Immunodeficiency Virus Type 1 Decay Following Entry into Resting CD4+ T Cells.” Journal of Virology, vol. 79, no. 4, 15 Feb. 2005, pp. 2199–2210, https://doi.org/10.1128/jvi.79.4.2199-2210.2005.