Our Dynamics Modeling section ran throughout the project. Corresponding to the four parts of the project design, we built four sub-models: reverse transcription isothermal amplification model, enzyme kinetic model for trans cleavage of CRISPR-Cas system, Model of starch hydrolysis by γ-amylase and Electrical signal conversion model for glucose concentration. Each sub-model was connected using upstream sub-model outputs as the downstream sub-model inputs. Finally, we simulated the kinetics behavior of the entire detection system, based on ODE equations and experimental data. Simulation output and wet lab experiment data reached a R-squared greater than 0.98. Mathematical models were solved using R as a software tool.
The simplified linc00857 reverse transcription process is shown in the following equation. The reaction process is rapid relative to the DNA isothermal amplification process and can be approximated as the amount of RNA added being reverse transcribed to the amount of DNA template in a very short time.
For this section, we drew on the RPA model from the dry experiments of the 2019 Thessaly team. Based on the kinetic equations in the literature [1], we construct models by simplifying the reaction mechanism and changing some parameter values in conjunction with the data from the experiment.
i. There is perfect symmetry between reactions involving forward and reverse primers.
ii. ATP depletion was negligible over the time span examined.
iii. The fluid is well mixed and all reactions that occur are kinetically controlled.
At a constant temperature of 37°C, single-stranded DNA-binding proteins (GP32s) first binds to the primer, and then recombinases replace them in the primer DNA to form a polymer. When the polymer searches for a sequence that is completely complementary to the primer on the template DNA, the template DNA melts with the help of GP32s. At the same time, under the action of DNA polymerase, the nucleotides in the system are consumed to form a new DNA complementary chain. This process is repeated continuously to achieve exponential amplification of DNA numbers.
Compared with the reaction formula given in the literature [1], we simplified the four-step reaction formula in which n recombinases(nR) sequentially substitute m GP32(mG) for primer binding to a one-step reaction formula (2). And the reaction rate constant for this step was adjusted according to the literature.
There are chemical equations as follows, in which substance names and parameter names can be found in Appendix 1.
Ordinary differential equations describing the RPA process are shown below, where the substance names and parameter names can be found in Appendix 1.
Based on the actual reaction time and RNA input of RT-RPA for the wet experiment, and the assumptions of the reverse transcription step, we took the RNA input of 1E-12 M as the input of DNA. The DNA concentration at 900s was taken as the output of the RPA step.
The following DNA amplification concentration profiles are shown. At 900s, the value of DNA is 2.97E-7.
Kinetics of the rest relevant variables in RPA model are shown below.
Since the Cas-gRNA-specific activation step is strongly influenced by the base sequence, and it is often performed together with the trans-cleavage activity verification experiment during the wet-lab reaction. Therefore, we treat this process in the same way as the reverse transcription process. It is assumed that the step is very fast, i.e. the amount of DNA input in the previous step can be approximated as the amount of Cas-gRNA-activator complex generated in the next step.
Based on the literature [4], we quantified the kinetic parameters of Cas12a and Cas14a assuming that the trans-cleavage process of Cas proteins conforms to the Mie equation for enzyme kinetics, and we treated the activated Cas-gRNA-activator complex as an enzyme.
Where V is the rate of production of the cleaved reporters, Vmax is the highest measured reaction velocity for the experiment, KM is the Michaelis constant, [S] is the concentration of uncleared substrate at equilibrium.
i. Fluorescence value unit conversion
To achieve comparability of data, we converted the fluorescence units (AU) measured to concentration units (nM) using the formula in the literature [5] for the calibration of fluorescence values to reporter gene concentrations.
ii. Fluorescence signal fitting Fitting the wet experimental data after concentration unit conversion, we obtained the changes in cleaved reporter concentrations for the Cas12a and Cas14a over 180 min.
a. Cas12a
b. Cas14a
iii.Fitting of enzyme kinetic equations
According to the above fitted trend of the concentration of the cleaved reporter substrate, the production rate of the cleaved reporter substrate at the beginning of the reaction is basically linear. Combining with previous literature, we selected the data of the first two cycles (Cas12a, about 4 min) and the first three cycles (Cas14a, about 6 min) for the calculation of the initial rate. And then we fitted the initial velocities of the Cas12a and Cas14a at the ssDNA reporter substrate concentration gradient by the Michaelis–Menten kinetics and Linear fitting, respectively. The following are the fitting results.
c. Cas12a for the first two cycles
d. Cas14a for the first three cycles
The results of fitting the above two sets of data are summarized as follows.
Combined with the above experimental data, we modeled the kinetics of Cas12a and Cas14a separately in this experiment. The corresponding reaction differential equations are shown below, and the reaction rate constants were obtained from the fitting results.
We assumed that the rate of γ-amylase production was proportional to the rate at which ssDNA was cleaved in the DNA nanosponge. And we determined this proportionality factor (k12) based on experimental data on the release of γ-amylase from the DNA nanosponge.
γ-amylase (EC 3.2.1.3) is an exoenzyme attacking both exo-(1-4) and branch-point (1-6)-linkages to produce glucose. The reaction of γ-amylase is inhibited by the accumulation of product glucose at higher substrate starch concentrations. According to the literature [3], we obtained a Kinetic model for inhibition of starch hydrolysis by competitive products by 10U/ml γ-amylase in high concentration soluble solution. We also converted the data from the literature for velocities based on the definition of enzyme activity. （PS. One unit of enzyme activity was defined as the amount of enzyme, which hydrolyzed 1 mg of starch per minute under specified conditions.）
As the concentration of the enzyme increased, the hydrolysis rate increased proportionally in the system containing soluble starches. So we assumed that:
The member responsible for the hardware section fitted the detected glucose concentration to the voltage signal using a modified blood glucose meter.
This gave us the differential equation for the change in glucose concentration on the change in electrical signal.
For the entire assay system, four corresponding kinetic models were constructed in parts. As the wet experimental design was implemented in the same reaction system the Cas-gRNA-activator complex cleaves the DNA nanosponge in trans, releases the amylase hydrolysis substrate to produce glucose and detects the blood glucose concentration by a modified blood glucose meter.
So we combined the CRISPR-Cas system trans-cleavage kinetic model, the γ-amylase hydrolysis starch model and the glucose concentration conversion to electrical signal model. The combined process was simulated with Cas12 as an example, and the curve of glucose concentration in the system with time was obtained. As shown below, the green and orange lines are the simulated data, while the red dots are the experimental data. The R-squared between the simulated output of Cas12a and the experimental data is 0.98872, which is calculated based on the following equation.
And according to the linear equation between the glucose concentration and the electrical signal obtained in the hardware section, the glucose concentration can be converted into an electrical signal to obtain a reading.
In summary, our project achieves a kinetic simulation of the entire assay system. In the RT-RPA part, we took the DNA concentration at 900s as the output by fitting the analysis to the relevant parameters in isothermal amplification, using the RNA of 1E-12M as input. In the Cas section, we took the DNA concentration output from the previous step as input, modelled the enzymatic kinetics of the Cas-gRNA-activator complex trans-cleavage activity and output the concentration of cleaved reporter DNA just reaching plateau. In the γ-amylase section, we converted the cleaved reporter DNA concentration output in the previous step to the γ-amylase release concentration and constructed a kinetic model of γ-amylase hydrolysis with competitive inhibition of the product, taking the glucose concentration just reaching plateau as the output. Finally, in the electrochemical section, we converted the glucose concentration to an electrical signal and obtained the glucometer readings.
We performed molecular dynamics simulation for six kinds of CasΦ mutants collected from literature and wet lab design. The RMSD of different CasΦ mutants were calculated using gmx rms module in the Gromacs package. We analyzed the changes in the skeleton carbon atom of the protein and the central structure of the entire protein.
The average structure after 400 ns of the simulation was extracted to check the RMSD and Rg of the simulated trajectory respectively. As shown in Figure 8, in the two mutant types of CasΦ (Mut-1 and Mut-4), the RMSD for the whole system is far away, with a large shift after 100 ns, after which it keeps converging. While the RMSD of the wild-type CasΦ and other mutants present the stability of about 0.32 nm after 150 ns, with less vibration of the system. This stability result may be related to the conformational changes induced by hydrogen bonds within the protein molecule.
As shown in the Figure 9, after 400 ns of folding, the molecular structure of the seven proteins has partially changed. After 350 ns, they begin to converge and maintain a relatively stable fluctuation range. Among them, the conformational changes of Mut-3 and Mut-4 after a period of equilibrium before 200 ns are small, the convergence degree is high, and the system stability is good. The wild-type protein has a large structural change at 60 and 150 ns, and the internal interaction force is unstable after three-dimensional folding. In the positive control group (Mut-1 and Mut-2), the basic fluctuation range was small after 100 ns, and some small fluctuations occurred at 3.66 nm.
As shown in Figure 10, in the simulated 400 ns, the energy fluctuation amplitude of the whole system is small, the wild-type protein (WT) is the lowest, and the structural stability is the best, followed by the stability of Mut-2 and Mut-3, and the energy of the whole system remains unchanged.
During the simulation process, the number of intramolecular hydrogen bonds in each group of proteins fluctuated to a certain extent, and changed at about 490. Among them, the vibration amplitude of Mut-4 and Mut-6 is smaller than that of other proteins, and the number of hydrogen bonds in the system has been maintained all the time, indicating that the intramolecular interaction is relatively stable, and it is easy to maintain a more rigid protein conformation. (Figure 11)
During the simulation, the displacement range of each atom is small, and it is basically stable. In the Mut-3 and Mut-5 protein groups, the amino acids at 202-218 shifted by 1 nm, indicating that the intramolecular force has changed. The positive control group Mut-1 and Mut-2 basically maintained the displacement below 0.5 nm, and the atomic vibration amplitude was small, indicating that the system was stable. (Figure 12)
[1] Moody, C. et al. “A mathematical model of recombinase polymerase amplificationunder continuously stirred conditions.’’Biochem. Eng. J. 2016, 112, 193– 201
[2] Huyke, Diego A et al. “Enzyme Kinetics and Detector Sensitivity Determine Limits of Detection of Amplification-Free CRISPR-Cas12 and CRISPR-Cas13 Diagnostics.” Analytical chemistry vol. 94,27 (2022): 9826-9834. doi:10.1021/acs.analchem.2c01670
[3] Mian Li et al.“Kinetic enhancement of starch bioconversion in thermoseparating aqueous two-phase reactor systems .”Biochemical Engineering Journal,Volume 11, Issue 1,2002,Pages 25-32,ISSN 1369-703X.
[4] Huyke, Diego A et al. “Enzyme Kinetics and Detector Sensitivity Determine Limits of Detection of Amplification-Free CRISPR-Cas12 and CRISPR-Cas13 Diagnostics.” Analytical chemistry vol. 94,27 (2022): 9826-9834. doi:10.1021/acs.analchem.2c01670
[5] Ramachandran, Ashwin, and Juan G Santiago. “CRISPR Enzyme Kinetics for Molecular Diagnostics.” Analytical chemistry vol. 93,20 (2021): 7456-7464. doi:10.1021/acs.analchem.1c00525
