Model

Abstract


Tear fluids contain various types of miRNA. Therefore, to quantify only the selected miRNA as a biomarker for disease, the sequence specificity of the amplification system is crucial. In the Dry Lab, we discovered multiple sequences similar to the miRNA chosen as the target in the POIROT, using the database obtained from a comprehensive search for miRNAs predicted to be present in tear fluids. It was found that identifying these as the target using SDA would be difficult, and it was favored to adopt a highly specific TWJ system for the amplification system.

Introduction


Sequence specificity is one of the most significant factors in the miRNA amplification and quantification system. In amplification systems with low specificity, amplification reactions can proceed due to not only the target miRNA but also miRNAs with similar sequence, resulting in output signals that do not accurately reflect the concentration of the target miRNA in the sample. Since numerous miRNAs are present in samples derived from humans, it is essential to avoid amplification reactions triggered by non-target miRNAs in order to accurately quantify the miRNA selected as a biomarker. In fact, it has been confirmed by the Wet Lab that when SDA without a TWJ complex is used, amplifications do occur even with miRNAs that have only a few differences in the bases compared to the target miRNA. (refer to Results)

On the other hand, the required level of sequence specificity should be determined considering the miRNAs contained in the sample. If miRNAs with similar sequence to the sequence of target miRNA were not contained in the sample, the required level of sequence specificity would be lower. When using complex reactions involving many more types of molecules, such as TWJ, which is known for its high sequence specificity, there is a possibility of drawbacks such as reduced amplification efficiency. Hence, it cannot be conclusively stated that a system with higher specificity is always better. It is important to ensure reliability by designing a system that can accurately quantify the biomarker even under the presence of miRNAs that can potentially exist in the sample and interfere with the system.

With this in mind, the Dry Lab aimed to contribute to the Wet Lab and the project by using data from previous studies to search for similar sequences among the miRNAs that may be present in the sample.

Methods


There are few studies that have comprehensively explored miRNAs in tear fluids, and there was no appropriate data available for searching similar sequences of the miRNAs selected as biomarkers. Therefore, we considered a Meta-analysis, which integrates and analyzes multiple research reports. However, Human Practices to Professor Ochiya, an expert in disease-related miRNAs, has proven that it is quite challenging to extract raw data from different studies, align the data formats, and reanalyze them. It was also revealed that a Meta analysis on miRNAs in tear fluids, an area that is still not researched enough yet, would be unable to ensure reliability. Although research using a Meta-analysis to identify disease-specific biomarker miRNAs is actively conducted, no established method exists yet, and it would be difficult to develop a highly reliable method within the one-year duration of the project. Target miRNAs identified by unreliable methods would be nothing less than abdicating the responsibility as a scientist. It was concluded that it would be more appropriate to refer to the data from a single, reliable paper that has comprehensively explored miRNAs. We referred to data from the study 1 introduced by Professor Ochiya, which conducted a comprehensive exploration of miRNAs in aqueous humor, and we researched the presence of similar sequences in the sample.

The study 1 reported the detection of 1623 types of miRNAs in aqueous humor. However, the data presented in the study does not include the specific sequences of each miRNA, making it impossible to directly search for sequences similar to the target. Therefore, a database containing 1493 miRNAs was created by referencing miRBase Release 22.1 2, which includes numerous reported miRNA sequences, and extracting the human miRNAs with the same names as those included in the data of the study along with their sequences.

General similarity search software, such as BLAST, aims to evaluate the homology between sequences and uses its own scoring system to assess the statistical significance of sequence occurrence frequency 3. While such software is very useful for searching similar sequences, it does not take the interactions between nucleotides during evaluation into consideration. Therefore, it may not accurately assess the concerns of this project regarding the binding of miRNAs similar to the target and templates, and the potential for unintended amplification reactions. Consequently, a program was developed to search for similar sequences using simpler and more intuitive scoring metrics, such as the percentage of matching nucleotides and the number of mismatches. Although these metrics may not accurately evaluate the interactions between miRNAs with similar sequences and templates, they are easier for users to understand. This makes it easier to set search conditions according to the expected performance of the amplification system and helps narrow down similar sequences that should be considered in the selection of target miRNAs and system design, as well as sequences that should be verified through wet experiments or simulations.

In this program, all miRNAs from the database that have a similarity score above a certain threshold when compared to a given miRNA sequence are output. The algorithm is as follows:

  1. Input the following three parameters: "The sequence of the miRNA for which similar sequences are to be searched", "The minimum length of the subsequence to be compared", "The minimum similarity score threshold for the output sequences".
  2. Extract miRNAs from the database.
  3. Compare the lengths of the input sequence and the extracted database sequences, and select an n such that (length of the shorter sequence) ≧ n ≧ (minimum subsequence length)
  4. For each n, extract consecutive subsequences of n nucleotides from both the input sequence and the database sequences, and calculate the similarity (either the percentage of matching nucleotides or the number of mismatches).
  5. Calculate the similarity for all possible n from step 3 and all possible subsequences from step 4, taking the maximum value as the similarity score between the two sequences.
  6. Execute steps 2 to 5 for all sequences in the database, and display all miRNAs with similarity scores above the given threshold.
hoeghoge
Figure 1: Behavior of similar sequence search program.
The program selects sequences from the database and evaluates their similarity to the target sequence by checking for base matches.

In this project, considering the concentration ratios of aqueous humor between glaucoma patients and the control group, hsa-miR-10b-5p, hsa-miR-375, and hsa-miR-30d-5p were selected as target miRNAs. Similarity searches for these three miRNAs were conducted using the aforementioned methods. Taking into account the results of specificity evaluation from wet experiments using the 22 nt miRNA hsa-let-7b and the let-7 family, the minimum subsequence length was set to 20, and similar sequences with a sequence match ratio of 80% or higher to the target miRNAs were extracted.

Results


hsa-miR-10b-5p

miRNAs that are similar in sequence to hsa-miR-10b-5p include hsa-miR-10a-5p, hsa-miR-100-5p, hsa-miR-99a-5p, and hsa-miR-99b-5p. miR-10a-5p consists of 23 nt, differing from miR-10b-5p at only the 12th nucleotide. miR-100-5p, miR-99a-5p, and miR-99b-5p consist of 22 nucleotides, with 4, 5, and 6 mismatches, respectively, when aligned with the 3' end of miR-10b-5p.

hsa-miR-375

As a result of the similarity sequence search, no miRNAs with sequences similar to hsa-miR-375 were found.

hsa-miR-30d-5p

miRNAs with sequences similar to hsa-miR-30d-5p were identified as hsa-miR-30a-5p and hsa-miR-30e-5p. All of these miRNAs are 22 nucleotides long, with hsa-miR-30a-5p having a mismatch of 1 nucleotide and hsa-miR-30e-5p having a mismatch of 2 nucleotides.

hoeghoge
Figure 2: Similar miRNA sequences found in similar sequence search of (a) hsa-miR-10b-5p and (b) hsa-miR-30d-5p.
Light green bases in similar sequences represent mismatches with the target sequence.

Conclusion and Future Prospects


As a result of the biomarker similarity sequence search for miRNAs expected to be present in tear fluid, similar sequences were found for miR-10b-5p and miR-30d-5p. Among these, there are miRNAs that differ by only 1 or 2 nucleotides from the biomarkers, making it challenging to distinguish them using simple SDA based on wet results. These similar miRNAs found in the aqueous humor are also thought to be present in tear fluid, potentially leading to false detections during quantitative analyses using tears. Therefore, it is supported that a more sequence-specific TWJ, rather than a simple SDA, should be employed in the amplification system, leading to improvements in the project.

The program created this time can search for similar sequences of miRNAs contained in various samples other than aqueous humor by changing the database being targeted. POIROT can be used for the quantification of any miRNA by modifying the template sequence. Therefore, when applied to samples other than tear fluid, it is expected to contribute to the appropriate selection of biomarkers, ensuring the specificity and quantitativeness of the system by searching for similar sequences of target miRNAs within a database that includes miRNAs present in the sample.

In the current search for similar sequences, a similarity evaluation metric focusing solely on the matches and mismatches of nucleotide sequences was used. However, it is also possible to evaluate the binding energy between miRNAs with similar sequences to the template and target, as well as assess the impact of these similar sequences on the amplification system through simulations using an ODE model. Since these evaluations take longer to compute compared to the similarity evaluation methods used in this analysis, extracting similar sequences from numerous miRNAs using this program enables efficient assessment.

The program and database used for the sequence search in this study have been incorporated into the software and are publicly available on GitLab. For more details, please visit: Software

References


  1. Tanaka, Y., Tsuda, S., Kunikata, H. et al. (2014). Aqueous Humor of Glaucoma Patients Assessed with a Microarray System. Sci Rep 4, 5089. https://doi.org/10.1038/srep05089

  2. Faculty of Biology, Medicine and Health, The University of Manchester. miRBase. https://www.mirbase.org

  3. The National Center for Biotechnology Information. The Statistics of Sequence Similarity Scores. https://www.ncbi.nlm.nih.gov/BLAST/tutorial/Altschul-1.html

Abstract


Our goal is to develop a home testing device for glaucoma. To this end, we explored the use of strand displacement amplification (SDA) methods that do not require complex temperature changes. In Dry lab, we constructed an ordinary differential equation (ODE) model to select appropriate amplification reactions. First, based on the experimental data, we fitted the ODE model parameters that were unknown and experimentally challenging to determine to replicate the experimental results in Wet lab. The ODE model revealed that three-way junction (TWJ-SDA) alone was insufficient to activate Cas due to low amplification efficiency. Based on previous researches, then, we proposed the alternative amplification systems, including TWJ-SDA > Exponential amplification reaction (EXPAR), TWJ-SDA > 2step-SDA and TWJ-SDA > 3step-SDA, and compared amplification efficiency, the positive-to-negative signal ratio (P/N), and robustness to enzyme activity and time using the ODE model. As a result, we proposed using TWJ-SDA > 3step-SDA as the best amplification system and suggested the optimal initial template concentration.

Fitting the parameters

Purpose

In constructing the ODE model, it is necessary to define the parameters used within the model. However, in this amplification system, the reaction rate constants for nucleic acid binding and nicking, and the Michaelis constant were unknown and difficult to determine experimentally. The goal of this section is to appropriately estimate these parameters based on wet lab experimental results and construct an ODE model that replicates a portion of these results.

Methods

In constructing the ODE model, the parameters like these were generally determined through fitting. We developed several ODE models throughout the project. As several ODE modellings of EXPAR had been constructed, the parameters on EXPAR were fitted first and the fitted parameters were used in the ODE models of other amplification systems.

In EXPAR, the process begins with the target binding to the template, followed by a polymerase extension reaction that produces double strands (Fig. 1). Next, nickase recognizes the nicking recognition sequence within the double strands, leading to nicking. The polymerase then initiates a new extension reaction from the nicking site, displacing the downstream single-strand DNA (ssDNA) encountered during synthesis and producing the full double-strands DNA (dsDNA). Since the displaced ssDNA has the same sequence as the target, it binds to the template again, initiating further extension reactions and achieving overall exponential amplification.

hogehoge
Fig. 1. The reaction mechanism of EXPAR.

The ODE model for EXPAR in this project was constructed with the references to multiple prior studies that made ODE models of EXPAR 1, 2 and the reaction mechanism of ultrasensitive DNA amplification reactions 3, which is similar to that of EXPAR (Fig. 2). As X could bind to the 5' end of T, ODE models were built in these studies based on different scenarios, such as X binding to both the 3' and 5' ends of T or only to the 3' end. These ODE models were constructed for each case by Dry lab. Still, since no significant difference was observed in the amplification behaviour of the reaction products, we opted to construct the simpler ODE model, assuming X binds only to the 3' end of T to reduce the number of parameters.

hogehoge
Fig. 2. The model of EXPAR.
X: target. T: template. XT: the complex of X and T. W: the complete dsDNA. Z: the nicked dsDNA.

This ODE model is formulated as shown in Eq. 1. The way to determine the reaction rates in this ODE model follows prior research 3, with the assumption that binding and dissociation reactions follow the law of mass action and enzyme reactions follow the Michaelis-Menten equation. For more assumptions in the construct of ODE model, see Assumption The experiments of EXPAR were conducted at 50 ℃, with an initial T concentration of 4.40e-2 µM and an initial X concentration ranging from 0 M (NC) to 1 pM-1 µM SYBR Green Ⅰ was used to detect fluorescence, which indicates the presence of dsDNA. Bst 2.0 polymerase, optimal at 65 ℃, was employed, and Nt.BstNBI nickase, optimal at 55 ℃, was used. Their initial concentrations were 1.50e-2 U/µL and 9.38e-2 U/µL, respectively. The parameters used in the ODE model are listed below in Tab. 1. Notably, \(k_{cat}P_{0}\) is expressed as such because \(k_{cat}\) and \(P_{0}\) appear only as a product and not individually. Given that 1 unit of Bst 2.0 polymerase is defined as the amount of the enzyme required to incorporate 25 nmol of dNTP into acid-insoluble material in 30 minutes at 65℃ 4, U is translated to mol as \(1U=25 nmol/1800 s/k_{cat}\), eliminating the need to measure the exact value of \(k_{cat}\)​. Since the fitting experiments were conducted at 50℃, where polymerase activity is reduced to 30-50% 5, we assumed 40% activity for the calculations. For the initial nickase concentration, we referred to prior studies 3 modelling the reaction in which the same enzyme was used and 0.2 U/µL corresponded to 2.6e-2 µM. Nt.BstNBI Nickase retains 100% activity at 50℃. Additionally 6, the experimental data suggested its activity remains constant during the reaction, though a prior research considered the heat inactivation of nickase. The association constants were calculated using NUPACK. For more detailed experimental methods, see Experiments_EXPAR

\begin{align} \dfrac{d[X]}{dt}&=-k_{1}[X][T]+\dfrac{k_{1}}{a}[XT]+\dfrac{k_{cat}}{17}\dfrac{P_{0}}{mc}[Y]\notag\\ \dfrac{d[T]}{dt}&=-k_{1}[X][T]+\dfrac{k_{1}}{a}[XT]\notag\\ \dfrac{d[XT]}{dt}&=k_{1}[X][T]-\dfrac{k_{1}}{a}[XT]-\dfrac{k_{cat}}{27}\dfrac{P_{0}}{mc}[XT]\notag\\ \dfrac{d[W]}{dt}&=\dfrac{k_{cat}}{27}\dfrac{P_{0}}{mc}[XT]-k_{2}[N][W]+\dfrac{k_{cat}}{17}\dfrac{P_{0}}{mc}[Y]\notag\\ \dfrac{d[Y]}{dt}&=k_{2}[N][W]-\dfrac{k_{cat}}{17}\dfrac{P_{0}}{mc}[Y]\notag\\ \end{align}

Eq. 1 The equations of EXPAR
.

うっひょー

\[ \begin{array}{|c|c|c|c|c|} \hline \textbf{Parameter} & \textbf{Definition} & \textbf{Value} & \textbf{Unit} & \textbf{Source} \\ \hline k_{1} & \text{Reaction rate constant} & 1.14 \times 10^{1} & \mu M^{-1}s^{-1} & \text{Fitting} \\ k_{2} & \text{Nicking rate constant} & 1.00 \times 10^{1} & \mu M^{-1}s^{-1} & \text{Fitting} \\ k_{cat} & \text{Turnover number} & k_{cat}P_{0} = 8.33 \times 10^{-2} & s^{-1} & \text{NEB}^{4,5} \\ P_{0} & \text{Polymerase initial concentration} & k_{cat}P_{0} = 8.33 \times 10^{-2} & \mu M & \text{NEB}^{4,5} \\ N_{0} & \text{Nickase initial concentration} & 1.22 \times 10^{-2} & \mu M & \text{Özay, B. et. al.}^{3}, \text{NEB}^{6} \\ a & \text{Association constant (X, T)} & 5.43 \times 10^{2} & \mu M^{-1} & \text{NUPACK} \\ m & \text{Michaelis–Menten constant} & 1.08 \times 10^{-1} & \mu M & \text{Fitting} \\ c & 1 + \frac{[XT] + [Z]}{m} & - & - & \text{Özay, B. et. al.}^{3} \\ \hline \end{array} \]

The mean fluorescence of the three experiments for each target concentration, excluding the background, was as follows (Fig. 3).

hogehoge
Fig. 3. The fluorescence of EXPAR (Wet).

The fluorescence was predicted by the simulation as follows (Fig. 4). Since the fluorescence is measured by using SYBR Green Ⅰ, which basically detects double-stranded DNA, it was assumed that the total concentration of W and Z was proportional to the fluorescence intensity, as well as the prior study 1 that also performed fitting of an ODE model of EXPAR. When the concentration graph reaches a plateau, the total concentration of W and Z was approximately equal to the initial template concentration of 4.4e-2 μM as there was almost no T or XT in the simulation. Therefore, this value 4.4e-2 was assumed to correspond to the plateau fluorescence value of 3.00e5 in Fig. 3. Both the fluorescence graph from the wet experiment and the concentration graph from the dry simulation reach a plateau. After reaching the plateau, however, the concentration remains stable in the simulation, while the experimental fluorescence continues to increase. According to advice from Professor Komiya in Human Practices, this might be due to SYBR Green Ⅰ emitting slight fluorescence in response to increasing amounts of X, which is much more abundant than dsDNA.

The parameters \(k_{1},k_{2},m\) was fitted using a genetic algorithm with a population size of 100 and 50 generations, and the used data was located between the dotted lines (0~2.80e5) in Fig. 4. In the simulation, the initial fluorescence was set to 0, but in the wet experiment, when the initial X concentration was high, the amplification began immediately, and before the solution was placed in the measurement device, the graph of the fluorescence started to rise. Hence, in Dry lab simulation, it was assumed that the reaction started, and that 1 minute later the fluorescence measurement began. The amplification of the negative control (NC) and of the high concentrations like 100 nM and 1 μM, where the fluorescences reached at 2.80e5 within the first minute was excluded from the data used in the fitting. The amplification of NC is thought to occur because of ab initio synthesis, which refers free dNTPs bind to the template without primers, triggering the extension reaction 7.

hogehoge
Fig. 4. The fluorescence of EXPAR (Dry).

A genetic algorithm is a method that explores solutions by representing data (a set of proposed solutions) as "individuals" with genes and repeatedly performing operations such as selection, crossover, and mutation, selecting individuals with high fitness. In this project, we used a genetic algorithm to fit the reaction rate constants of nucleic acid binding and nicking, and the Michaelis constant.
First, in each generation, 100 individuals were determined with genes representing the three constants, and the fitness was defined by how closely the simulated results matched the experimental data based on the mean squared error (MSE). Individuals with smaller MSEs were considered fitter. The algorithm was designed so that, for each generation, individuals were defined to undergo crossover with probability of 50 %, mutation with that of 20%, and copied with that of 30 %. For crossover, in which two individuals need to be selected, mutation and copying, in which one individual needs to be selected, a selection process known as "tournament selection" was used. The tournament selection involves randomly picking several individuals and selecting the best one among them. In this fitting process, two individuals were randomly chosen, and the individual the fitness of which is the higher was selected. The tournament selection was occured twice for crossover, and once for mutation and copying. By limiting the number of participants in each tournament, the randomness of this algorithm increases, promoting the diversity in the selected individuals.

For crossover, the "Blend Crossover" method was applied. Blend Crossover generates the new genes of the child individual by mixing the genes of two parent individuals. The gene of the child individual is randomly selected from an extended range of the parent individuals' genes, controlled by a parameter called \(\alpha\). Let's assume the two parent individuals \(P_{1},P_{2}\) have genes \(x_{1},x_{2}\), respectively. When performing Blend Crossover, the gene c of the child individual is randomly chosen from the following range.

\begin{align} c \in \left[ \min(x_{\text{min}}) - \alpha \cdot (\max(x_{\text{max}}) - \min(x_{\text{min}})), \max(x_{\text{max}}) + \alpha \cdot (\max(x_{\text{max}}) - \min(x_{\text{min}})) \right]\notag\\ \end{align}

In these equation, \(x_{\text{min}}, x_{\text{max}}\) were defined as follows.

\begin{align} x_{\text{min}} = \min(x_1, x_2)\notag\\ x_{\text{max}} = \max(x_1, x_2)\notag\\ \end{align}

\(\alpha\) is the parameter controlling the degree of expansion in the range for the child individual's genes and typically takes values in the range \(0 \leq \alpha \leq 1\). In this way, Blend Crossover allows the child individual's genes to explore beyond the range of the parent genes, expanding the search space and reducing the risk of getting stuck in a local optimum.

Next, in the mutation process for this fitting, each gene of an individual had a 20% probability of being modified by adding a random value drawn from a Gaussian distribution. The Gaussian distribution has a mean of 0 and a standard deviation of 0.01. The child individual's gene is calculated as follows

\begin{align} x_{\text{new}} = x + N(0, 0.01)\notag\\ \end{align}

The small standard deviation allows for fine-tuning of the genes, and by gradually changing the genes, the search range is expanded, preventing the algorithm from settling into a local optimum. Mutation is essential for exploring solutions that cannot be found by crossover alone.

Additionally, copying occurs with a 30% probability, and these three processes were repeated until the number of the child generation reaches 100. After that, the child generation became the current generation, and the same process was repeated for 50 generations. Finally, the individual with the highest fitness was selected and its genes were used as the fitted constants.

Result

The parameters obtained through fitting were listed in Tab. 1, and the fluorescence of the simulation results compared with that of the experimental results was shown below (Fig. 5).

hogehoge
Fig. 5. The fluorescence of EXPAR (Fitting).

Conclusion

The purpose of this section was to determine the unknown parameters in the ODE model of the amplification system in this project, specifically the reaction rate constants for nucleic acid binding, the nicking reaction, and the Michaelis constant, which are tough to obtain experimentally. By employing a genetic algorithm, we were able to fit the experimental results of EXPAR with the simulation results from the ODE model to evaluate these unknown constants. The determined values are summarised in Tab. 1.

The amplification rate of TWJ-SDA


Purpose

In the previous section discussing specificity, the high specificity of TWJ-SDA (Fig. 6) was demonstrated. For more precise information, see Model. However, the amplification efficiency of TWJ-SDA is limited by the nature of the reaction, where the target binds to the template and proceeds through strand displacement, resulting in a maximum yield of amplification products determined by the target concentration. The purpose of this section was to evaluate whether TWJ-SDA alone had sufficient effeciency as the amplification system in this project by using simulations based on the ODE model of TWJ-SDA. While the presence of amplification can be examined by observing the fluorescence, the actual concentration of the amplification products cannot be measured, making this modelling essential.

hogehoge
Fig. 6. The reaction mechanism of TWJ-SDA.

Method

We constructed the ODE model for TWJ-SDA based on the reaction mechanism of TWJ-SDA (Fig. 7). In this model, the binding of the target to the template and that of primer to the template were not considered, as they did not contribute to the overall amplification efficiency and the robustness of the amplification system, which is a key criterion for evaluating the amplification system. By excluding these complexes, it was also aimed to reduce the number of parameters in this ODE model.

hogehoge
Fig. 7. TThe model of TWJ-SDA.
C: the complex of target and primer. T: template. CT: the complex of C and T. W: the complete dsDNA. Z: the nicked dsDNA. X: the product.

\[ \begin{array}{|c|c|} \hline \textbf{Complex} & \textbf{Concentration (nM)} \\ \hline \text{CT} & 8.34 \times 10^{2} \\ \text{T} & 1.66 \times 10^{2} \\ \text{C} & 1.65 \times 10^{2} \\ X_{1}T & 4.10 \times 10^{-2} \\ \text{PT} & 2.83 \times 10^{-2} \\ \hline \end{array} \]

Tab. 2. The concentrations of TWJ-SDA complexes
.

Based on the above model diagram (Fig. 7), the following ODE was formulated (Eq. 2). Considering practical device implementations, the simulation was conducted setting the temperature, which was required to determine enzyme activity, to 37 ℃. Since the amplification efficiency did not change significantly when there was an ample amount of template, excess initial template concentration of 1 M and the initial target concentration of 100 fM were assumed in the simulation. The enzyme concentrations were optimised based on Wet experiments. For Polymerase, we used Bst. Large Fragment, which has an optimal temperature of 65 ℃, with an initial concentration of 6.25e-2 U/µL, and for Nickase, we used Nb.BbvCI, which has an optimal temperature of 37 ℃, with an initial concentration of 1.20e-1 U/µL. As for the sequences, the target was designated as a biomarker for glaucoma, with the template and primer being complementary to it.

The constants used in the ODE were referenced from fitting, NEB, and NUPACK (Tab. 3). The definition of 1 unit of the Polymerase, Bst Large Fragment, is the amount of enzyme required to incorporate 10 nmol of dNTP into acid-insoluble material in 30 minutes at 65 ℃ 9, expressed as \(1U=10 nmol/1800s/k_{cat}\). The Polymerase activity at 37 ℃ is calculated to be 12.5%, as it ranges from 10% to 15% 10. For the initial concentration of Nickase, we performed unit conversion from U to mol in the same manner as done in EXPAR, assuming the Nickase, Nt. BbvCI, operated at 100% activity at the optimal temperature of 37 ℃ 6.

\begin{align} \frac{d[C]}{dt}=-k_{1}[C][T]+\frac{k_{1}}{a}[CT]\notag\\ \frac{d[T]}{dt}=-k_{1}[C][T]+\frac{k_{1}}{a}[CT]\notag\\ \frac{d[CT]}{dt}=k_{1}[C][T]-frac{k_{1}}{a}[CT]-\frac{k_{cat}}{26}\frac{P_{0}}{mc}[CT]\notag\\ \frac{d[W]}{dt}=\frac{k_{cat}}{26}\frac{P_{0}}{mc}[CT]-k_{2}[N][W]+\frac{k_{cat}}{18}\frac{P_{0}}{mc}[Z]\notag\\ \frac{d[Z]}{dt}=k_{2}[N][W]-\frac{k_{cat}}{18}\frac{P_{0}}{mc}[Z]\notag\\ \frac{d[X]}{dt}=\frac{k_{cat}}{18}\frac{P_{0}}{mc}[Z]\notag\\

Eq. 1 The equations of EXPAR
.

\[ \begin{array}{|c|c|c|c|c|} \hline \textbf{Parameter} & \textbf{Definition} & \textbf{Value} & \textbf{Unit} & \textbf{Source} \\ \hline k_{1} & \text{Reaction rate constant} & 1.14 \times 10^{1} & \mu M^{-1}s^{-1} & \text{Fitting} \\ k_{2} & \text{Nicking rate constant} & 1.00 \times 10^{1} & \mu M^{-1}s^{-1} & \text{Fitting} \\ k_{cat} & \text{Turnover number} & k_{cat}P_{0}= 4.34 \times 10^{-2} & s^{-1} & \text{NEB}^{9,10} \\ P_{0} & \text{Polymerase initial concentration} & k_{cat}P_{0}= 4.34 \times 10^{-2} & \mu M & \text{NEB}^{9,10} \\ N_{0} & \text{Nickase initial concentration} & 1.56 \times 10^{-2} & \mu M & \text{Özay, B. et. al.}^{3}, \text{NEB}^{6} \\ a & \text{Association constant (C, T)} & 3.03 \times 10^{1} & \mu M^{-1} & \text{NUPACK} \\ m & \text{Michaelis–Menten constant} & 1.08 \times 10^{-1} & \mu M & \text{Fitting} \\ c & 1 + \frac{[\text{CT}] + [\text{Z}]}{m} & - & - & \text{Özay, B. et. al.}^{3} \\ \hline \end{array} \]

Results

hogehoge
Fig. 8. TWJ-SDA (Dry).

When considering the reaction time upon device implementation, it was estimated to be around 40 minutes. However, the simulation results (Fig. 8) indicated that after 40 minutes, the concentration of the amplification product reached approximately only 3 pM. This concentration is far from the required 2.5 nM concentration of the PAM sequence needed for Cas activation, which was shown by the Wet result, making it clear that achieving this efficiency with TWJ-SDA alone is virtually impossible. For more information about the result of Cas activation, see results

Conclusion

From the above, it is evident that in our project to develop a glaucoma testing device, the amplification efficiency of TWJ-Amp is insufficient. The ODE model suggests that it is necessary to combine TWJ-Amp with a more efficient amplification reaction to achieve the desired results.

Optimise Amplification


Purpose

The TWJ-SDA model indicated that it was insufficient for the amplification of PAM-containing dsDNA to the approximately required 2.5 nM for Cas recognition from the result of experiment. For more information about the result of Cas activation, see results
Therefore, it was necessary to use the amplified product of TWJ-SDA as a target and conduct other reactions with high amplification efficiency. Based on a review of prior studies, linking EXPAR and Multi step-SDA 11 were proposed, suggesting TWJ-SDA > EXPAR, TWJ-SDA > 2step-SDA and TWJ-SDA > 3step-SDA as the overall amplification system (Fig. 9, 10, 11). However, experimentally testing these three amplification systems was challenging due to cost constraints (money, effort, and time). Consequently, the objective of this section was to propose the optimal amplification system and initial concentrations using simulations based on the ODE model.

hogehoge
Fig. 10. TThe reaction mechanism of TWJ-SDA > 2step-SDA.
hogehoge
Fig. 11. The reaction mechanism of TWJ-SDA > 3step-SDA.

Methods

We constructed the ODE models of the three amplification systems (Fig. 12, 13, 14; Eq. 3, 4, 5; Tab. 4, 5, 6). In the previous ODE models for EXPAR and TWJ-SDA, the binding of amplification products to the 5' end of the template was not considered. However, since when combining multiple amplification reactions, this binding can inhibit subsequent amplification reactions, the ODE models that account for this factor were constructed. Since there was no existing ODE model for Multi step-SDA in prior studies, the ODE models were constructed based on the reaction mechanisms (Fig. 9, 10, 11) as those of EXPAR and TWJ-SDA. Assuming practical device implementations, the simulation was conducted setting the temperature, which is required to determine enzyme activity, to 37 ℃, and using enzyme concentrations optimised through actual Wet experiments. For Polymerase, Bst. Large Fragment, which has an optimal temperature of 65 ℃, was used with an initial concentration of 8.00e-2 U/µL, and for Nickase, Nb.BbvCI, which has an optimal temperature of 37 ℃, was used with an initial concentration of 2.00e-1 U/µL. As for sequences, they could be freely changed for the amplification products, including TWJ-SDA, and were discussed in other sections. Thus, they were omitted from these ODE models discussion as they were irrelevant to the overall amplification efficiency and the robustness of the amplification system. The constants used in the ODE were determined in the same way as for TWJ-SDA. For binding constants, we used the same constants if the base pairs near the binding sites remain unchanged (Tab. 4, 5, 6).

hogehoge
Fig. 12. The model of TWJ-SDA > EXPAR.

\begin{align} \frac{d[C]}{dt}&=-k_{1}[C][T_{1}]+\frac{k_{1}}{a_{1}}[CT_{1}]-k_{1}[C][X_{2}T_{1}]+\frac{k_{1}}{a_{1}}[X_{2}CT_{1}]\notag\\ \frac{d[T_{1}]}{dt}&=-k_{1}[C][T_{1}]+\frac{k_{1}}{a_{1}}[CT_{1}]-k_{1}[X_{2}][T_{1}]+\frac{k_{1}}{a_2}[X_{2}T_{1}]\notag\\ \frac{d[CT_{1}]}{dt}&=k_{1}[C][T_{1}]-frac{k_{1}}{a_{1}}[CT_{1}]-k_{1}[x_{2}][CT_{1}]+\frac{k_{1}}{a_{2}}[X_{2}CT_{1}]-\frac{k_{cat}}{28}\frac{P_{0}}{mc}[CT_{1}]\notag\\ \frac{d[W_{1}]}{dt}&=\frac{k_{cat}}{28}\frac{P_{0}}{mc}[CT_{1}]-k_{2}[N][W_{1}]+\frac{k_{cat}}{22}\frac{P_{0}}{mc}[Z_{1}]\notag\\ \frac{d[Z_{1}]}{dt}&=k_{2}[N][W_{1}]-\frac{k_{cat}}{22}\frac{P_{0}}{mc}[Z_{1}]+\frac{k_{cat}}{6}\frac{P_{0}}{mc}[X_{2}CT_{1}]\notag\\ \frac{d[X_{2}CT_{1}}]&=k_{1}[CT_{1}][Z_{2}]-\frac{k_{1}}{a_{2}}-\frac{k_{cat}}{6}\frac{P_{0}}{mc}[X_{2}CT_{1}]+k_{1}[C][X_{2}T_{1}]-\frac{k_{1}}{a_{1}}[X_{2}CT_{1}]\notag\\ \frac{d[X_{2}T_{1}}]&=k_{1}[X_{2}][T_{1}]-\frac{k_{1}}{a_{2}}[X_{2}T_{1}]-k_{1}[C][X_{2}T_{1}]+\frac{k_{1}}{a_{1}}[X_{2}CT_{1}]\notag\\ \frac{d[X_{2}]}{dt}&=-k_{1}[X_{2}][CT_{1}]+\frac{k_{1}}{a_{2}}[X_{2}CT_{1}]-k_{1}[X_{2}][T_{1}]+\frac{k_{1}}{a_{2}}[X_{2}T_{1}]+\frac{k_{cat}}{22}\frac{P_{0}}{mc}[Z_{1}]+\frac{k_{cat}}{22}\frac{P_{0}}{mc}[Z_{2}]-k_{1}[X_{2}][T_{3}]+\frac{k_{1}}{a_{2}}[X_{2}T_{3}]\notag\\ \frac{d[T_{2}]}{dt}&=k_{1}[X_{2}][T_{2}]-\frac{k_{1}}{a_{2}}[X_{2}T_{2}^{3}]+k_{1}[X_{2}][T_{2}]-\frac{k_{1}}{a_{2}}[X_{2}T_{2}^{5}]\notag\\ \frac{d[X_{2}T_{2}^{5}]}{dt}&=k_{1}[X_{2}][T_{2}]-\frac{k_{1}}{a_{2}}[X_{2}T_{2}^{5}]-k_{1}[X_{2}][X_{2}T_{2}^{5}]+\frac{k_{1}}{a_{2}}[X_{2}T_{2}]-\frac{k_{cat}}{22}\frac{P_{0}}{mc}[X_{2}T_{2}^{5}]\notag\\ \frac{d[X_{2}T_{2}^{3}]}{dt}&=k_{1}[X_{2}][T_{2}]-\frac{k_{1}}{a_{2}}[X_{2}T_{2}^{3}]-k_{1}[X_{2}][X_{2}T_{2}^{3}]+\frac{k_{1}}{a_{2}}[X_{2}T_{2}]\notag\\ \frac{d[X_{2}T_{2}]}{dt}&=k_{1}[X_{2}][X_{2}T_{2}^{5}]-\frac{k_{1}}{a_{2}}[X_{2}T_{2}]+k_{1}[X_{2}][X_{2}T_{2}^{3}]-\frac{k_{1}}{a_{2}}[X_{2}T_{2}]-\frac{k_{cat}}{4}\frac{P_{0}}{mc}[X_{2}T_{2}]\notag\\ \frac{d[Z_{2}]}{dt}&=k_{2}[N][W_{2}]-\frac{k_{cat}}{22}\frac{P_{0}}{mc}[Z_{2}]+\frac{k_{cat}}{4}\frac{P_{0}}{mc}[X_{2}T_{2}]\notag\\ \frac{d[W_{2}]}{dt}&=\frac{k_{cat}}{26}\frac{P_{0}}{mc}[X_{2}T_{2}^{5}]-k_{2}[N][W_{2}]+\frac{k_{cat}}{22}\frac{P_{0}}{mc}[Z_{2}]\notag\\ \frac{d[T_{3}]}{dt}&=-k_{1}[X_{2}][T_{3}]+\frac{k_{1}}{a_{2}}[X_{2}T_{3}]\notag\\ \frac{d[X_{2}T_{3}]}{dt}&=k_{1}[X_{2}][T_{3}]-\frac{k_{1}}{a_{2}}[X_{2}T_{3}]-\frac{k_{cat}}{26}\frac{P_{0}}{mc}[X_{2}T_{3}]\notag\\ \frac{d[W_{4}]}{dt}&=\frac{k_{cat}}{26}\frac{P_{0}}{mc}[X_{2}T_{3}]\notag\\ \end{align}

Eq. 3 The equations of TWJ-SDA > EXPAR
.

\[\begin{table}[h!] \centering \begin{tabular}{|l|l|l|l|l|} \hline \textbf{Parameter} & \textbf{Definition} & \textbf{Value} & \textbf{Unit} & \textbf{Source} \\ \hline $k_{1}$ & Reaction rate constant & 1.14e1 & $\mu M^{-1}s^{-1}$ & Fitting \\ \hline $k_{2}$ & Nicking rate constant & 1.00e1 & $\mu M^{-1}s^{-1}$ & Fitting \\ \hline $k_{cat}$ & Turnover number & $k_{cat}P_{0} = 5.56e-2$ & $s^{-1}$ & NEB \cite{amp9} \cite{amp10} \\ \hline $P_{0}$ & Polymerase initial concentration & $k_{cat}P_{0} = 5.56e-2$ & $\mu M$ & NEB \cite{amp9} \cite{amp10} \\ \hline $N_{0}$ & Nickase initial concentration & 2.60e-2 & $\mu M$ & Özay et al. \cite{amp3}, NEB \cite{amp6} \\ \hline $a_{1}$ & Association constant (C, $T_{1}$), (C,$X_{2}T_{1}$) & 2.84e1 & $\mu M^{-1}$ & NUPACK \\ \hline $a_{2}$ & Association constant ($CT_{1}$, $X_{2}$), ($X_{1}$,$T_{1}$), ($X_{2}$,$T_{2}$), ($X_{2}$,$T_{3}$), ($X_{2}$,$X_{2}T_{2}$) & 1.77e7 & $\mu M^{-1}$ & NUPACK \\ \hline m & Michaelis–Menten constant & 1.08e-1 & $\mu M^{-1}$ & Fitting \\ \hline c & $1 + \left( \left[ CT_{1} \right] + \left[ Z_{1} \right] + \left[ X_{2}CT_{1} \right] + \left[ X_{2}T_{2} \right] + \left[ Z_{2} \right] + \left[ X_{3}X_{2}T_{2} \right] + \left[ X_{2}T_{3} \right] + \left[ X_{2}T_{2} \right] \right) / m$ & - & - & \cite{amp3} \\ \hline \end{tabular} \caption{Parameter definitions and values.} \end{table} \]

hogehoge
Fig. 13. The model of TWJ-SDA > 2step-SDA.

\begin{align}\frac{d[C]}{dt} $= - k_{1} [T_{1}] [C] + \frac{k_{1}}{a_{1}} [CT_{1}] - k_{1} [C] [X_{2}T_{1}] +\frac{k_{1}}{a_{4}} [X_{2}CT_{1}]\notag\\ \frac{d[T_{1}]}{dt} $= - k_{1} [T_{1}] [C] + \frac{k_{1}}{a_{1}} [CT_{1}] - k_{1}[T_{1}][X_{2}]+\frac{k_{1}}{a_{3}} [X_{2}T_{1}]\notag\\ \frac{d[CT_{1}]}{dt} $= k_{1} [C] [T_{1}] - \frac{k_{1}}{a_{1}} [CT_{1}] - k_{1} [CT_{1}] [X_{2}] + \frac{k_{1}}{a_{2}} [X_{2}CT_{1}] - \frac{k_{cat}}{28 m c} [CT_{1}]\notag\\ \frac{d[Z_{1}]}{dt} $= k_{2N} [W_{1}] - \frac{k_{cat}}{22 m c} [Z_{1}] + \frac{k_{cat}}{6 m c} [X_{2}CT_{1}]\notag\\ \frac{d[W_{1}]}{dt} $= \frac{k_{cat}}{28 m c} [CT_{1}] - k_{2N} [W_{1}] + \frac{k_{cat}}{22 m c} [Z_{1}]\notag\\ \frac{d[X_{2}CT_{1}]}{dt} $= k_{1} [CT_{1}] [X_{2}] - \frac{k_{1}}{a_{2}} [X_{2}CT_{1}] - \frac{k_{cat}}{6 m c} [X_{2}CT_{1}] + k_{1} [C] [X_{2}T_{1}] - \frac{k_{1}}{a_{4}} [X_{2}CT_{1}]\notag\\ \frac{d[X_{2}T_{1}]}{dt} $= - \frac{k_{1}}{a_{3}} [X_{2}T_{1}] + k_{1} [T_{1}] [X_{2}] - k_{1} [C] [X_{2}T_{1}] + \frac{k_{1}}{a_{4}} [X_{2}CT_{1}]\notag\\ \frac{d[X_{2}]}{dt} $= \frac{k_{cat}}{22 m c} [Z_{1}] - k_{1} [T_{1}] [X_{2}] + \frac{k_{1}}{a_{3}} [X_{2}T_{1}] - k_{1} [CT_{1}] [X_{2}] + \frac{k_{1}}{a_{2}} [X_{2}CT_{1}] + \frac{k_{1}}{a_{5}} [X_{2}T_{2}] - k_{1} [X_{2}] [T_{2}] - k_{1} [X_{2}] [X_{3}T_{2}] + \frac{k_{1}}{a_{8}} [X_{3}X_{2}T_{2}]\notag\\ \frac{d[T_{2}]}{dt} $= - k_{1} [T_{2}] [X_{2}] + \frac{k_{1}}{a_{5}} [X_{2}T_{2}] - k_{1} [T_{2}] [X_{3}] + \frac{k_{1}}{a_{7}} [X_{3}T_{2}]\notag\\ \frac{d[X_{2}T_{2}]}{dt} $= k_{1} [X_{2}] [T_{2}] - \frac{k_{1}}{a_{5}} [X_{2}T_{2}] - \frac{k_{cat}}{26 m c} [X_{2}T_{2}] - k_{1} [X_{2}T_{2}] [X_{3}] + \frac{k_{1}}{a_{6}} [X_{3}X_{2}T_{2}]\notag\\ \frac{d[Z_{2}]}{dt} $= k_{2N} [W_{2}] - \frac{k_{cat}}{24 m c} [Z_{2}] + \frac{k_{cat}}{2 m c} [X_{3}X_{2}T_{2}]\notag\\ \frac{d[W_{2}]}{dt} $= \frac{k_{cat}}{26 m c} [X_{2}T_{2}] - k_{2N} [W_{2}] + \frac{k_{cat}}{24 m c} [Z_{2}]\notag\\ \frac{d[X_{3}X_{2}T_{2}]}{dt} $= k_{1} [X_{2}T_{2}] [X_{3}] - \frac{k_{1}}{a_{6}} [X_{3}X_{2}T_{2}] - \frac{k_{cat}}{2 m c} [X_{3}X_{2}T_{2}] + k_{1} [X_{2}] [X_{3}T_{2}] - \frac{k_{1}}{a_{8}} [X_{3}X_{2}T_{2}]\notag\\ \frac{d[X_{3}T_{2}]}{dt} $= k_{1} [T_{2}] [X_{3}] - \frac{k_{1}}{a_{7}} [X_{3}T_{2}] - k_{1} [X_{2}] [X_{3}T_{2}] + \frac{k_{1}}{a_{8}} [X_{3}X_{2}T_{2}]\notag\\ \frac{d[X_{3}]}{dt} $= \frac{k_{cat}}{24 m c} [Z_{2}] - k_{1} [T_{2}] [X_{3}] + \frac{k_{1}}{a_{7}} [X_{3}T_{2}] - k_{1} [X_{2}T_{2}] [X_{3}] + \frac{k_{1}}{a_{6}} [X_{3}X_{2}T_{2}] + \frac{k_{1}}{a_{9}} [X_{3}T_{3}] - k_{1} [X_{3}] [T_{3}] - k_{1} [X_{3}] [X_{4}T_{3}] + \frac{k_{1}}{a_{12}} [X_{4}X_{3}T_{3}]\notag\\ \frac{d[T_{3}]}{dt} $= - k_{1} [T_{3}] [X_{3}] + \frac{k_{1}}{a_{9}} [X_{3}T_{3}] - k_{1} [T_{3}] [X_{4}] + \frac{k_{1}}{a_{11}} [X_{4}T_{3}]\notag\\ \frac{d[X_{3}T_{3}]}{dt} $= k_{1} [X_{3}] [T_{3}] - \frac{k_{1}}{a_{9}} [X_{3}T_{3}] - \frac{k_{cat}}{26 m c} [X_{3}T_{3}] - k_{1} [X_{3}T_{3}] [X_{4}] + \frac{k_{1}}{a_{10}} [X_{4}X_{3}T_{3}]\notag\\ \frac{d[Z_{3}]}{dt} $= k_{2N} [W_{3}] - \frac{k_{cat}}{24 m c} [Z_{3}] + \frac{k_{cat}}{2 m c} [X_{4}X_{3}T_{3}]\notag\\ \frac{d[W_{3}]}{dt} $= \frac{k_{cat}}{26 m c} [X_{3}T_{3}] - k_{2N} [W_{3}] + \frac{k_{cat}}{24 m c} [Z_{3}]\notag\\ \frac{d[X_{4}X_{3}T_{3}]}{dt} $= k_{1} [X_{3}T_{3}] [X_{4}] - \frac{k_{1}}{a_{10}} [X_{4}X_{3}T_{3}] - \frac{k_{cat}}{2 m c} [X_{4}X_{3}T_{3}] + k_{1} [X_{3}] [X_{4}T_{3}] - \frac{k_{1}}{a_{12}} [X_{4}X_{3}T_{3}]\notag\\ \frac{d[X_{4}T_{3}]}{dt} $= k_{1} [T_{3}] [X_{4}] - \frac{k_{1}}{a_{11}} [X_{4}T_{3}] - k_{1} [X_{3}] [X_{4}T_{3}] + \frac{k_{1}}{a_{12}} [X_{4}X_{3}T_{3}]\notag\\ \frac{d[X_{4}]}{dt} $= \frac{k_{cat}}{24 m c} [Z_{3}] - k_{1} [T_{3}] [X_{4}] + \frac{k_{1}}{a_{11}} [X_{4}T_{3}] - k_{1} [X_{3}T_{3}] [X_{4}] + \frac{k_{1}}{a_{10}} [X_{4}X_{3}T_{3}] - k_{1} [T_{4}] [X_{4}] + \frac{k_{1}}{a_{13}} [X_{4}T_{4}]\notag\\ \frac{d[T_{4}]}{dt} $= - k_{1} [T_{4}] [X_{4}] + \frac{k_{1}}{a_{13}} [X_{4}T_{4}]\notag\\ \frac{d[X_{4}T_{4}]}{dt} $= k_{1} [X_{4}] [T_{4}] - \frac{k_{1}}{a_{13}} [X_{4}T_{4}] - \frac{k_{cat}}{36 m c} [X_{4}T_{4}]\notag\\ \frac{d[W_{4}]}{dt} $= \frac{k_{cat}}{36 m c} [X_{4}T_{4}]\notag\\