Mathematical modeling is a powerful tool to validate and model the kinetics of different systems in synthetic biology. Herein, we are using mathematical modeling to detect the increased level of phenylalanine (phe) in phenylketonuria patients in our diagnostic platform. It depends on a whole cell-based biosensor through a cascade of reactions to finally end by formation of β-galactosidase that turns the color into blue once bound to its substrate (X-gal) as mentioned in figure (1). We used populations of transcription activator (TA), transcription activator DNA (TADNA) ,RNA ,DNA , and beta galactosidase (B) in this model.

Figure (1) represents the cascade of reactions in whole cell-based biosensor model.

dTA

dt

= K1*TA - K2*TADNA - K4*TA + K3*TADNA

This equation represents that phenylalanine enters the cell through permease ,therefore activating transcription activator (TA) that forms at a rate of K1 and a rate of formation of K3 when TADNA dissociates. In addition to, degradation at a rate of K4 and bind to DNA to form TADNA at of rate K2 as mentioned in figure (2)

Figure (2) illustrates TA formation resulted from phenylalanine in the media and TADNA formation through combination of DNA and TA

d TADNA

dt

= K2*TADNA - K5*TADNA - K3*TADNA

The second equation is about the formation of TADNA depending on a formation rate of K2 by association of transcription activator (TA) and DNA and disappears at a rate of K5 (formation of RNA from DNA by transcription) and at a rate of K3 when TADNA dissociates into transcription activator (TA) and DNA as mentioned in figure (3)

Figure (3) represents the formation of TADNA and transcription of DNA into RNA

dDNA

dt

= K3*TADNA - K2*TADNA

This equation represents the formation of DNA at a rate of K3 when TADNA dissociates into transcription activator (TA) and DNA. It disappears at a rate K2 in which TADNA forms as mentioned in figure (4)

Figure (4) shows combination of (TA) and DNA to form TADNA and its dissociation to from DNA

dRNA

dt

=K5*TADNA - K6*RNA - K7*RNA

Equation (4) is about formation of RNA depending on a formation rate of K5 (formation of RNA from DNA by transcription) and disappears at a rate of K7 (RNA degradation) and K6 for RNA to be translated into beta galactosidase as mentioned in figure (5)

Figure (5) illustrates Formation of RNA and translation of RNA into beta-galactosidase

dB

dt

= K6*RNA - K8*B

This equation represents Beta-galactosidase formation (B) depending on a formation rate of K6 in which RNA would be translated into beta galactosidase and degrades at a rate of K8 as mentioned in figure (6)

Figure (6) illustrates formation and degradation of beta-galactosidase.

Graph (1) represents level of the biomarker (phenylalanine) as it enters the whole cell-based biosensor untill it reaches the steady state after about 30 time units

Graph (2) represents level of the beta-galactosidase as it releases in consistance with phenylalanine level untill it reaches the steady state after about 30 time units

Graph (3) represents the relationship between the biomarker and the released beta-galactosiase

Graph(3) illustrates a direct relation between biomarker and beta-galactosidase ,so as the biomarker increases, the released amount of beta-galactosidase increases till it reaches constant value after about 30 time units. Therefore, the maximum amount of the biomarker releases the maximum amount of beta-galactosidase.

This model is to simulate the kinetics of the riboswitch (L7Ae with kink turns) that is used in our circuit. The designed circuit is to detect if the increasing substance is either phenylalanine or tyrosine via TyrR. So if phenylalanine level is elevated, L7Ae is formed as it is downstream TyrR as shown in figure (7) that forms a complex via binding with its kink-turn on another circuit; that complex inhibits expression of cas12g as shown in figure (8), so the circuit will be able to express lacZ alpha (beta-galactosidase) in diagnostic circuit or PAH in the therapeutic circuit. If tyrosine level is elevated with a decreased level of phenylalanine, it activates tyrR inhibitory promoter so no L7Ae would be expressed resulting in cas12g expression to control the circuit. we used these set of populations in our model :mRNA1on (O), L7Ae (A) , mRNA2off (B) , mRNA2on (T) , PP7 (R) , mRNA3on (Q) , mRNA3off (M) , and mRNA1off (V) as shown in figure (9).

Figure(7) The first diagnostic circuit design including L7Ae downstream TyrR

Figure(8) The second diagnostic circuit design including kink-turns and cas12g

Figure (9) illustrates the kinetics of all reactions in riboswitch model

dO

dt

= cr1*O - dmRNA1on * O

The first equation represents the formation of this mRNA (O) depending on a formation rate of cr1 and its dissociation rate is dmRNA1on as shown in figure (10)

dA

dt

= PL7Ae * O - dL7Ae * A

The rate of formation of L7Ae (A) depends on a formation rate of PL7Ae and its dissociation rate is dL7Ae as shown in figure (10)

Figure (10) describes the kinetics of mRNA (O)

dB

dt

= cr2*B - dmRNA2off * B

The rate of formation of mRNA (B) depends on a formation rate of cr2 and its dissociation rate is dmRNA2off as shown in figure (11)

Figure (11) illustrates formation and degradation of (B) mRNA

dT

dt

= K2 * A *B - dmRNA2on * T

The rate of formation of (T) depends on a formation rate of K2 and its dissociation rate is dmRNA2on as shown in figure (12)

dR

dt

= pP * T - dP * P

The rate of formation of (P) depends on a formation rate of pP and it dissociates at a rate of dP as shown in figure (12)

Figure (12) formation of mRNA (T) and its translation into protein (R)

dQ

dt

= cr3*Q - dmRNA3on * Q

The rate of formation of mRNA (Q) depends on a formation rate of cr3 and this mRNA dissociates at a rate of dmRNA3on as shown in figure (13)

Figure (13) illustrates formation and degradation of mRNA (M)

dM

dt

= K1 * A * Q - dmRNA3off * M

The rate of formation of mRNA (M) depends on a formation rate of K1 and this mRNA dissociates at a rate of dmRNA3off as shown in figure (14)

Figure (14) illustrates formation and degradation of mRNA (M)

dV

dt

= K3 * P * O - dmRNA1off * V

The rate of formation of mRNA (V) depends on a formation rate of K3 and this mRNA dissociates at a rate of dmRNA1off as shown in figure (15)

Figure (15) formation of mRNA (v) from complex (RO)

dP1

dt

= PP1*P1 - d1 * p1

dP2

dt

=PP2*P2 - d2 * p2

The rates of formation of protein 1 (P1) and protein 2 (P2) depend on formation rates of PP1 and PP2 ,respectively and they dissociate at rates of p1 and p2 ,respectively.

Graph (4) represents the kinetics of the used riboswitch

Graph (4) illustrates riboswitch kinetics in which Q represents the condition where L7Ae is expressed and bound to its kink-turns ,therefore inhibiting the expression of cas12g. However, M represents no expression of L7Ae in which cas12g would be expressed to control the circuit if the phenylalanine is absent.

In order to simulate the kinetics of aptamer binding with phenylalanine on lateral flow assay, we have used a set of ordinary differential equations (ODEs) using different populations including: free analyte concentration (A), initial analyte concentration (Ao), free detector concentration (P) , initial detector concentration (Po) , analyte diffusivity (DA), detector diffusivity (DP), flow rate (U), length of test strip (L), position of test strip(dx) ,initial free aptamer concentration (Ro) , and free aptamer concentration (R) as shown in figure (16).

Figure (16) illustrates the kinetics of the aptamer binding with its target on LFA

dA

dt

=DA*(A+A-1-2*A)/dx) - (vc/dx) * (A-A-1) / 2*dx) - (PA-RA)

The first equation represents the kinetics of the analyte that is dependant on flow rate U and its value decreases along test strip by . The concentration of analyte decreases until via binding to the detector and the Aptamer to form PA and RA , as shown in figure (17).

Figure (17) represents the kinetics of the analyte (A), aptamer-analyte complex (RA), and aptamer-detector-analyte complex (RPA)

dP

dt

DP*(P+P-1-2*A / dx) - (vc/dx) * (P-P-1) / 2*dx) - (PA-RA)

This equation is about the detector in which the concentration of detector decreases as it forms PA and RA and depends on flow rate along test strip by , as shown in figure (18)

dPA

dt

=DP*(PA+PA-1-2*A)/dx) - (vc/dx) * (PA-PA-1) / 2*dx) - (PA-RA)

The equation represents the binding between the detector and the analyte as the concentration of the detector-analyte complex decreases until disappears to form P and A. It also depends on flow rate along test strip by, as shown in figure (18)

Figure (18) represents binding of the analyte to the detector (PA)

dRA

dt

=-(fPA-fRPA1)

The rate of formation of the aptamer-analyte complex decreases until it disappears to form R and A. It depends on flow rate along test strip by , as shown in figure (17)

dPRA

dt

=fPA-fRPA2

The rate of formation of the detector-aptamer-analyte complex decreases until it disappears to form A, R, and P. It depends on flow rate along test strip by , as shown in figure (16)

Graph (5) illustrates the relation between analyte (A) and aptamer analyte complex (RA)

Graph (5) represents plotting both analyte in free state (A) and aptamer analyte complex (RA) after the binding occurs between them. It’s shown that the analyte (phenylalanine) is loaded to the test strip through the sample at the beginning. (A) shows gradual decrease of analyte over the time units. On the other hand, There’s a gradual increase in the aptamer-analyte complex over the time units via binding between the aptamer with its target.This model was beneficial in the optimization of aptamer to analyte concentration that is essential in constructing the consumption line in our lateral flow assay (LFA).

Sequence :CTCTCGGGACGACCGCGTTTCCCAAGAAAGCAAGTATTGGTTGGTCGTCCC

G value:-23.19kcal/mol

Sequence :CTCTCGGGACGACCGGTGGGGGTTCTTTTTCAGGGGAGGTACGGTCGTCCC

G value:-25.67 kcal/mol

Sequence :CTCTCGGGACGACGAGGCTGGATGCATTCGCCGGATGTTCGATGTCGTCC

G value:-19.00 kcal/mol