Dynamic Reaction Analyzer -ODEs

To simulate the detection system built in this project, we established a series of ordinary differential equations(ODEs). By comprehensively and systematically simulating the reactions in the system, we obtained the change tendency of the substances in the system and estimated its operating time, which is beneficial for understanding the mechanism of our system. Meanwhile, our ODEs also helped the hardware team estimate the order of time required for the reaction system to reach equilibrium.

The construction of ODEs can be divided into three parts: AvrRpt2 input part, COR input part, and GFP signal part. We present the reactions, involved substances, their corresponding assumptions, and kinetic equations for each part as follows.

Note: To simplify the kinetic description in our ODE system, the single-step enzymatic cleavage of a linker is represented by two sequential reactions: enzymatic hydrolysis first produces an intermediate dimer (without the linker), which then undergoes dissociation into the two constituent peptide chains.

In addition, based on wetlab analyses indicating that peptides A and B share high sequence and structural homology, we have assumed functional equivalence in their binding kinetics. Accordingly, in our model, their respective binding and dissociation rate constants with other peptide chains are set to be equal.

a. AvrRpt2 Input

AvrRpt2, a protein specifically secreted by Pst, is difficult to detect at low concentrations extracted from a leaf. We therefore designed a positive feedback circuit to amplify its signal, as shown in the figure below. The system is based on two key components: the TEV protease N-terminus fused to peptide A and its inactive C-terminus fused to peptide B, connected by a cleavable linker containing recognition sites for both AvrRpt2 and TEV protease; and the complete TEV C-terminus fused to peptide A′, which is free in solution and exhibits high affinity for peptide A. When AvrRpt2 cleaves the linker, peptide A is released and binds to A′, reconstituting an active TEV protease. This active protease then cleaves additional linker molecules, initiating an autocatalytic amplification cycle that significantly enhances detection sensitivity. The reaction equations and corresponding ODEs for this process are derived as follows. See more details in description part.

Reactions

Assumptions

1. The reaction follows first-order kinetics with respect to both the substrate and the product.
2. The activities of AvrRpt2 and TEV protease do not affect each other.
3. The two peptides and the linker at the middle end of the generated TEV protease do not affect its function (kinetic constant). At the same time, the binding strength of the peptides linked to TEV protease does not take into account the linked part of TEV protease either.
4. Temperature change does not affect the kinetic parameters of each reaction.
5. The substrate of the reaction catalyzed by AvrRpt2 and TEV protease is sufficient to satisfy the Michaelis-Menten equation.
6. The kinetic constants (kcat and Km) of TEV protease for cutting linkerAB and linkerCD are considered as the same (Connected with part C) because the cutting sites of TEV preotease of linkerAB and linkerCD are the same.
7. The binding rate constants of all peptide binding reactions are assumed to be equal, all being 10⁵ M^-1 s^-1.
8. When TEVp cuts both linkerAB and linkerCD simultaneously, the Michaelis-Menten equations of the two parts do not affect each other.
9. The lysis or inactivation of AvrRpt2 and TEV protease in the system is not taken into account.

Variables

Parameters

Note: Kinetic parameters labeled as "MPEK's prediction" were obtained using the MPEK deep learning framework ², which predicts enzymatic reaction kinetic parameters from protein sequences. Parameters labeled as "Prodigy's prediction" were derived from the PRODIGY web server ³, which estimates protein-protein binding affinity. The same below.

Odes

$$ \frac{dc_{\mathrm{linkerAB}}}{dt} = -\frac{k_{\mathrm{cat,AvrRpt}} \cdot c_{\mathrm{AvrRpt}} \cdot c_{\mathrm{linkerAB}}}{K_{\mathrm{m,AvrRpt}} + c_{\mathrm{linkerAB}}} - \frac{k_{\mathrm{cat,TEVp}} \cdot c_{\mathrm{TEVp}} \cdot c_{\mathrm{linkerAB}}}{K_{\mathrm{m,TEVp}} + c_{\mathrm{linkerAB}}} $$

$$ \frac{dc_{\mathrm{AB}}}{dt} = \frac{k_{\mathrm{cat,AvrRpt}} \cdot c_{\mathrm{AvrRpt}} \cdot c_{\mathrm{linkerAB}}}{K_{\mathrm{m,AvrRpt}} + c_{\mathrm{linkerAB}}} + \frac{k_{\mathrm{cat,TEVp}} \cdot c_{\mathrm{TEVp}} \cdot c_{\mathrm{linkerAB}}}{K_{\mathrm{m,TEVp}} + c_{\mathrm{linkerAB}}} - k_{1,\mathrm{AB}} \cdot c_{\mathrm{AB}} + k_{2,\mathrm{AB}} \cdot c_A \cdot c_B $$

$$ \frac{dc_{\mathrm{TEVp}}}{dt} = k_{2,\mathrm{AA'}} \cdot c_A \cdot c_{\mathrm{A'}} - k_{1,\mathrm{AA'}} \cdot c_{\mathrm{TEVp}} $$

$$ \frac{dc_A}{dt} = k_{1,\mathrm{AB}} \cdot c_{\mathrm{AB}} + k_{1,\mathrm{AA'}} \cdot c_{\mathrm{TEVp}} - k_{2,\mathrm{AA'}} \cdot c_A \cdot c_{\mathrm{A'}} - k_{2,\mathrm{AB}} \cdot c_A \cdot c_B $$

$$ \frac{dc_{\mathrm{A'}}}{dt} = k_{1,\mathrm{AA'}} \cdot c_{\mathrm{TEVp}} - k_{2,\mathrm{AA'}} \cdot c_A \cdot c_{\mathrm{A'}} + k_{1,\mathrm{A'B}} \cdot c_{\mathrm{A'B}} - k_{2,\mathrm{A'B}} \cdot c_B \cdot c_{\mathrm{A'}} $$

$$ \frac{dc_B}{dt} = k_{1,\mathrm{AB}} \cdot c_{\mathrm{AB}} - k_{2,\mathrm{AB}} \cdot c_A \cdot c_B + k_{1,\mathrm{A'B}} \cdot c_{\mathrm{A'B}} - k_{2,\mathrm{A'B}} \cdot c_B \cdot c_{\mathrm{A'}} $$

$$ \frac{dc_{\mathrm{A'B}}}{dt} = -k_{1,\mathrm{A'B}} \cdot c_{\mathrm{A'B}} + k_{2,\mathrm{A'B}} \cdot c_B \cdot c_{\mathrm{A'}} $$

b. COR Input

In the other detection pathway, we targeted coronamycin (COR), a compound specifically secreted by Pseudomonas syringae. The underlying principle relies on COR's ability to mediate the interaction between JAZ1 and COI1 proteins. We fused JAZ1 to the C-terminal and COI1 to the N-terminal of the plum pox virus (PPV) protease, respectively. The presence of COR brings JAZ1 and COI1 into proximity, enabling the reconstitution of an active PPV protease. Conversely, in the absence of COR, the protease remains inactive. This COR-dependent mechanism forms the basis of the part described below.

Reactions

$$\text{COR} + \text{COI1} + \text{JAZ1} \underset{k_{1,\mathrm{PPVp}}}{\overset{k_{2,\mathrm{PPVp}}}{\rightleftharpoons}} \text{PPVp}$$

Assumptions

Suppose the reaction follows first-order kinetics with respect to both the substrate and the product.

Variables

Parameters

Odes

$$\frac{dc_{\mathrm{PPVp}}}{dt} = k_{1,\mathrm{PPVp}} \cdot c_{\mathrm{COR}} \cdot c_{\mathrm{COI1}} \cdot c_{\mathrm{JAZ1}} - k_{2,\mathrm{PPVp}} \cdot c_{\mathrm{PPVp}}$$

$$\frac{dc_{\mathrm{COI1}}}{dt} = -k_{1,\mathrm{PPVp}} \cdot c_{\mathrm{COR}} \cdot c_{\mathrm{COI1}} \cdot c_{\mathrm{JAZ1}} + k_{2,\mathrm{PPVp}} \cdot c_{\mathrm{PPVp}}$$

$$\frac{dc_{\mathrm{COR}}}{dt} = -k_{1,\mathrm{PPVp}} \cdot c_{\mathrm{COR}} \cdot c_{\mathrm{COI1}} \cdot c_{\mathrm{JAZ1}} + k_{2,\mathrm{PPVp}} \cdot c_{\mathrm{PPVp}}$$

$$\frac{dc_{\mathrm{JAZ1}}}{dt} = -k_{1,\mathrm{PPVp}} \cdot c_{\mathrm{COR}} \cdot c_{\mathrm{COI1}} \cdot c_{\mathrm{JAZ1}} + k_{2,\mathrm{PPVp}} \cdot c_{\mathrm{PPVp}}$$

c. the GFP signal

In this part, the TEV and PPV proteases cleave their target linkers to release peptides C, C′, D, and E. The highest-affinity pair C and C′, fused to GFP10 and GFP11, subsequently drives the self-assembly of the GFP10-11 fragment, which then combines with the resident GFP1-9 to reconstitute functional GFP and emit fluorescence. The ODEs of this part can be described as follows.

Reactions

$$\text{linkerCD} \xrightarrow{\text{TEVp}} \text{C} + \text{D}$$

$$\text{linkerC'E} \xrightarrow{\text{PPVp}} \text{C'} + \text{E}$$

$$\text{C} + \text{C'} \underset{k_{1,\mathrm{CC'}}}{\overset{k_{2,\mathrm{CC'}}}{\rightleftharpoons}} \text{CC'}$$

$$\text{C} + \text{D} \underset{k_{1,\mathrm{CD}}}{\overset{k_{2,\mathrm{CD}}}{\rightleftharpoons}} \text{CD}$$

$$\text{C} + \text{E} \underset{k_{1,\mathrm{CE}}}{\overset{k_{2,\mathrm{CE}}}{\rightleftharpoons}} \text{CE}$$

$$\text{D} + \text{C'} \underset{k_{1,\mathrm{C'D}}}{\overset{k_{2,\mathrm{C'D}}}{\rightleftharpoons}} \text{C'D}$$

$$\text{E} + \text{C'} \underset{k_{1,\mathrm{C'E}}}{\overset{k_{2,\mathrm{C'E}}}{\rightleftharpoons}} \text{C'E}$$

$$\text{D} + \text{E} \underset{k_{1,\mathrm{DE}}}{\overset{k_{2,\mathrm{DE}}}{\rightleftharpoons}} \text{DE}$$

$$\text{GFP1-9} + \text{CC'} \underset{k_{1,\mathrm{GFP}}}{\overset{k_{2,\mathrm{GFP}}}{\rightleftharpoons}} \text{GFP}$$

Assumptions

1. The kinetic constants (kcat and Km) of TEV protease for cutting linkerAB and linkerCD are considered the same. (Connected with part a)
2. The activities of PPV protease and TEV protease do not affect each other.
3. The substrate of the reactions catalyzed by TEV protease and PPV protease is sufficient to satisfy the Michaelis-Menten equation.
4. The binding rate constants of all peptide binding reactions are assumed to be equal, all being 10⁵ M^-1 s^-1.
5. Substances that TEV protease or PPV protease connect do not affect their function (kinetic constant). At the same time, the binding strength of the peptides linked to TEV protease and PPV protease does not take into account the linked part of TEV protease/PPV protease either.
6. When TEVp cuts both linkerAB and linkerCD simultaneously, the Michaelis-Menten equations of the two parts do not affect each other.
7. The lysis or inactivation of PPV protease and TEV protease in the system is not taken into account.

Variables

Parameters

Odes

$$\frac{dc_{\mathrm{linkerCD}}}{dt} = -\frac{k_{\mathrm{cat,TEVp}} \cdot c_{\mathrm{TEVp}} \cdot c_{\mathrm{linkerCD}}}{K_{\mathrm{m,TEVp}} + c_{\mathrm{linkerCD}}}$$

$$\frac{dc_{\mathrm{linkerC'E}}}{dt} = -\frac{k_{\mathrm{cat,PPVp}} \cdot c_{\mathrm{PPVp}} \cdot c_{\mathrm{linkerC'E}}}{K_{\mathrm{m,PPVp}} + c_{\mathrm{linkerC'E}}}$$

$$\frac{dc_{\mathrm{CC'}}}{dt} = -k_{1,\mathrm{CC'}} \cdot c_{\mathrm{CC'}} + k_{2,\mathrm{CC'}} \cdot c_C \cdot c_{C'}$$

$$\frac{dc_{\mathrm{CD}}}{dt} = \frac{k_{\mathrm{cat,TEVp}} \cdot c_{\mathrm{TEVp}} \cdot c_{\mathrm{linkerCD}}}{K_{\mathrm{m,TEVp}} + c_{\mathrm{linkerCD}}} - k_{1,\mathrm{CD}} \cdot c_{\mathrm{CD}} + k_{2,\mathrm{CD}} \cdot c_C \cdot c_D$$

$$\frac{dc_{\mathrm{CE}}}{dt} = -k_{1,\mathrm{CE}} \cdot c_{\mathrm{CE}} + k_{2,\mathrm{CE}} \cdot c_C \cdot c_E$$

$$\frac{dc_{\mathrm{C'D}}}{dt} = -k_{1,\mathrm{C'D}} \cdot c_{\mathrm{C'D}} + k_{2,\mathrm{C'D}} \cdot c_D \cdot c_{C'}$$

$$\frac{dc_{\mathrm{C'E}}}{dt} = -k_{1,\mathrm{C'E}} \cdot c_{\mathrm{C'E}} + k_{2,\mathrm{C'E}} \cdot c_E \cdot c_{C'}$$

$$\frac{dc_{\mathrm{DE}}}{dt} = -k_{1,\mathrm{DE}} \cdot c_{\mathrm{DE}} + k_{2,\mathrm{DE}} \cdot c_D \cdot c_E$$

$$\frac{dc_C}{dt} = k_{1,\mathrm{CC'}} \cdot c_{\mathrm{CC'}} - k_{2,\mathrm{CC'}} \cdot c_C \cdot c_{C'} + k_{1,\mathrm{CD}} \cdot c_{\mathrm{CD}} - k_{2,\mathrm{CD}} \cdot c_C \cdot c_D + k_{1,\mathrm{CE}} \cdot c_{\mathrm{CE}} - k_{2,\mathrm{CE}} \cdot c_C \cdot c_E$$

$$\frac{dc_D}{dt} = k_{1,\mathrm{CD}} \cdot c_{\mathrm{CD}} - k_{2,\mathrm{CD}} \cdot c_C \cdot c_D + k_{1,\mathrm{C'D}} \cdot c_{\mathrm{C'D}} - k_{2,\mathrm{C'D}} \cdot c_D \cdot c_{C'} + k_{1,\mathrm{DE}} \cdot c_{\mathrm{DE}} - k_{2,\mathrm{DE}} \cdot c_D \cdot c_E$$

$$\frac{dc_{C'}}{dt} = k_{1,\mathrm{CC'}} \cdot c_{\mathrm{CC'}} - k_{2,\mathrm{CC'}} \cdot c_C \cdot c_{C'} + k_{1,\mathrm{C'E}} \cdot c_{\mathrm{C'E}} - k_{2,\mathrm{C'E}} \cdot c_E \cdot c_{C'} + k_{1,\mathrm{C'D}} \cdot c_{\mathrm{C'D}} - k_{2,\mathrm{C'D}} \cdot c_D \cdot c_{C'}$$

$$\frac{dc_E}{dt} = k_{1,\mathrm{DE}} \cdot c_{\mathrm{DE}} - k_{2,\mathrm{DE}} \cdot c_D \cdot c_E + k_{1,\mathrm{CE}} \cdot c_{\mathrm{CE}} - k_{2,\mathrm{CE}} \cdot c_C \cdot c_E + k_{1,\mathrm{C'E}} \cdot c_{\mathrm{C'E}} - k_{2,\mathrm{C'E}} \cdot c_E \cdot c_{C'}$$

$$\frac{dc_{\mathrm{GFP}}}{dt} = k_{2,\mathrm{GFP}} \cdot c_{\mathrm{GFP1to9}} \cdot c_{\mathrm{CC'}} - k_{1,\mathrm{GFP}} \cdot c_{\mathrm{GFP}}$$

$$\frac{dc_{\mathrm{GFP1to9}}}{dt} = -k_{2,\mathrm{GFP}} \cdot c_{\mathrm{GFP1to9}} \cdot c_{\mathrm{CC'}} + k_{1,\mathrm{GFP}} \cdot c_{\mathrm{GFP}}$$

Final Result: System Reaction Kinetics

By integrating the aforementioned system of ordinary differential equations and assigning practically relevant initial concentrations and kinetic parameters, we simulated the temporal evolution of key molecular species. The resulting concentration profiles are shown below:

Fig.4 reveals that the complete reaction system reaches equilibrium within approximately 800 seconds, with GFP fluorescence intensity peaking around 700 seconds. This relatively short response time confirms the operational efficiency of our design and supports its suitability for rapid detection applications.

Hardware Response Time Estimation

The hardware implementation employs a modified pathway in Part c: instead of linking GFP10 and GFP11 to peptides C and C′, we fuse them to the N- and C-termini of trehalase (TreA), enabling CC′ dimerization to reconstitute active TreA.

The corresponding reaction is:

$$ \text{Trehalose (Treh)} \xrightarrow{\text{CC'}\,(\mathrm{TreA}),\,\mathrm{H}_2\mathrm{O}} 2\,\text{Glucose} $$

We describe this step using Michaelis–Menten kinetics:

$$ \frac{dc_{\mathrm{Treh}}}{dt} = - \frac{k_{\mathrm{cat,TreA}} \cdot c_{\mathrm{CC'}} \cdot c_{\mathrm{Treh}}}{K_{m,\mathrm{TreA}} + c_{\mathrm{Treh}}} $$

$$ \frac{dc_{\mathrm{Glucose}}}{dt} = 2 \cdot \frac{k_{\mathrm{cat,TreA}} \cdot c_{\mathrm{CC'}} \cdot c_{\mathrm{Treh}}}{K_{m,\mathrm{TreA}} + c_{\mathrm{Treh}}} $$

Using literature-derived values k_{cat, TreA} = 720 s^-1 and K_{m, TreA} = 0.3 mM (pH 7.0, 30°C)⁴, we updated the ODE system and obtained the following dynamics:

The full system equilibrates within about 900 seconds. Notably, glucose concentration peaks rapidly around 60 seconds, attributable to the high catalytic turnover rate (k_cat) of TreA. This result provides a quantitative basis for estimating the operational timeline of the hardware.

Sensitivity Analysis: Identifying Key Parameters for Acceleration

During the simulation process, we found that Part A has a significant influence on the overall system response time. We performed a local sensitivity analysis on this module, with results visualized below:

Fig.7 displays the sensitivity of concentration variables (vertical axis) to parameter variations (horizontal axis). Areas with higher values of the corresponding color indicate strong sensitivity. Two parameters—k_{cat, TEVp} and K_m,TEVp—exhibit pronounced effects on system dynamics, suggesting that optimizing these catalytic parameters offers the most effective route to further reducing the total response time.

References

1. Raran-Kurussi, S.; Toezser, J.; Cherry, S.; Tropea, J.E.; Waugh, D.S. The catalytic and other residues essential for the activity of the midgut trehalase from Spodoptera frugiperda (2010), Insect Biochem. Mol. Biol., 40, 733-741.↩
2. Wang, J., Yang, Z., Chen, C., Yao, G., Wan, X., Bao, S., ... & Jiang, H. (2024). MPEK: a multitask deep learning framework based on pretrained language models for enzymatic reaction kinetic parameters prediction. Briefings in Bioinformatics, 25(5), bbae387. ↩
3. Xue, L. C., Rodrigues, J. P., Kastritis, P. L., Bonvin, A. M., & Vangone, A. (2016). PRODIGY: a web server for predicting the binding affinity of protein–protein complexes. Bioinformatics, 32(23), 3676-3678. ↩
4. Silva, M.C.; Terra, W.R.; Ferreira, C. Differential temperature dependence of tobacco etch virus and rhinovirus 3C proteases (2013), Anal. Biochem., 436, 142-144.↩