Introduction to Modelling SHERLOCK Assay Dynamics

---

Accurate modelling of the SHERLOCK detection assay is essential for optimizing its performance and interpreting experimental data effectively. This assay, which uses CRISPR-Cas13a technology, benefits from two primary modelling approaches:

1. Detailed Law of Mass Action Model: This model accounts for all relevant elements and their interactions within the SHERLOCK assay. It uses ordinary differential equations to describe how the Cas13a-crRNA complex binds to the target RNA, becomes activated, and cleaves fluorescent reporters. This comprehensive approach provides a detailed view of the assay's dynamics, crucial for validating experimental results.
2. Simplified Michaelis-Menten Model: To streamline analysis and reduce computational complexity, we use the Michaelis-Menten kinetics framework. This model focuses on the rate-limiting trans-cleavage activity, simplifying the system to key parameters like the catalytic turnover rate ($ k_{\text{cat}} $) and Michaelis constant ($ K_M $). It helps predict the assay's timescale and optimize conditions for accurate and efficient detection.

The model is validated through back-of-the-envelope checks, ensuring that kinetic parameters align with expected physical limits. These models not only enhance our understanding of the SHERLOCK assay but also guide its practical application and refinement.

Table 1. A short summary of modelling and calculations that were done.

Developing an effective and accurate model

---

Goal: To develop a model that reliably represents the SHERLOCK reaction kinetics to hone some practical aspects of conducting the test.

Model based on the Law of Mass Action

At the beginning a model involving all the elements involved in the SHERLOCK detection assay was considered, as illustrated in Figure 1.

Figure 1.The dynamics of the SHERLOCK detection assay action. The Cas13a-crRNA complex binds to the target RNA sequence with a constant $ k_1 $. The activated complex exhibits non-specific (“collateral”) endonuclease activity, cleaving the RNA probes with a constant $ k_2 $ which induces fluorescence. Cas13a can also cut specifically at the target sequence with a constant $ k_3 $. The backward reaction rates are omitted due to being negligible.

At first, the Cas13a-crRNA complex binds to the target RNA sequence from Prymnesium parvum (which is previously transcribed from DNA) with a kinetic constant of $ k_1 $. Subsequently, the Cas13a starts exhibiting endonuclease activity. Typically for itself, the Cas13a cuts both the target RNA and non-specifically the RNAs in the environment. Upon being cut, the fluorescent RNA reporter (in case of our test – RNaseAlert) exhibits fluorescence. This event proceeds with a kinetic constant $ k_2 $. The specific cut at the target sequence proceeds with a kinetic constant $ k_3 $. The kinetic constants of the backward reactions are non-zero, however, in each case the equilibrium is far-shifted towards the right hand side which makes the constants of the backward reactions negligible and computationally unnecessary. Consequently, the ordinary differential equations describing the dynamics of the SHERLOCK detection assay are as follows:

Where $\text{GC}$ is the Cas13a-crRNA complex, $\text{GCT}$ is the Cas13a-crRNA complex bound to the target sequence, $ P_a $ is the cut reporter which exhibits fluorescence, $ P_i $ is the uncut reporter, $ T_c $ is the cut target sequence, and $ T_{\text{uc}} $ is the uncut target sequence. Using the example of equation (1) it can be seen that the concentration of the Cas13a-crRNA complex not bound to the target RNA sequence decreases as the reaction with kinetic constant $ k_1 $ proceeds, and increases with a progress of the reaction with kinetic constant $ k_3 $.

These dynamics can be plotted as shown in Figure 2, where the Equations (1)-(6) were numerically solved for specified time points. This operation was performed by using the solve_ivp function which is a part of the SciPy Python library. 97 time points were selected to be within the range of 0-90 minutes. The initial conditions for the simulations and their sources are stated in Table 1.

Table 2. The initial conditions and kinetic constants used to solve the system of ordinary differential equations describing the dynamics of the SHERLOCK detection assay action (1-6).

The figure below shows the plots of the solved equations.

Figure 2. The solutions of equations (1) - (6) using the kinetic constant values and initial conditions from Table 2. Abbreviations as in the equations.

Below, we show the plot of the experimental data from our test for synDNA and SynCrRNA sequences provided by Kellner et al. [2]

Figure 3. Test Sherlock assay fluorescence readings for synDNA + syncrRNA sequences from Kellner et al. [2].

The solution to equation (3) represents the change in the cleaved fluorescent probe which corresponds to an increase in fluorescence. However, there are two main challenges in converting the fluorescence data to the concentration of cleaved RNaseAlert: first, the composition of the RNaseAlert probe is unknown, making it difficult to express its concentration in moles; second, there's no established relationship between fluorescence intensity and the concentration of cleaved probe. Despite this, the similar shape of the charts included in Figures 2 (3) and 3 suggest that the approach used to model the SHERLOCK assay dynamics is likely correct.

Overall, the model based on the law of mass action allows to relate the chemical reaction rate to concentration. It gives us a comprehension of the changes in concentration of all reactanst present in the reaction which is useful in gaining an understanding of the overall dynamics of the system.

Michaelis-Menten kinetics

Although the model described in the previous section provides valuable insights into the dynamics of SHERLOCK detection assay activity, it involves many parameters, which can obscure key interactions. Some of these interactions have a direct impact on the assay's response time and the optimization of test sample and fluorescent reporter quantities. Additionally, the model comes with a significant computational cost. Therefore, simplifying the model is a reasonable approach.

Following the line of assumptions presented in [3], it is possible to reduce the model to a simple form of a Michaelis-Menten equation. Reframing the understanding of the dynamics into the Michaelis-Menten kinetics, a model for trans-cleavage activity in CRISPR diagnostics involves two main reactions: a cis-cleavage of the target RNA and a trans-cleavage reaction of the RNAs in the environment. In the cis-cleavage reaction, the target nucleic acid activates the CRISPR enzyme. Studies indicate that cis-cleavage happens rapidly even at low target concentrations, with reaction times around 100 seconds [3]. Therefore, the assumption is made that cis-cleavage is not the rate-limiting step in diagnostic applications. As a result, the initial enzyme concentration, $ [\text{GCT}_0] $ now denoted as $ E_0 $, is assumed to be approximately equal to the concentration of the target nucleic acid, under the premise that all target molecules have reacted quickly with the abundantly available enzyme molecules.

The trans-cleavage activity (non-specific endonuclease activity), which is the focus of the model, is a multiple-turnover enzymatic reaction catalysed by the target-activated enzyme. This reaction is typically the rate-limiting step for generating a detectable signal in CRISPR diagnostics. The trans-cleavage reaction can be modelled using the Michaelis-Menten equation, which describes the kinetics of enzyme-mediated reactions. In this model, $ E $ represents the target-activated Cas13a-crRNA complex, $ S $ (for substrate) – uncleaved reporter, $ C $ is the enzyme-reporter complex, and $ P $ represents the cleaved reporter molecules. The rate constants for the reaction include $ k_f $ (forward rate), $ k_r $ (reverse rate), and $ k_{\text{cat}} $ (catalytic turnover rate). Notice that the endonuclease activity of enzyme at the target site is omitted. This is due to it being slow compared to non-specific activity [1].

To simplify the model, the quasi-steady-state assumption is used. This assumption states that the concentration of the intermediate enzyme-substrate complex remains nearly constant over time. Consequently, the concentration of this complex can be approximated, allowing the product formation rate (reporter cleavage) to be described by the Michaelis-Menten equation. Graphically the reduction of the model can be represented by Figure 4.

Figure 4. Illustration of the SHERLOCK detection assay activity reduced to Michaelis-Menten kinetics.
$ k_{\text{cat}} $ is equal to $ k_2 $.

Overall, the refined model assumes that cis-cleavage is not the rate-limiting step and focuses on the trans-cleavage activity, applying the Michaelis-Menten kinetics framework to simplify and describe the reaction dynamics. This allows for effective prediction of the kinetics of the SHERLOCK detection assay activity. Following the line of assumptions presented in this section we ensure that such a model can be a valid closed-form approximation for reporter cleavage versus time.

Fitting the model with kinetic data

The most suitable literature data seem to be the kinetic constants reported in [1] for the following set of reaction reactants: Cas13a, mature crRNA, “Target 1”, RNaseAlert trans-substrate. This is due to the presence the same CRISPR enzyme and reporter as used in our experiments, as well as the use of mature (not precursor) crRNA. The $ k_{\text{cat}} $ is equal to 27.1 s^-1, $ K_M $ of 294 nM, and $ E_0 $ can be approximated by the target concentration which is equal to 10 pM. Figure 5 is a Michaelis-Menten plot for these values.

Figure 5. Michaelis-Menten kinetics plot using kinetic constants from [1]. $ k_{\text{cat}} $ = 27.1 s-1, $ K_M $ = 294 nM

Using the literature data presented above Equation 8 can be solved for the concentration of cleaved reporter concentration over time. Figure 6 represents the results of the calculations.

Figure 6. The model based on Michaelis-Menten kinetics solved for cleaved reporter (RNaseAlert) concentration over time. The left plot is for the initial reporter concentration which is higher than the $ K_M $, plot on the right – lower than $ K_M $.

The left plot is for the initial reporter concentration of 2 μM, typically used in our experiments. It is higher than the $ K_M $ = 294 nM. This causes the reaction rate to approximate to zero-order kinetics. In this regime, the rate of substrate depletion (cleavage of the reporter) becomes approximately constant.

The plot on the right represents the solution for the initial reporter concentration of 2 nM which is lower than $ K_M $.

Comparing these results with Figure 3 which contains wet-lab results, it can be deducted that the $ K_M $ for the crRNA/target pair utilized in the experiments must have been higher than the initial reporter concentration and therefore different than the value from [1]. It is important to remember that even for the same CRISPR enzyme, the $ k_{\text{cat}} $ and $ K_M $ catalytic parameters may differ depending on the sequence of crRNA and the target sequence.  It is thus recommended that kinetic measurements should be performed for various target/crRNA target pairs in order to select a pair with the highest catalytic efficiency which may lead to the improvement of a low limit of detection (LOD) and increase the assay speed for the same LOD [3].

Model check

---

Back-of-the-envelope checks for kinetic constants

In our mathematical model, we applied back-of-the-envelope checks described in [3] to assess the accuracy of published Michaelis−Menten kinetics data for the CRISPR-Cas13a enzyme. These checks serve as quick validations to ensure that the reported enzyme parameters do not violate basic physical principles. The checks focus on three key parameters: $ \alpha $, $ \beta $ and $ \gamma $.

$ \alpha $ is a dimensionless parameter that compares the amount of cleaved reporters at time $ t_{\text{lin}} $ to the initial concentration of the uncleaved reporters $ S_0 $. It is calculated as:

where $ v $ is the reported reaction velocity and $ t_{\text{lin}} $ is time, when fluorescence intensity increase to time is linear. If $ \alpha > 1 $, this suggests that product formed during the initial linear portion of the kinetics measurements exceeds the initial amount of substrate, which is not realistic.

$ \beta $ compares the enzyme’s maximum reaction velocity to its theoretical maximum. It is calculated as:

where $ v_{\max} $ is the theoretical maximum rate of product formation when the substrate concentration is significantly greater than $ K_M $ , and is given by $ k_{\text{cat}} \cdot E_0 $. If $ \beta > 1 $ , the reported values suggest that the enzyme can reach catalytic efficiency beyond realistic limits for the concentration of substrate provided. This means that the enzyme seems to perform better than physically possible, implying an error in the reported kinetic parameters.

$ \gamma $ tests if the time scale of linear portion of the kinetics data $ t_{\text{lin}} $ is at most on the same order as the reaction time scale $ \tau $.

Results for gathered data

From data available in [1], $ v $ can be estimated to be equal to 12 nM/min, $ t_{\text{lin}} $ is approximately 10 min, and $ S_0 $ is approximately equal to 1000 nM.

For these values, and previously cited $ k_{\text{cat}} $ and $ K_M $, the $ \alpha $ and $ \beta $ parameters were less than one meaning the first two back-of-the-envelope checks were passed.

However, the calculated $ \gamma $ was of the order of magnitude of thousands compared to the expected order of one. This suggests that the enzyme is operating in a way that is inconsistent with the expected kinetics or alternatively, some procedural errors or incorrect estimations were made, necessitating a need for additional or improved measurements.

Practical aspects

---

Timescale of the SHERLOCK reaction

Goal: To assess time scale needed to obtain reliable results of SHERLOCK assay

The time scale $ \tau $ for significant product formation is equal to:

$ \tau $ governs the characteristic time to approximately complete the reaction.

It is:

• directly proportional to the $ K_M $ which represents a higher affinity of the activated enzyme to the reporter
• inversely proportional to $ k_{\text{cat}} $ related to turnover rate and the start SHERLOCK assay complex concentration, $ E_0 $.

Using the example of kinetic constants for Cas13a, mature crRNA, Target 1, RNaseAlert Trans-substrate from [1] which seem the most applicable for our project, we can approximate the timescale for the SHERLOCK reaction.

With $ k_{\text{cat}} $ of 27.1 s^-1, $ K_M $ of 294 nM, and target concentration of 10 pM approximately equal to $ E_0 $, $ \tau $ is approximately 18 minutes, which seems to be at a right order of magnitude for a 1-hour reaction.

It is, however, important to remember that since the $ K_M $ for the experiments on the PrymDetect system was most likely higher than 294 nM (mentioned in the section above), the reaction time scale for our reaction should be longer, which seems to be reflected in our wet-lab experiments. It is also important to point out that the assumption that the target sequence concentration is equal to $ E_0 $ is not entirely applicable since not all of the RNA target is utilized to activate the Cas13a-crRNA complex [4]. This lowers the value of $ E_0 $ used for calculations further elongating the reaction time scale.

The dynamics of RPA and IVT processes

Goal: Better understand the connections between RPA and IVT.

Recombinase Polymerase Amplification (RPA) is a rapid and efficient DNA amplification technique. Unlike PCR (Polymerase Chain Reaction), RPA operates at a constant, moderate temperature, typically between 37°C and 42°C, making it more accessible for field applications, which is why it was chosen for our project. Despite the advantages of this technique, its quick saturation hinders accurate quantitation [5]. Although it is claimed that at low primer concentrations the reaction does not saturate [5], our experiments do not support this notion suggesting a power-law increase for a wide range of both primer and target concentrations (see the “Measurement” page for a detailed description of our experiments). These experiments, however, can be subjected to further optimization, in particular by changing the dilution of the RPA kit [6]. In the particular case of our project, a T7 polymerase promoter was also added during the RPA process (by utilizing specifically designed forward primer) in order to enable in vitro transcription of the amplified DNA of Prymnesium parvum parallelly to the SHERLOCK reaction.

In vitro transcription (IVT) is the synthesis of RNA from a DNA template outside of the living cells. In this process, a DNA template with a specific promoter sequence, such as the T7 promoter, is used to guide RNA polymerase, like T7 RNA polymerase, to initiate transcription. This reaction proceeds in a linear fashion [7].

The dynamics of the RPA and IVT processes are illustrated in Figure 7. The amount of the initial Prymnesium parvum DNA remains unchanged (is a constant function) since it serves as a template in the RPA reaction. As the amplification progresses a new DNA sequence, DNA2, consisting of the initial Prymnesium parvum DNA with added T7 polymerase promoter is synthesized. This process follows an exponential growth with saturation [5]. Gompertz function, modelling slow initial growth, rapid middle growth, and a decelerating final phase that asymptotes to a maximum value was used to visualize it. This function is often used to describe biological growth [8].

Where: $ B(t) $ is the concentration of DNA2 at time t, $ B_{\text{max}} $ is the maximum possible concentration of DNA2, $ \mu $ is the growth rate, $ \lambda $ is the lag time or shift along the time axis (to account for delayed start), $ e $ is Euler's number (approximately 2.718).

The progression of IVT i.e. the increase in RNA concentration, can be modelled using a simple linear function. The RPA stage lasts for 30 minutes while the SHERLOCK reaction combined with IVT is incubated for 60 minutes, hence the increase in RNA concentration begins at the 30-min timepoint.

Figure 7. A plot visualizing the progress of RPA and IVT reactions over time. There is no scale on the y-axis since the plot serves solely as an overview of the progression.

Conclusions

Although sufficient data was not publicly available nor collected during the project, Figure 7 demonstrates an overview of progression the reactions, enabling for a clearer understanding of these processes which are an essential part of the PrymDetect. Since CRISPR-based systems differ in the type of Cas enzyme, target sequence (ssDNA or RNA), preamplification methods, and the necessity for in vitro transcription of the target sequence, such visualization aids the comprehension of this particular system designed to fit the needs of Prymnesium parvum detection.

It is important to remember that not all of the RNA target is utilized to activate the Cas13a-crRNA complex [4]. Hence, the RNA concentration does not directly translate into the magnitude of the fluorescence caused by the endonuclease activity of the Cas13a-crRNA-target RNA assay. Furthermore, the increase in fluorescence due the endonuclease activity and increase in RNA concentration due to the IVT reaction happen simultaneously and at different rates making them difficult to disentangle and mathematically model.

Our approach is focused on qualitative intuition enabling for grasping the key concepts of the enzyme dynamics in the PrymDetect system. However, additional and extensive experiments involving determination of multiple kinetic constants and calibration curves, as well as optimization of experimental parameters would need to be performed in order to determine the behaviour of the system in a fully quantitative manner.

References

[1] Nalefski, E. A., Patel, N., Leung, P. J., Islam, Z., Kooistra, R. M., Parikh, I., ... & Madan, D. (2021). Kinetic analysis of Cas12a and Cas13a RNA-Guided nucleases for development of improved CRISPR-Based diagnostics. Iscience, 24(9).

[2] Kellner, M. J., Koob, J. G., Gootenberg, J. S., Abudayyeh, O. O., & Zhang, F. (2019). SHERLOCK: nucleic acid detection with CRISPR nucleases. Nature protocols, 14(10), 2986-3012.

[3] Ramachandran, A., & Santiago, J. G. (2021). CRISPR enzyme kinetics for molecular diagnostics. Analytical chemistry, 93(20), 7456-7464

[4] Fozouni, P., Son, S., de León Derby, M. D., Knott, G. J., Gray, C. N., D’Ambrosio, M. V., ... & Ott, M. (2021). Amplification-free detection of SARS-CoV-2 with CRISPR-Cas13a and mobile phone microscopy. Cell, 184(2), 323-333.

[5] Gootenberg, J. S., Abudayyeh, O. O., Kellner, M. J., Joung, J., Collins, J. J., & Zhang, F. (2018). Multiplexed and portable nucleic acid detection platform with Cas13, Cas12a, and Csm6. Science, 360(6387), 439-444.

[6] Valloly, P., & Roy, R. (2022). Nucleic acid quantification with amplicon yield in recombinase polymerase amplification. Analytical Chemistry, 94(40), 13897-13905.

[7] Marras, S. A., Gold, B., Kramer, F. R., Smith, I., & Tyagi, S. (2004). Real‐time measurement of in vitro transcription. Nucleic acids research, 32(9), e72-e72.

[8] Winsor, C. P. (1932). The Gompertz curve as a growth curve. Proceedings of the national academy of sciences, 18(1), 1-8.