Math Modeling

Overview:

Modeling in Biosciences is an important tool that can enable us to gain profound insights into our biological systems. It played a pivotal role in guiding and facilitating the development of our detection system, ultimately saving us significant time by steering us away from unproductive design paths. We harnessed basic models to simulate the kinetics of our enzyme cascade in the Cas12-a process and the RCA process. We also combined modeling the process of RCA and Cas12-a, which is a more accurate representation of our actual experiment. To better understand our reaction with the combined mechanism. We also combined reactions of the one-plot assay that can better capture our experiment by linking the Rfree of RCA model with the S in our Cas12a model.

Assumptions:

1. We assume all rate constants remain constant due to a lack of external environmental factors.
2. We assume the total quantity of padlock probes, SplintR ligase, phi29 polymerase, fluorophores, and quenchers is known initially and remains constant over time.
3. We assume that DNA polymerase activity is expected to remain the same regardless of the primer, miRNA, or dNTP concentration.
4. For our calculation of the expansion rate of phi29, we assume that the unit of phi29 is approximately 1 mol in number. 1 U = 1 fmol. We are not sure how much the number is actually due to the unavailability of data from the company. So we can calculate that 0.1 fmole/uL = 1 nM = 10^6 fM. With more data in the future, we can make a better calculation of the k values.
5. For the equilibrium equations of enzyme-substrate binding, we assume that E+SES and that the forward and backward reaction rate is approximately 1. With more data in the future, we can make a better assumption of the k values.

RCA Modeling

One of the important reactions in our model is RCA:, as shown in the figure below.

Based on the characteristics of phi29 polymerase that :

1. elongation will only happen when miRNA primer is present
2. the DNA template will not dissociate
3. phi29 will reach its best efficiency once activated.[1]

We used three kinetics equations to model the cis reaction of E1 and both the cis and trans reaction of E2, as shown below in figure 2.

Together, the three equations map the approximate behavior of the enzymes during the trans-cleavage and cic-cleavage.

To simplify our model due to the limited research available of the RCA mechanism of our specific reaction. We assume that the rate at which MP is produced from M and P and MPSplR produced from MP and SplR, and MPcircphi29 produced from MPcirc and phi29 are of similar value with neglectable differences.

We also calculated the expansion rate of phi29 to determine our k8 value. Based on our WebLab experiment, we added 20ul of phi29. The phi29 we used come from NEB, and from their product description, we found that the concentration is 0.1U/ul and that one unit is defined as the amount of enzyme that will incorporate 0.5 pmol of dNTP into acid insoluble material in 10 minutes at 30°C [2]

20ul x 0.1U/ul=2U

From here we can calculate the rate of phi29 expansion. Based on assumption, one padlock is approximately 60nt in length.

(0.5pmol*1NA/60)/(20*10^-6)=416.65 pM/U*s

2U*416.65 pM/U*s=833.3 pM/s= 8.33x10^6 fM/s

We derived an ODE system based on these equations.

There are many methods and softwares that can be used to solve the system of ODE, such as MatLab, Mathematica, Python Libraries, Scipy etc. These tools provide different numerical methods and solvers for approximating solutions to ODE systems. After careful consideration and extensive research into the available options, we made the decision to utilize the ode15s solver in MATLAB for our study. The ode15s solver is known for its balance between accuracy and computational efficiency. It uses a variable-step, variable-order algorithm that adapts to the characteristics of the ODE system, ensuring that it can handle a wide range of problems efficiently.

Based on our assumption that the Rfree will not degrade, is it reasonable that there is no plateau in our modeled Rfree in RCA process.

Cas12a

One of the important reactions in our model is Cas12a:, as shown in the figure below.

In our model, we used three Michaelis-Menten kinetics equations to model the cis reaction of E1 and both the cis and trans reaction of E2, as shown below.

Based on these three equations, we derived the system of ODE:

dy(1)=k2*y(3)-k1*y(1)*y(2);

dy(2)=-k1*y(1)*y(2)+k2*y(3)-k4*y(5)*y(2)+k5*y(6);

dy(3)=k1*y(1)*y(2)-k2*y(3)-k3*y(3);

dy(4)=k3*y(3)+k6*y(6);

dy(5)=k3*y(3)-k4*y(5)*y(2)+k5*y(6)+k6*y(6)+k8*y(8)-k7*y(5)*y(7)+k9*y(8);

dy(6)=k4*y(5)*y(2)-k5*y(6)-k6*y(6);

dy(7)=-k7*y(5)*y(7)+k8*y(8);

dy(8)=k7*y(5)*y(7)-k8*y(8)-k9*y(8);

dy(9)=k9*y(8);

dy(10)=k9*y(8);

Again, we decided to use the MatLab ode15s solver to solve the system of equations, an example plot of the P3 fluorescent is shown below (Figure 3). In this figure, we can see that there is a clear plateau around 0.2 seconds at a concentration of 5 fM. This is logical because, at lower concentrations, the reactions reach equilibrium relatively quickly. As the reaction progresses, the concentration of the fluorescent product, P3, steadily increases until it levels off at 5 fM.

The plateau indicates that the system has reached a dynamic equilibrium where the rate of formation of P3 equals the rate of its degradation or consumption. In other words, at 0.2 seconds and a concentration of 5 fM, the system has settled into a state where the forward and reverse reactions occur at the same rate. This stability is a fundamental concept in chemical kinetics and is often used to determine reaction mechanisms and rate constants.

RCA-Cas12a Combined Model

To better understand our reaction to the combined mechanism. We were able to produce the combined reaction of the one-plot assay that can better capture our experiment. The Rfree in our RCA model was the same as the S substrate in our cas12a model. Even though in our assumption for the RCA model that there is no upper limit, our combined model will still reach a limit because S2 is the limiting reactant in our case.

Thus, the ODE of our combined model is as shown below:

dy(1)=k2*y(3)-k1*y(1)*y(2);

dy(2)=-k1*y(1)*y(2)+k2*y(3);

dy(3)=k1*y(1)*y(2)-k3*y(3)*y(4)-k2*y(3)+k4*y(5);

dy(4)=k4*y(5)-k3*y(3)*y(4)+k5*y(5);

dy(5)=k4*y(3)*y(4)-k3*y(5)-k5*y(5);

dy(6)=k5*y(5)-k6*y(6)*y(7)+k7*y(8);

dy(7)=-k6*y(6)*y(7)+k7*y(8);

dy(8)=k6*y(6)*y(7)-k7*y(8);

dy(9)=k8*y(8)-k9*y(9)*y(10)+k10*y(11)-k12*y(13)*y(9)+k13*y(14); dy(10)=k10*y(11)-k9*y(9)*y(10);

dy(11)=k9*y(9)*y(10)-k10*y(11)-k11*y(11);

dy(12)=k11*y(11)+k14*y(14);

dy(13)=k11*y(11)-k12*y(13)*y(9)+k13*y(14)+k14*y(14)+k16*y(16)-k15*y(13)*y(15)+k17*y(16);

dy(14)=k12*y(13)*y(9)-k13*y(14)-k14*y(14);

dy(15)=-k15*y(13)*y(15)+k16*y(16);

dy(16)=k15*y(13)*y(15)-k16*y(16)-k17*y(16);

dy(17)=k17*y(16);

dy(18)=k17*y(16);

After using the ode15s solver in MatLab, we were able to produce the reaction curve as shown below. To verify our model with experimental data, we plotted a few of our key reactants with varying miRNA concentrations as experimented in WetLab. All our components reacted in such a way that the rate increases as the miRNA concentration increases.

We also simulated the change in phi29 to see if we can get the similar trend of optimization as we have observed in WetLab results. Based on our wet lab results, 0.2U/ul is the optimal concentration of phi29. Thus, we decided to change the concentration of our phi29 to observe the trend.

However, interestingly, our model did not show the optimal concentration at 0.2U/ul. Instead, it showed that as the concentration of phi29 increases, the reaction increases. This disparity in results between WetLab and Modeling should be further researched.

Future Direction

Our model will provide important insights into confirmation of our experimental results. The reaction rate constants can be better determined through verifying with experimental data. We hope to gather even more data to better calibrate our model. We also hope to quantitatively assess the model sensitivity by measuring how much the model will be affected by changing a small parameter due to potential experimental errors.

In our future modeling study, we would also like to confirm if a change in concentration of SplR, RNP, PLP has a similar effect as we have observed in our wet lab experiment.

There are many ways in which our model affects our project. One most notable contribution of our model is that it allows us to confirm any experimental trends as well as see any discrepancies in our experimental and theoretical results. For example, in our modeling, we have found a trend of increasing concentration of phi29 leading to increasing product, which is different from our experimentally found optimized concentration on 0.2U/ul. While the increasing concentration of miRNA leading to increasing product is confirmed by our experimental results. So in the future, we would like to further experiment with phi29 in a lab setting to find reasons for why the 0.2U/ul concentration showed the greatest product in our current experiment. Once we obtain a more calibrated model, our model can provide a more comprehensive test of initial concentration and a more accurate detection indicator for COPD. By finding the optimal concentration for our detection tool kit, we can improve our detection efficiency and accuracy, while saving necessary resources to bring the price down as much as possible.

Citation

[1]https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2944734/

[2]https://www.neb.com/en/products/m0269-phi29-dna-polymerase#Product%20Information

[3]https://doi.org/10.3390/molecules23123178

[4]https://doi.org/10.1093/nar/gkt1032

[5]https://www.nature.com/articles/s41551-023-01033-1#Abs1