We are dedicated to develop a model specific to our project, allowing the wet lab to adjust experimental conditions based on predictions generated by the model. Our model primarily consists of two main components, each further divided into four parts: RT-RPA, RCA, PCRD, and the gold nanoparticle diffusion model. We are concurrently developing all four models and will subsequently combine the RT-RPA model with the PCRD model and the RCA model with the gold nanoparticle diffusion model, resulting in the comprehensive models.

Before entering the model part, we first need to analyze the stability of linear and circular structure by diving in the free energy aspect.

    Purposes

By comparing the free energy of linear and circular RNA, we can know the stability of two different structures. Furthermore, we can know the most ideal target to build up the further investigation.

    Methods

We use the RNAeval Webserver in the Vienna RNA to estimate the free energy from the known sequences. In addition, we also separate the linear and the circular structure by using the advanced options in it. We list out the result as follows:

		hsa_circ_0101802	hsa_circ_0004771
Linear	Circular	Linear	Circular
-51.4	-43.3	-42.6	-37.2

Sheet 1. Free energy of different sequences and structures (Units: kcal/mol)

    In addition, the structure of each is given as follows:

hsa_circ_0101802	Linear	Circular

We have designed two methods for detecting the target cicRNAs: “Gold nanoparticle-based colorimetric assay” and “lateral flow rapid screening assay”.

Sheet 2. The structure of hsa_circ_0101802

hsa_circ_0004771	Linear	Circular

Sheet 3. The structure of hsa_circ_0004771

    Results

Based on the aforementioned, we can know that the free energy of two different structures do not have significant difference, which mean they are both stable in the free energy aspect. By taking the features of the 3 end and the 5 end will not easily be separated in the circular form, we choose circular RNA as our target to use in the following investigation.

    Purpose

The establishment of the RT-RPA model is to assist Wetlab in determining whether the current experimental design can yield the expected results through simulation before the actual experiment.

    Assumptions

a. UvsY forms a stable hexamer and binds to UvsX prior to the start of the test.
b. Gp32 binding must proceed with recombinase binding.
c. all primers exist as this complex at t= 0s with negligible error.
d. perfect symmetry between the reactions involving forward and reverse primers.
e. fluid is well mixed and all reactions are kinetically controlled.
f. RPA is a 50 ul reaction

    Equation

We use SimBiology to construct the RT-RPA model. Therefore, we need to make equations for the crucial reactions in RT-RPA and input them into SimBiology. In order to construct the model more effectively, we divide the RT-RPA process into several steps. Firstly, we need to make an equation for the binding of the forward primer (F) and the Single-stranded binding protein (Gp32), here represented by G.

G + F ⇆ FG (1)

Next, the recombinase complex (UvsX.UvsY complex), here represented by R. Displacing Gp32 from the primer in a four-step process:

FG + R ⇆ FGR ' (2a)

FGR ' → FGRi + G (2b)

FGRi + R ⇆ FGRc (2c)

FGRc → FRc + G (2d)

Once recombinase has bound a primer, it seeks out homology in the DNA template. It then inserts the primer into the DNA and dissociates from both by hydrolyzing ATP.

In Reaction (3), the recombinase filament escorts the primer to the region of homology in the DNA template, here represented by D.

FRc + D ⇆ FRcD (3)

Next, hydrolysis of ATP drives the dissociation of the recombinase from the primer/DNA complex, FD. According to the paper, we find that UvsX has two modes of ATP hydrolysis. It can hydrolyze ATP to either AMP or ADP. So we use two equations for two modes of ATP hydrolysis.

FRcD + ATP -> FD + R + AMP + PPi (4a)

FRcD + ATP -> FD + ADP + H3PO4 (4b)

After the yielding of the primer/DNA complex, the Polymerase can bind reversibly to the DNA primer complexes.

FD + P ⇆ PFD (5)

Once the tertiary complex has been formed polymerase can extend the primer to form a new strand of DNA using nucleotides (dNTPs).

PFD + dNTPs → P + PPi + D (6)

Up to this point, we have demonstrated how we constructed the RPA model. Now, by incorporating the step of using reverse transcriptase to reverse transcribe circRNA into DNA templates, we can further refine it into the RT-RPA model.

CircRNA + E → D (7)

This is what the complete RT-RPA model looks like in simbiology:

Figure 1. RT-RPA reaction pathway in MATLAB Simbiology

    Parameter

Variable	Reaction	Value	Unit
k1f	Forward rate constant for (1)	500,000,000	1/(molarity*second)
k1r	Reverse rate constant for (1)	5,000	1/second
keq2af	Forward rate constant for (2a)	60,460,000	1/(molarity*second)
keq2ar	Reverse rate constant for (2a)	1.47	1/second
k2b	Rate constant for (2b)	0.047	1/second
keq2cf	Forward rate constant for (2c)	60,460,000	1/(molarity*second)
keq2cr	Reverse rate constant for (2c)	33	1/second
k2d	Rate constant for (2d)	0.0046	1/second
k3f	Forward rate constant for (3)	60,460,000	1/(molarity*second)
k3r	Reverse rate constant for (3)	59.37	1/second
k4acat	Forward rate constant for (4a)	4.22	1/second
k4bcat	Forward rate constant for (4b)	8.32	1/second
k5f	Forward rate constant for (5)	12,000,000	1/(molarity*second)
k5r	Reverse rate constant for (5)	0.06	1/second
k6f	Forward rate constant for (6)	87	1/second

Sheet 4. the parameter of RT-RPA model

    Initial Values of Species

Variable	Species	Value	Unit
R	Recombinase complex (UvsX.6*UvsY)	5.9e-6	molarity
F	Forward primer	4.8e-7	molarity
G	Single stranded DNA binding protein (Gp32)	9.4e-7	molarity
FG	Forward primer/Gp32 complex	0	molarity
FGR'	Unstable Forward primer/Gp32/ Recombinase complex	0	molarity
FGRi	Forward primer/Gp32/ isolated Recombinase complex	0	molarity
FGRc	Forward primer/Gp32/ cluster Recombinase complex	0	molarity
FRc	Forward primer/ cluster Recombinase complex	0	molarity
FRcD	Forward primer/ cluster Recombinase/DNA template complex	0	molarity
D	DNA template	1.7391e-8	molarity
FD	Forward primer/DNA complex	0	molarity
PFD	Polymerase/forward primer/DNA complex	0	molarity
PPi	Pyrophosphate	0	molarity
P	Polymerase	1.34e-6	molarity
dNTPs	Nucleotides	2.00e-4	molarity

Sheet 5. the initial value of RT-RPA model

    Result

Due to the lack of reaction constants for reverse transcriptase, we are currently only able to simulate a portion of the RPA in the RT-RPA model. We use the initial values provided by Wetlab for simulation and set the reaction time to 20 minutes. The simulation results show the quantity of cDNA rapidly increases in the early stages and levels off in the later stages of the reaction, which is consistent with the result in the paper. Therefore, we proceed to conduct further simulations.

Figure 2. Simulation Results of cDNA after a 20 minutes RPA Reaction

According to the paper, similar concentrations should be reached at the same reaction time, even if the initial values of cDNA are different. Therefore, we set the quantity of cDNA to 1/10 and 1/100 of the original for simulations. The simulation results indicate that, even with different initial cDNA amounts, similar concentrations are reached after 20 minutes of reaction. This not only reaffirms the result in the paper but also instills our confidence in the accuracy of this model.

Figure 3. Simulation Results of different initial values of cDNA after a 20 minutes RPA Reaction

After gaining confidence in this model, we want to apply this model in practical scenarios. Wetlab suggests that the quantity of the forward primer may influence the lower limit of the initial values of cDNA. Therefore, we decided to use this model to verify this proposition, setting the quantity of the forward primer to 1/10 and the reaction time to 40 minutes for simulations. The simulation results indicate that when the quantity of the forward primer is reduced to 1/10 of the original, the lower limit of the initial values of cDNA changes from 10^-33 to 10^-15, confirming our proposition.

Figure 4. simulation results of cDNA after a 40 minutes RPA reaction when the forward primer quantity is 4.8e-7

Figure 5. simulation results of cDNA after a 40 minutes RPA reaction when the forward primer quantity is 4.8e-8

    Conclusion

Although the lack of reaction constants for reverse transcriptase restricts our simulations to only a portion of the RPA in the RT-RPA model, through the three simulations mentioned above, we have gained considerable confidence in the RT-RPA model. In the future, we are committed to obtaining the reaction constants for reverse transcriptase through experiments to complete the RT-RPA model. After confirming its accuracy, we will integrate it with the PCRD model.

    Purpose

The establishment of the RCA model is to assist Wetlab in determining whether the current experimental design can yield the expected results through simulation before the actual experiment.

    Assumptions

a. The degradation rate of cDNA is negligible. b. Once the probe binds to cDNA, it won 't dissociate. c. The ratio of inner binding sites to outer binding sites is 1:1

We use SimBiology to construct the RCA model. Therefore, we need to make equations for the crucial reactions in RCA and input them into SimBiology. In order to construct the model more effectively, we divide the RCA process into several steps. First, We need to make an equation for the reverse transcription of circRNA (C) to cDNA (D) for Protoscript II(E).

C+E→D (1)

Next, the ssDNA probe(D) on the gold nanoparticles binds to the binding sites on cDNA. According to the paper, we have determined that the hybridization of ssDNA is a second-order reaction. Therefore, we can construct the following equation:

P+D→H (2)

However upon further investigation of other papers, we discovered that the hybridization of ssDNA is influenced by steric hindrance, and the reaction constants can vary up to 42 times. Therefore, we have divided equation (2) into (2a) and (2b) to distinguish between outer binding sites and inner binding sites.

Probe+D(o) →H (2a)

Probe+D(i) →H (2b)

Finally, we need to consider the dissociation of cDNA and circRNA, leading to the following equation:

C→degradation (3)

D→degradation (4)

This is what the complete RCA model look like in simbiology:

Figure 6. RCA reaction pathway in MATLAB Simbiology

    Result

Due to a lack of crucial parameters and initial value for some reactants, we can only present our ordinary differential equations (ODEs) at this point. However, once we collect the relevant data and parameters we need in the future, we will execute our model and complement it with experiments to validate its accuracy.

The ODEs of RCA:

    Purposes

As to predict the color and further diagnosis the symptoms, we have to analyze the reason of changing color and to quantify the changing extent.

    Assumptions

a. Every particle is a sphere
b. Every particle has same magnitude against each other

    Working principles

The extent of changing color is determined by coupling intensity between particles, which can be quantified as:

where the schematic diagram can be given as:

Figure 7. Schematic diagram of nano gold particles

    Results

Before calculation, based on the refer study, the given equation only applied when the radius of the particle is 2 to 27nm. Since that the diameter of the nano gold particle is 13nm and the absorption wavelength is 520nm, we can know that the possible color change is:

Figure 8. Changing color range

However, the connection of cDNA will increase the magnitude of the particle, which also change the absorption wavelength, the more it connects the more the wavelength will increase. Hence, there may be some deviation between the prediction and the actual condition.

    Purposes

In order to simulate the diffusion of nano gold complexes or even further predict the color change, we use the method below to observe how does the particles diffuse , change under different concentration and the whole diffusion process, etc.

    Assumptions

a. Particles are placed in the middle of Eppendorf
b. No pressure difference in both ends
c. Each particle occupy same volume
d. cDNA will 100% connect to particles
e. Turbulent kinetic energy is proportional to the amount of particles

    Method

We use the function in Simscale to simulate the particles diffusion in the liquid under equilibrium flow condition. To have a precise simulation, we first draw the geometric characteristics in CAD to show that the diffusion only happen in a small portion of Eppendorf since the liquid volume is quite small.

Figure 9. Geometric illustration of small portion of Eppendorf

We then input the necessary physics variables to retrieve the simulation results, the variables are listed as follows:

Parameter	Notation	Value	Units
Kinematic viscosity	v	9.338 × 10-7	m2/s
Density	ρ	1.0	kg/m3
Velocity inlet	Vz	-0.01t + 0.01	m/s
Velocity outlet	Uz	0	m/s
Pressure inlet	P	1.013 × 105	Pa

Sheet 6. Physics variables in Simscale simulation

Since the concentration difference of the two ends may lead to the diffusion be influenced by osmotic pressure, therefore we take that into consideration and have the further derivation:

From the given equations, we can know that the relationship between velocity is proportional to the coefficient of time. Base on the input value and the geometric characteristics, we can run the simulation of diffusion as follows:

Animation 1. Diffusion process of the particles

In addition, we can list out the change of value along with time:

Figure 10. Physics value change along with time

From the animation above, we can take a closer look at the cutting plane to analyze the distribution of the particles by the extent of turbulent kinetic energy. The graph is given as:

Figure 11. Diffusion under time-depending velocity

For the distance calculation, we assume that every particle is homogeneously distributed in the liquid, thus occupy for the same volume of the liquid. On the other hand, since the turbulent kinetic energy is proportional to the amount of particles, we then use the turbulent energy of each data to represent the amount of particles account for every part of the area. Furthermore, we can also calculate the volume of the area by using the amount of particles. In the end, we assume that every area is a cubic, which the volume is side length cubed.

From the previous description, we can set up the maximum and minimum distance of the particles under different concentration in the below graph:

Figure 12. The trend of distance under different concentration

From Figure 10., we can know that when the concentration increase, then the amount of bonding cDNA will increase. Therefore, we can know that when the particles connect to more cDNA, the distance between particles will increase, which makes the particles in the liquid become relatively not segregated. Overall, it further indicates that when the amount of virus increase, base on the condition of relatively not segregated, the color will become red. The result fit with the experiment, vice versa.

Moreover, the difference between the given concentration are given as the below animations:

0.1	1	10	100

Animation 2. Diffusion under different concentration

    Discussion

The reason why the distance is not 27nm can be contributed to the following reasons:

1. Simscale can only one diffusion repeatedly, which differs from the actual condition that the diffuse last for a period of time.
2. We only simulate one fluid which cannot consider the other that also contributed to the diffusion condition.
3. The gravity is neglected, which may influence the diffusion along y axis.
4. The area is not fully occupied by the particles.
5. The boundary condition of the fluid will influence the particles that near the boundary.

    Improvement

In order to make the diffusion model able to determine which part of the fluid may have color after diffusion, we may use the solid-fluid interaction module in Comsol to simulate the whole procedure and further determine the distance between particles.

Geometric characteristics:

Shape	Upper radius (m)	Lower radius (m)	Height (m)
Cone	0.025	0.005	0.0279

Sheet 7. Geometric characteristics of modeling

Afterward, we input the above information into Comsol we can know that:

Figure 13. Modeling of fluid part in the container

In the end, we need to input the characteristics of both fluid and the solid to simulate the diffusion result to attain the distance between every particle.

    Purposes

To understand how each parameter affects the signal at the test line and further optimize the design computationally, we design a Later Flow Assay model. This model is used in iGEM team Rochester(2020), which adapted from the model created by Liu et al (2017). In addition, since we have two test lines and several different ingredients in our experiment, we add a factor of the composition of the buffer to determine the following procedures. We also add the RT-RPA model and combine it with PCRD model which generated in Symbiology to make it more workable. More assumptions are elaborated on as follows:

    Assumptions

a. Analytes and detectors do not react before entering the test line.
b. Kinetic rate constant are constant with concentration, which do not alter with binding.
c. No loss of analytes and detectors during the process.
d. The solution contain analytes, detectors and receptors are well-mixed.
e. One analyte can only be captured by one detector.
f. Two complexes will not influence each other.

    Glossary

Variable	Definition	Units
𝐴_𝐷	Free analyte (DIG) concentration	nM
𝐴_𝐷0	Initial analyte concentration (DIG) in sample
𝐴_𝐹	Free analyte (FAM) concentration	nM
𝐴_𝐹0	Initial analyte concentration (FAM) in sample
𝑃	Free detector concentration
𝑃_0	Initial free detector concentration
𝑅	Free receptor concentration
𝑅_0	Initial free receptor concentration
〖𝐷_𝐴〗_𝐷	Analyte (DIG) diffusivity	mm²/s
〖𝐷_𝐴〗_𝐹	Analyte (FAM) diffusivity	mm²/s
𝐷_𝑃	Detector diffusivity	mm²/s
𝑈	Flow rate	mm/s
𝐿	Length of test strip	mm
𝑥	Position on test strip
𝑡_1	Time since sample entered test line 1	s
𝑡_2	Time since sample entered test line 2	s
𝐾_(𝑎_𝐷1 )	Detector-analyte (DIG) association constant	nM⁻¹s⁻¹
𝐾_(𝑑_𝐷1 )	Detector-analyte (DIG) dissociation constant	s⁻¹
𝐾_(𝑎_𝐹1 )	Detector-analyte (FAM) association constant	nM⁻¹s⁻¹
𝐾_(𝑑_𝐹1 )	Detector-analyte (FAM) dissociation constant	s⁻¹
𝐾_(𝑎_𝐷2 )	Receptor-analyte (DIG) association constant	nM⁻¹s⁻¹
𝐾_(𝑑_𝐷2 )	Receptor-analyte (DIG) dissociation constant	s⁻¹
𝐾_(𝑎_𝐹2 )	Receptor-analyte (FAM) association constant	nM⁻¹s⁻¹
𝐾_(𝑑_𝐹2 )	Receptor-analyte (FAM) dissociation constant	s⁻¹

Sheet 8. Correspondence of symbol

    Analysis:

Unfortunately, since most of the value cannot be found in any given information, we refer to the numeric value from the igem 2020 Rochester team. Due to the reason of the related value is similar to the refer study, we assume that the error range is 10X~0.1X. Then we use this assumption to do the following sensitivity analysis:

Figure 14. Sensitivity analysis of parameters

To be more specific, the comparison between two major factors-the forward reaction rate and the backward reaction rate can be elaborated as:

Figure 15. Sensitivity comparison between forward and backward reaction rate

Based on the result shows in Figure 15, we can know that the most influential factor in the PCRD analysis is the forward reaction rate. Hence, we change forward reaction rate and lists below:

Figure 16. RPA under different forward reaction rate (no backward reaction rate)

From the graph we can clearly observe the significant influence of how forward reaction rate change the RPA along with time. It is clear that the larger the forward reaction rate, the faster the RPA will increase, which make the target concentration increase and further affect the whole experiment. Then, we take backward reaction rate into consideration:

Figure 17. RPA under different forward reaction rate

From the graph we can see that the backward reaction rate has a relatively larger influence on the case which has smaller forward reaction rate. This is because when the magnitude of two major factors are close, then it is hard to produce a net reaction which help the experiment can proceed. Worth mentioning, after reaching the peak value of the RPA, the backward reaction rate will also decrease the concentration.

    Analysis:

In order to calculate the concentration of different components along with time, we have to write out the overall equations as follows. Worth mentioning, the reaction is different than the normal PCRD model since we have 2 test lines. Therefore, we have to calculate the reactions twice and assume that the reactants will not react with each other during the flow.

*Note: Ka1& Kd1 correspond to the reaction that “+P” ; whereas Ka2& Kd2 correspond to the reaction that “+R”

Similarly, the rate of formation of each is the association of the molecules minus the dissociation of complexes. Still, we consider two different kinds of conditions since we have 2 test lines.

By taking the advection and the diffusion into consideration, the listed partial differentiation equations can help us to simulate the concentration of each component along with the time change:

**Note: The “A” below include A_D & A_F

In conclusion, the 2 test lines will influence each other result. Therefore, during the calculation we need to add some items in the functions(as the above has shown) to get the final result. By taking the influence of 2 kinds of component, we can predict the parameters that we want.

Reference

[1] Liu, Jie ; Berger, Christopher L ; Morrical, Scott W (2013). Kinetics of Presynaptic Filament Assembly in the Presence of Single-Stranded DNA Binding Protein and Recombination Mediator Protein. Biochemistry (Easton), 2013, Vol.52 (45), p.7878-7889 https://pubs.acs.org/doi/full/10.1021/bi401060p

[2] Moody, Clint ; Newell, Heather ; Viljoen, Hendrik (2016). A mathematical model of recombinase polymerase amplification under continuously stirred conditions. Biochemical engineering journal, 2016, Vol.112, p.193-201 https://www.sciencedirect.com/science/article/abs/pii/S1369703X16301152?casa_token=nv3pO9UHxvYAAAAA:9kA9Wc27RzRSbEKsF5uC_q44f_6Z-g4MGT-T63K7gns8uvCOv2T9ttnD-UmuxCmZKDGIWak_roM#bib0090

[3] Farb, Joshua N. ; Morrical, Scott W. (2009). Role of Allosteric Switch Residue Histidine 195 in Maintaining Active-Site Asymmetry in Presynaptic Filaments of Bacteriophage T4 UvsX Recombinase. Journal of molecular biology, 2009, Vol.385 (2), p.393-404 https://www.sciencedirect.com/science/article/pii/S0022283608014125

[4] Takeda, Yoshihiro ; Kondow, Tamotsu ; Mafuné, Fumitaka (2008). Hybridization of ssDNA with a Complementary DNA Probe Tethered to a Gold NanoparticleEffect of Steric Hindrance Caused by Conformation. Journal of physical chemistry. C, 2008, Vol.112 (1), p.89-94 https://pubs.acs.org/doi/10.1021/jp0665980

[5] Xu, Jun ; Craig, Stephen L (2005). Thermodynamics of DNA Hybridization on Gold Nanoparticles. Journal of the American Chemical Society, 2005, Vol.127 (38), p.13227-13231 https://pubs.acs.org/doi/full/10.1021/ja052352h

[6] Chen, Chunlai ; Wang, Wenjuan ; Ge, Jing ; Zhao, Xin Sheng (2009). Kinetics and thermodynamics of DNA hybridization on gold nanoparticles. Nucleic acids research, 2009, Vol.37 (11), p.3756-3765 https://academic.oup.com/nar/article/37/11/3756/1087825?login=false

[7] Prashant K. Jain, Wenyu Huang, and Mostafa A. El-Sayed (2007). On the Universal Scaling Behavior of the Distance Decay of Plasmon Coupling in Metal Nanoparticle Pairs: A Plasmon Ruler Equation. Nano Letters. Vol. 7, No. 7 2080~2088. https://pubs.acs.org/doi/10.1021/nl071008a

[8] R. Tsekov, P. Georgiev, S. Simeonova and K. Balashev (2017). Quantifying the Blue Shift in the Light Absorption of Small Gold Nanoparticles. Sci. 70 (2017) 1237-1246. https://arxiv.org/ftp/arxiv/papers/1702/1702.04513.pdf

[9] Wu R, Peng H, Zhu J-J, Jiang L-P and Liu J (2020). Attaching DNA to Gold Nanoparticles With a Protein Corona. Front. Chem. 8:121. doi: 10.3389/fchem.2020.00121

[10] Cheng-Chu Wu (2008). Application of Metal Nanoparticles ( MNPs ) Labeling to DNA Microarray. https://tdr.lib.ntu.edu.tw/bitstream/123456789/8832/1/ntu-98-1.pdf

[11] Lewis GN (1908). The Osmotic Pressure of Concentrated Solutions and the Laws of the Perfect Solution. Journal of the American Chemical Society. https://pubs.acs.org/doi/abs/10.1021/ja01947a002

[12] Laurence, D. (2002). Applications of Reynolds Averaged Navier Stokes Equations to Industrial Flows. In van Beeck, J. P. A. J.; Benocci, C. (eds.). Introduction to Turbulence Modelling, Held March 18-22, 2002 at Von Karman Institute for Fluid Dynamics.

[13] Mott, P. H.; Roland, C. M. (2012). Limits to Poisson's ratio in isotropic materials—general result for arbitrary deformation. Physica Scripta. Chemistry Division, Naval Research Laboratory. [doi](https://en.wikipedia.org/wiki/Doi_(identifier)):https://doi.org/10.1088/0031-8949/87/05/055404

[14] iGEM team Rochester (2020). Available from: https://static.igem.org/mediawiki/2020/3/33/T--Rochester--lfam.pdf

[15] Liu, Z., Hu, J., Li, A., Feng, S., Qu, Z., & Xu, F. (2017). The effect of report particle properties on lateral flow assays: A mathematical model. Sensors and Actuators. B, Chemical, 248, 699-707. doi:10.1016/j.snb.2017.04.024

[16] Berli, C. L. A., & Kler, P. A. (2016). A quantitative model for lateral flow assays. Microfluidics and Nanofluidics, 20(7), 1-9. doi:10.1007/s10404-016-1771-9

[17] Linge S., Langtangen H.P. (2017) Advection-Dominated Equations. In: Finite Difference Computing with PDEs. Texts in Computational Science and Engineering, vol 16. Springer, Cham. https://doi.org/10.1007/978-3-319-55456-3_4

[18] U.S. Food & Drug Administration (2013). Sylvant Pharmacology Review. Available from: https://www.accessdata.fda.gov/drugsatfda_docs/nda/2014/125496Orig1s000PharmR.pdf