Model

Introduction

Pesticide contamination poses a significant threat to public health and the environment, demanding safe and innovative detection methods. As a result, the PestiGuard biosensor is designed to meet this need by utilizing a molecular interaction: the binding of a pesticide to an RNA aptamer triggers a conformational change that silences a GFP reporter gene, resulting in the loss of fluorescence.

While SELEX is a traditional method for producing aptamers with high specificity and affinity, the process is costly, time-consuming, and impractical for a secondary school laboratory setting [1]. Therefore, to overcome these limitations, we designed a genetic engineering shortcut by inserting the DNA sequence of the aptamer into an E. coli plasmid, allowing the cell's machinery to produce the RNA aptamer.

We realised that one of the fundamental challenges in our design is that the aptamer sequences we obtain from research papers are typically selected as DNA aptamers, rather than RNA aptamers. The sequence of this DNA is inserted into our plasmid, but the functional ligand-binding unit is the transcribed mRNA strand. Given the structural differences between DNA and RNA, the RNA version of the DNA-derived sequence may fold into a structure with a different binding affinity. Computational folding and dynamics simulations are crucial for evaluating whether the RNA aptamer can bind to the pesticides and eventually lead to a conformation change.

To ensure our biosensor would function as intended, we conducted a three-part computational study: (1) structure generation, (2) molecular docking, and (3) molecular dynamics simulations of the aptamer-ligand complex. This approach enables us to gain deeper insights and also allows us to address two critical questions before conducting wet-lab experiments:

  1. Structural and Conformational Validation: A core design risk was whether the RNA aptamer, despite the DNA-to-RNA transition, would fold into a functional structure and maintain a stable binding interaction with the pesticide. Molecular dynamics (MD) simulations enabled us to validate our design premise before moving to the lab.
  2. Mechanistic Validation: Visually and quantitatively observe the dynamic conformational change the RNA molecule undergoes immediately after the pesticide binds. This step directly validates the entire premise of the biosensor—that binding does cause the necessary change to the aptamer itself, which may lead to silencing of the GFP reporter gene.

Applying MD to RNA aptamer-pesticide complexes is not a very common practice, especially when compared to the highly mature field of protein-ligand modeling. The primary impediment stems from a fundamental lack of established computational and structural resources: insufficient high-resolution experimental data for RNA structures, and limited comprehensive databases [2]. RNA's inherent structural flexibility further complicates accurate simulation. Consequently, our objective is to determine the feasibility and reliability of this computational approach within an under-explored chemical system.

PestiGuard Computational Design Workflow

The following workflow was implemented to validate the RNA aptamer's function and binding mechanism, moving from structure prediction to dynamic simulation.

Tools We Utilized and Their Significances

Phase 1: Structure Generation and Visualization

The initial phase converts the aptamer's linear sequence into a functional three-dimensional model, which is essential for predicting binding behavior.

1. mFold: From RNA sequence to 2D structure

mFold is the crucial starting point of our design, as it predicts the secondary structure of the RNA sequence by calculating the minimum free energy. This initial modeling enables us to assess the structure's stability and ensures that our subsequent 3D modeling and Molecular Dynamics simulations are based on a biologically plausible architecture.

2. RNAComposer: From RNA sequence to 3D structures

To construct the three-dimensional models of our RNA aptamers, we turned to RNAComposer. RNAComposer enables the generation of (.pdb) files for RNA aptamers. These structures serve as the starting models for the coming molecular docking in CHARMM-GUI and molecular dynamics simulations in GROMACS. This is the critical first step, as the unique 3D fold of each aptamer determines its ability to bind the pesticide and undergo the necessary conformational shift.

3. PyMOL: Visualization and Analysis

We utilize PyMOL, a molecular visualization tool, to visualize the three-dimensional structural files (.pdb) of RNA generated by RNAComposer. PyMOL enables us to create clear visual representations that effectively convey our design rationale.

Phase 2: Molecular Docking

This phase predicts the initial binding event and provides an empirical score for aptamer selection.

4. CHARMM-GUI: Ligand Docker

Using CHARMM-GUI's Ligand Docker function is a crucial step in predicting the interaction between the pesticide and RNA aptamer. We provide the software with the (.pdb) file of our aptamer and the (.sdf) file of the target pesticide downloaded from PubChem. The program then searches for the possible conformations to find the best-fit binding pose, scoring each one based on the estimated binding energy. This process is important as it allows us to rank the aptamers by predicted docking scores and identify the specific location and geometry of the initial binding.

Phase 3: Molecular Dynamics Simulations and Analysis of Aptamer-Ligand Complex

This is the final, most rigorous validation step, which models the dynamic behavior of the entire system over time.

5. CHARMM-GUI: MD Setup

CHARMM-GUI's second function serves as a crucial tool for preparing systems for Molecular Dynamics (MD) simulations. It constructs the simulation environment. This involves surrounding the structure with a solvated box and adding ions to neutralize the overall charge. We need to select an appropriate force field at this stage to model the molecular interactions. By performing this setup, CHARMM-GUI ensures the dynamic simulation is performed under physically realistic conditions.

6. GROMACS: Simulating the Functional Mechanism

Finally, GROMACS is used to perform molecular dynamics (MD) simulations. After CHARMM-GUI sets up the basis of the system, we utilize GROMACS to carry out energy minimization, equilibration (NVT and NPT), before running the MD, where the physical movements of the solvated system over a period of nanoseconds are stimulated. This is the final validation step of our dry lab work. By running MD on the aptamer alone and then on the aptamer-pesticide complex, we can directly observe the dynamic conformational change triggered by binding. This comparison conclusively demonstrates that pesticide binding leads to a sustained structural change, thereby proving the functional principle of our biosensor.

Stage 1: Testing on glyphosate aptamers

We began by testing glyphosate aptamers, using two different sequences. The sequence we chose for our plasmid construction, based on our computational "dry lab" work, is designated Glyphosate Y, and the one we didn't pick is Glyphosate X.

Phase 1: Structure Generation and Visualization

Constructing the 3D structures with mFold and RNAComposer

We initially modeled the RNA sequence using RNAComposer alone to generate the three-dimensional (3D) structure directly from the sequence. However, comparative docking studies revealed that these models show suboptimal stability and reduced binding affinity for the pesticide ligand. This finding motivated a revised methodology: predicting the structure's two-dimensional secondary folding pattern using mFold to create an energy-minimized blueprint, which then guides RNAcomposer's 3D modeling. The resulting mFold-guided structures demonstrated significantly improved stability, resulting in a stronger and more favorable binding interaction with the pesticide in the docking simulation. This confirmed that incorporating the thermodynamically predicted secondary structure is crucial for generating a biologically relevant 3D conformation, as the correct folding pattern is necessary to form the specific binding pocket required for ligand recognition.

We used mFold to predict the secondary structures of RNA aptamers. For each sequence, the structure exhibiting the lowest ΔG value was selected as the most thermodynamically favored conformation.

Predicted secondary structures generated by mFold

Figure 1a. 2D structure of Glyphosate Y RNA aptamer generated by mFold
ΔG = -29.28 kcal/mol
Figure 1b. 2D structure of Glyphosate X RNA aptamer generated by mFold
ΔG = -11.9 kcal/mol

Both predicted aptamer conformations demonstrate thermodynamic stability, as confirmed by their universally negative ΔG values. Furthermore, the RNA structures are characterized by the presence of stem-and-loop motifs. These exposed, single-stranded regions are structurally significant, as they frequently constitute the potential binding pockets for small-molecule ligands [3]. However, the definitive evaluation of binding competence cannot rely solely on secondary structure prediction. The binding ability and orientation of the aptamers relative to the pesticide ligand require rigorous validation through subsequent advanced computational and experimental methods, including molecular docking, molecular dynamics (MD) simulations, and wet-lab experiments. This stable architecture provides the essential framework for RNAComposer to accurately generate plausible and functional 3D models, ensuring that our subsequent docking simulations for binding affinity prediction focus only on the most viable molecular scaffolds.

Phase 2+3: Molecular Docking & Molecular Dynamics

What is Molecular Docking

Molecular docking is a key computational technique we used to predict precisely how our pesticide interacts with and binds to our RNA aptamer. Essentially, we use algorithms to find the most stable, energetically favorable 3D position where the pesticide will sit within the aptamer's binding pocket. This simulation provides a crucial output: a predicted binding score, which indicates how tightly these two molecules are expected to bind.

What is Molecular Dynamics

Molecular Dynamics (MD) is a computational simulation technique that models the movement of atoms and molecules over time by applying Newton's laws of motion. By simulating interatomic forces with a force field and integrating the equations of motion, MD provides a dynamic view of a system, offering insights into protein folding, ligand binding, and other time-dependent structural and functional processes.

We integrated molecular dynamics into our dry lab workflow because the performance of our biosensor highly depends on the aptamer's ability to bind the pesticide and undergo conformational changes. Therefore, we used MD simulations, as they are the only computational tool that can capture molecular behavior over time. We utilized MD for two purposes:

  1. Finding out the stability of the structure
    Molecular Dynamics simulations assess structural stability by generating a trajectory, while the Root Mean Square Deviation (RMSD) is the key metric used to analyze this trajectory quantitatively. To determine stability, RMSD values are plotted against time. A low, flattened, or plateaued RMSD trace after initial equilibration indicates that the structure is stable.
  2. Observing the Dynamic Mechanistic Switch
    The core principle of our biosensor is that binding triggers a conformational change that silences GFP translation. MD is essential for visualizing this dynamic event:
    • Record Molecular Interactions: We capture the continuous motion and record specific and time-dependent interactions between the pesticide and the aptamer.
    • Observe Conformational Change: MD directly models the structural change that the aptamer undergoes upon binding to the pesticide. This provides mechanistic validation that the RNA aptamer functions properly as a molecular switch.

Workflow

The 3D structure was transferred to CHARMM-GUI for system preparation. We selected the AMBER force field due to its high suitability and widespread use in simulating RNA and small molecules [4]. To align the modeling conditions with the wet lab's incubation temperature of 30°C, the temperature setting was set to 303.15 K. The alignment cross-check between our computational predictions and experimental results. The workflow proceeded in three stages:

  1. Running MD with aptamer only: We first performed an MD simulation on the aptamer alone using GROMACS. We then analyzed the trajectory to extract the lowest energy conformation observed, which served as the highly reliable receptor structure.
  2. Molecular Docking: Using this stabilized apo-aptamer conformation, we performed Molecular Docking via CHARMM-GUI to predict the initial binding pose of the pesticide within the pocket. Thirty different docked complexes, along with their estimated scores, were initially generated using the Vina, Smina, and Rxdock programs. Given Rxdock's specific optimization for RNA systems, the top-scoring complex from Rxdock (rank 1) was chosen as the representative model [5] [6]. (Please refer to Appendix 1 for our MD script)
  3. Running MD with aptamer-ligand complex: We ran MD simulations on the aptamer-ligand complexes to generate metrics that allow us to assess their binding interaction.

Results of the MD runs of the Glyphosate Aptamer only

We performed MD runs on the aptamer only and obtained RMSD graphs to analyze the stability of the aptamers.

What is RMSD

Root Mean Square Deviation (RMSD) is the gold-standard metric in structural biology and molecular docking used to quantify the difference between two molecular structures. Specifically, RMSD measures the average distance between the corresponding atoms of a predicted structure and a reference structure. A low RMSD indicates that the structure is stable and similar to the reference. High RMSD indicates a significant deviation from the reference. When plotted over time, the RMSD curve helps determine if a simulation has reached a stable state. Once the curve flattens out, it suggests the structure is fluctuating around an average equilibrium conformation.

We utilized the Root Mean Square Deviation (RMSD) analysis to track the conformational stability of the aptamer-only structure. RMSD quantifies the average atomic displacement of the RNA backbone over time relative to the initial structure, and we analyzed this graph for two critical reasons.

  1. First, for Quality Control, the primary goal was to confirm the system's thermodynamic convergence, as evidenced by the RMSD plot reaching a sustained plateau within the 50 ns window, which indicates that the system has achieved a reliable dynamic state.
  2. Second, for Stabilization, the analysis allows us to identify the equilibrium period and extract the most representative conformation. This step is crucial, as it ensures the structure used for docking is a refined structure that reflects its native flexibility in solution.

We can observe from the RMSD graphs (Figures 2a to 2b) that the plots reached a plateau after 15 nanoseconds. This plateau confirms that the complex had settled into a stable, energy-minimized binding conformation and was no longer undergoing large-scale structural drift. Critically, by then extracting the final structure for analysis at the 50 ns—well into this established plateau—we ensured that the conformation we used to analyze the binding and the resulting conformational change was the most reliable, stable, and functionally representative state of the complex, providing high confidence in our entire computational premise.

RMSD graphs of Aptamer only MD runs

Figure 2a. RMSD graph of Glyphosate Y aptamer Figure 2b. RMSD graph of Glyphosate X aptamer

Molecular Docking Analysis

Following the Aptamer-only MD run, we extracted the most stable state of the aptamer for docking. Our molecular docking analysis was conducted in three stages to assess the quality of our RNA aptamer structures: first, using a structure generated solely by RNAcomposer; second, using a structure refined by mFold before being processed by RNAcomposer; and third, using the most stable conformation extracted from the preliminary Molecular Dynamics (MD) simulation of the aptamer alone (which also incorporated mFold refinement).

The variation in docking scores is a direct result of the different mathematical scoring functions used by AutoDock Vina, Smina, and RxDock. AutoDock Vina and Smina employ fast, empirical scoring functions often trained via machine learning to correlate interaction terms with experimental binding data, and they are commonly used for protein-ligand docking [7] [8]. In contrast, RxDock utilizes a scoring function rooted in a force-field approach, incorporating more detailed physics-based terms, leading to quantitative differences in absolute scores [9]. For aptamer-ligand complexes involving highly polar molecules, RxDock is preferred because its force-field hybrid approach incorporates more detailed, explicit physical terms that better account for these complex polar and electrostatic forces. This justified our decision to rely on the RxDock scores to predict the binding stability ranking for subsequent MD simulations.

Based on the RxDock docking scores (Graph 1), we predicted that Glyphosate Y would exhibit higher binding stability in our subsequent MD runs. We acknowledge that the RNA-pesticide complex data are sparse, making docking scores merely rough predictions. Therefore, we must confirm the binding ability of aptamers using molecular dynamics simulations and, most importantly, validate the entire process with real-world wet-lab data.

Graph 1: Comparison of the docking scores of the Aptamer-Ligand Complexes across different platforms

Data analysis of the MD runs of aptamer and aptamer-ligand complexes

We ran MD simulations on the aptamer-ligand complexes to generate metrics—including Radius of gyration, SASA, and aptamer-ligand distance—that allow us to assess their binding interaction.

Radius of gyration (gmx gyrate)

The radius of gyration is the imaginary distance from the centroid at which the area of cross-section is imagined to be focused at a point to obtain the same moment of inertia. As its name suggests, this command computes the radius of gyration of a molecule and the radii of gyration about the x-, y-, and z-axes, as a function of time. The atoms are explicitly mass weighted. The axis components correspond to the mass-weighted root-mean-square of the radii components orthogonal to each axis, for example:

Where w_i is the weight value in the given situation (mass, charge, unit)

We utilise the radius of gyration to determine whether the aptamer is stably folded. The less fluctuation shown in the graph, the less change in shape the aptamer in the complex is. A drastic change in the radius of gyration indicates a conformational change in the aptamer. If the change in the radius of gyration is negligible, it means that the shape and center of mass of the aptamer have not changed substantially, remaining the same shape as the original shape.

Not only should we examine its stability, but its level is also crucial for our understanding of the state of the aptamer-ligand complex. A large radius of gyration indicates that the masses are spatially distributed apart. This is mainly due to reasons such as the aptamer not successfully making a conformational change and cannot bind with the ligand tightly. Therefore, we would like the radius of gyration to be as small as possible, meaning that the aptamer has a high chance of successfully binding with the ligand.

Solvent accessible surface area (gmx sasa)

As the name suggests, the Solvent Accessible Surface Area (SASA) measures the surface area of a biomolecule, in this case, the aptamer and the ligand complex, that is accessible by the solvent molecules. It is defined as the extent to which atoms on the surface of a molecule can form contact with the solvent, and it is generally measured in squared nanometers (nm²). We measured the Solvent-Accessible Surface Area (SASA) of the aptamer-ligand complex throughout the MD run trajectory. A drop in the total SASA of the aptamer-ligand complex over the course of an MD simulation indicates that a significant portion of the ligand is effectively interacting with the aptamer's surface.

Aptamer-ligand distance graph

To examine the effectiveness of the binding between the aptamer and the ligand, we generated an aptamer-ligand distance graph to visualize the interaction dynamics. This process began with the extraction of all frames from the molecular dynamics (MD) trajectory produced during the simulation run, which captures the conformational changes of the aptamer-ligand complex over time. To accurately calculate the distance between the center of the aptamer and that of the ligand, we modified a pre-existing Bash script found in a GROMACS tutorial written by Justin A. Lemkul, Ph.D. Ensuring that it suited our specific molecular system. This script enabled us to systematically compile a .dat file containing the distances for each frame, thereby transforming the raw trajectory data into a format suitable for analysis. After generating the .dat file, we utilized Xmgrace, a powerful plotting tool, to create a graph that visually represents the distance fluctuations throughout the simulation. In this graph, each point corresponds to a specific frame, enabling us to track the distance between the aptamer and ligand over time. We labeled the axes, indicating the distance in nanometers and frames in the MD run. Through this analysis, we can gain a deeper understanding of the binding affinity between the aptamer and the ligand.

Results of the MD runs of Glyphosate Aptamer and the Aptamer-Ligand Complex

Glyphosate Y

Figure 3a. 3D structure of the Glyphosate Y docked-complex

Video 1. The simulation of the Glyphosate Y docked-complex (25ns)

Figure 3b. RMSD graph of Glyphosate Y aptamer

Figure 3c. Radius of gyration graph of Glyphosate Y aptamer

Both the radius of gyration and Root Mean Square Deviation (RMSD) graphs confirm the high structural stability of the Glyphosate Y aptamer. The minimal fluctuations observed in the radius of gyration graph indicates that the aptamer maintains a consistent shape, while the stable profile shown in the RMSD graph further reinforces that the molecule undergoes very few conformational changes.

Figure 3d. Radius of gyration graph of Glyphosate Y docked complex

Figure 3e. SASA graph of Glyphosate Y docked complex

Figure 3f. Aptamer-ligand distance graph of Glyphosate Y docked complex

The glyphosate Y shows success in its trajectory, successfully binding together. From the graph of the radius of gyration of Glyphosate Y docked complex, we can see that there's a low level of fluctuations(within 0.25nm), meaning that the shape of the aptamer doesn't change drastically, meaning that the aptamer has stably bound with the ligand or that the aptamer does not bind with the ligand at all. However, in the case of glyphosate Y, the aptamer has clearly bound to the ligand throughout the entire 25 ns. First, we can see that the solvent-accessible surface area remains stable at a relatively low level throughout the whole 25 ns, indicating that the ligand and the aptamer have strong interactions. Finally, this is further supported by the fact that the distance between the aptamer and the ligand remains constant at a low level, indicating that the aptamer is successfully bound to the ligand.

Glyphosate X

Figure 4a. 3D structure of the Glyphosate X docked-complex

Video 2. The simulation of the Glyphosate X docked-complex (25ns)

Figure 4b. RMSD graph of Glyphosate X aptamer

Figure 4c. Radius of gyration graph of Glyphosate X aptamer

Both the radius of gyration and Root Mean Square Deviation (RMSD) graphs confirm the high structural stability of the Glyphosate X aptamer. The minimal fluctuations observed in the radius of gyration graph indicate that the aptamer maintains a consistent shape, while the stable profile shown in the RMSD graph further reinforces that the molecule undergoes very few conformational changes.

Figure 4d. Radius of gyration graph of Glyphosate X docked complex

Figure 4e. SASA graph of Glyphosate X docked complex

Figure 4f. Aptamer-ligand distance graph of Glyphosate X docked complex

The glyphosate X shows success in its trajectory, indicating that the aptamer and the ligand have a strong interaction, specifically within the first 15 ns. From the graph of the radius of gyration of Glyphosate X docked complex, we can see that there's a low level of fluctuations(within 0.25nm), meaning that the shape of the aptamer doesn't change drastically, and that the aptamer has stably bound with the ligand or that the aptamer does not bind with the ligand at all. However, in the case of glyphosate X, the aptamer has clearly bound to the ligand for the first 15 ns. First, we can see that the solvent accessible surface area remains stably at a relatively low level(66nm^2) for the first 15ns, and then increases to 68nm^2 till 25ns. This means that the ligand and the aptamer have strong interactions for the first 15ns and less strong interactions for the remaining 10ns. Finally, this is further proved by the point that the distance between the aptamer and the ligand also remains constant at a low level for the first 15ns, meaning that the aptamer is bound with the ligand successfully for the first part.

Results from computational modeling, including mFold, docking, and molecular dynamics (MD) simulations, suggest that the designed Glyphosate Aptamers may not achieve the high binding strength of the initial DNA aptamer. However, these in silico findings confirmed a reasonable likelihood of binding to the glyphosate molecule. Given this binding potential, the team decided to move forward with experimental (wet lab) validation. Specifically, Glyphosate Aptamer Y exhibited better simulated performance compared to Aptamer X, since it can bind to the ligand stably for a longer period of time, leading to its selection for use in constructing the experimental plasmid.

Stage 2: Testing on other aptamers

After we completed testing the glyphosate aptamers (Glyphosate X and Y), we proceeded to test the other aptamers in the study.

Phase 1: Structure Generation and Visualization

We used mFold to predict the secondary structures of RNA aptamers. Again, all of the predicted aptamer conformations are confirmed to be thermodynamically stable, as indicated by their consistently negative ΔG values. Additionally, these RNA structures share a key characteristic: they all possess stem-and-loop motifs.

Predicted secondary structures generated by mFold

Figure 5a. 2D structure of Malathion RNA aptamer generated by mFold
ΔG = -5.0 kcal/mol
Figure 5b. 2D structure of Acephate RNA aptamer generated by mFold
ΔG = -18.97 kcal/mol
Figure 5c. 2D structure of Acetamiprid RNA aptamer generated by mFold
ΔG = -13 kcal/mol
Figure 5d. 2D structure of Chlorpyrifos RNA aptamer generated by mFold
ΔG = -29 kcal/mol

Phase 2+3: Molecular Docking & Molecular Dynamics

Similar to the glyphosate aptamers, we performed the MD run on the aptamer alone and extracted the lowest energy conformation observed for docking.

We can observe from the RMSD graphs (Figures 6a to 6d) that the plots reached a plateau after 15 nanoseconds. We then extracted the final structure for analysis at 50 ns.

RMSD graphs of Aptamer only MD runs

Figure 6a. RMSD graph of Acephate aptamer Figure 6b. RMSD graph of Chlorpyrifos aptamer
Figure 6c. RMSD graph of Acetamiprid aptamer Figure 6d. RMSD graph of Malathion aptamer

Molecular Docking Analysis

Following the Aptamer-only MD run, we extracted the most stable state of the aptamer for docking. The results (Graph 2) clearly show that incorporating mFold to predict the secondary structure improves the docking scores, and using the most stable conformation from the preliminary MD run yields structures with the highest predicted stability for nearly all aptamers.

Graph 2: Comparison of the RXDock docking scores between the aptamer and ligand for structures generated solely by RNAcomposer, using a structure refined by mFold before being processed by RNAcomposer, and using the most stable conformation extracted from MD simulation of the aptamer alone (which also incorporated mFold refinement)

However, the docking results indicate that these aptamers exhibit a weaker binding ability compared to the glyphosate aptamers we have tested.

Graph 3: Comparison of the docking scores of the Aptamer-Ligand Complexes across different platforms

Lastly, we performed the MD runs on the Aptamer-Ligand Complexes, and the results are as follows:

Results of the MD runs of Aptamer and the Aptamer-Ligand complex

Acephate

Figure 7a. 3D structure of the Acephate docked-complex

Video 3. The simulation of the Acephate docked-complex

Figure 7b. RMSD graph of Acephate aptamer

Figure 7c. Radius of gyration graph of Acephate aptamer

Both the radius of gyration and Root Mean Square Deviation (RMSD) graphs confirm the high structural stability of the Acephate aptamer. The minimal fluctuations observed in the radius of gyration graph indicates that the aptamer maintains a consistent shape, while the stable profile shown in the RMSD graph further reinforces that the molecule undergoes very few conformational changes.

Figure 7d. Radius of gyration graph of Acephate docked complex

Figure 7e. SASA graph of Acephate docked complex

Figure 7f. Aptamer-ligand distance graph of Acephate docked complex

According to the above graphs, we can determine that the binding between acephate and the ligand is not particularly successful and can only occur within the first 5,000 ps. First, let's take a look at the graph of the radius of gyration. The radius fluctuates vigorously, meaning that its shape keeps changing and cannot stabilise when bound to the ligand. Next, we can see that the SASA graph actually decreases over time, meaning that the contact area between the aptamer and ligand has decreased or the aptamer is folding on its own. However, the outcome should be the latter. As the distance between the aptamer and the ligand fluctuates significantly, meaning their distance changes constantly, and they cannot bind together for a prolonged period.

The results mentioned above may be attributed to the transition from DNA to RNA, which renders the original DNA-binding pocket inaccessible, resulting in the poor outcomes observed.

Chlorpyrifos

Figure 8a. 3D structure of the Chlorpyrifos docked-complex

Video 4. The simulation of the Chlorpyrifos docked-complex

Figure 8b. RMSD graph of Chlorpyrifos aptamer

Figure 8c. Radius of gyration graph of Chlorpyrifos aptamer

Both the radius of gyration and Root Mean Square Deviation (RMSD) graphs confirm the low structural stability of the Chlorpyrifos aptamer. The rigorous fluctuations observed in the radius of gyration graph indicate that the aptamer changes its shape constantly, while the unstable profile shown in the RMSD graph further reinforces that the molecule undergoes lots of conformational changes.

Figure 8d. Radius of gyration graph of Chlorpyrifos docked complex

Figure 8e. SASA graph of Chloropyrifos docked complex

Figure 8f. Aptamer-ligand distance graph of Chlorpyrifos docked complex

Based on the above graphs, we can infer that the binding of the aptamer and ligand is not very ideal. We can deduce the aptamer has bound to the ligand for the first 10 ns since the distance between the aptamer and the ligand remains constantly at a relatively low level for the first 10ns. Then for the remaining 15ns, the distance between the aptamer and the ligand increases and the solvent accessible surface area decreases. This is most likely due to the ligand's unbinding resulting from the aptamer's conformational change.

Malathion

Figure 9a. 3D structure of the Malathion docked-complex

Video 5. The simulation of the Malathion docked-complex

Figure 9b. RMSD graph of Malathion aptamer

Figure 9c. Radius of gyration graph of Malathion aptamer

Both the radius of gyration and Root Mean Square Deviation (RMSD) graphs confirm the low structural stability of the Malathion aptamer. The rigorous fluctuations observed in the radius of gyration graph indicates that the aptamer changes its shape constantly, while the unstable profile shown in the RMSD graph further reinforces that the molecule undergoes lots of conformational changes.

Figure 9d. Radius of gyration graph of Malathion docked complex

Figure 9e. SASA graph of Malathion docked complex

Figure 9f. Aptamer-ligand distance graph of Malathion docked complex

Based on the above graphs, we can infer that the binding of the aptamer and ligand is not very ideal. We can deduce the aptamer has bound to the ligand for the first 10 ns since the distance between the aptamer and the ligand remains constantly at a relatively low level for the first 10ns. Then for the remaining 15ns, the distance between the aptamer and the ligand increases and the solvent accessible surface area decreases. This is most likely due to the ligand's unbinding resulting from the aptamer's conformational change.

Acetamiprid

Figure 10a. 3D structure of the Acetamiprid docked-complex

Video 6. The simulation of the Acetamiprid docked-complex

Figure 10b. RMSD graph of Acetamiprid aptamer

Figure 10c. Radius of gyration graph of Acetamiprid aptamer

Both the radius of gyration and Root Mean Square Deviation (RMSD) graphs confirm the high structural stability of the Acetamiprid aptamer. The minimal fluctuations observed in the radius of gyration graph indicates that the aptamer maintains a consistent shape, while the stable profile shown in the RMSD graph further reinforces that the molecule undergoes very few conformational changes.

Figure 10d. Radius of gyration graph of Acetamiprid docked complex

Figure 10e. SASA graph of Acetamiprid docked complex

Figure 10f. Aptamer-ligand distance graph of Acetamiprid docked complex

Based on the above graphs, we can infer that the binding of the aptamer and ligand is not ideal. The radius of gyration is fluctuating rigorously for the whole 25ns, meaning that the shape of the aptamer keeps on changing, therefore, the aptamer cannot form a stable complex with the ligand. The solvent accessible surface area and the distance between the aptamer and the ligand are also fluctuating rigorously, meaning that the ligand cannot consistently bind with the aptamer.

Modeling the RGB Lambda Law Using eGFP Reporter Gene with aptamer in Response to Varying Pesticide Concentrations

Abstract

This model investigates the RGB Lambda Law in the expression of a eGFP reporter gene with aptamer against to the pesticide concentrations.[1] We assess various concentrations (0, 0.3, 1, 3, 10, 30, 100, 300, 1000 mg/L) to determine the range where a linear relationship exists. Linear regression is used to establish a standard curve and evaluate the model with R-squared values.

Introduction

The RGB Lambda Law provides a framework for quantifying the fluorescence emitted by a eGFP reporter gene with aptamer in response to varying pesticide concentrations. This study aims to identify the concentration range that maintains a linear relationship with GFP expression, facilitating the development of a biosensor for pesticide detection.

Background

  • eGFP Reporter Gene: We transformed E.coli BL21 with a plasmid containing promoter, lac operator, pesticide-binding aptamer, reporter gene EGFP and T7 terminator. When a pesticide binds to the aptamer that is transcribed onto the mRNA, the pesticide triggers a conformational change that blocks the ribosome binding site and thus blocks the reporter gene translation. The reporter proteins therefore cannot be produced as effectively and the bacteria will glow dimmer. In the absence of pesticides, the RNA aptamer adopts a default conformation that leaves the RBS accessible. The ribosome can bind, initiate translation, and the cell produces the fluorescent protein, emitting a strong signal.
  • RGB Color Model: The RGB model blends red, green, and blue light to represent colors visible to the human eye. The exact RGB values can be calculated using spectroscopic data, where the intensity of the emitted light at various wavelengths is mapped to the RGB channels.[2] The digital images obtained from our device are analyzed. The intensity of green fluorescence can be quantified pixel by pixel to measure eGFP expression levels in cells. Describe how the RGB model applies to the fluorescence emitted from eGFP and its quantitative measurement.
  • Standard Curve: To ensure accuracy, the concept of a standard curve in quantitative analysis is introduced and the significance of linear regression in establishing a relationship between GFP fluorescence and pesticide concentration.[3]

Model Description

Mathematical Representation

The model can be expressed as:

RGB = R + G + B

Where:

  • RGB is the resulting color intensity of the emitted fluorescence.
  • R, G, and B represent the intensities of red, green, and blue light, respectively.

Variables

  • R: Intensity of red light
  • G: Intensity of green light
  • B: Intensity of blue light
  • RGB: Resulting color output from GFP fluorescence
  • [P]: Concentration of pesticide (0, 0.3, 1, 3, 10, 30, 100, 300, 1000 mg/L)

Implementation

Synthetic Biology Applications

By capturing a photo of the GFP produced from our system and uploading on the software, it offers comprehensive colour analysis to correlate brightness to pesticide concentration. The software also stores previous results for easy tracking and allows users to delete data quickly when needed.

Experimental Setup

  1. Inoculate a single colony of each transformed BL21(DE3) into 5 mL of LB medium containing kanamycin.
  2. Incubate the starter cultures overnight at 37°C with shaking at 220 RPM.
  3. Dilute 0.5 mL of each overnight starter culture into 50 mL of fresh, kanamycin-supplemented LB medium.
  4. Incubate the cultures at 37°C with shaking at 220 RPM.
  5. Monitor the OD600 every 30 minutes. When the OD600 reaches 0.5–0.6, add IPTG to reach suitable concentration and proceed to the next step.
  6. Aseptically split each culture into 9 separate tubes, each containing 5 mL culture, using 15 mL conical centrifuge tubes.
  7. Label the tubes from 1 to 9. Add corresponding pesticide to each tube according to the concentrations specified in Table 1 below.
Table 1: The concentration of pesticide in each set-up

  1. Induce protein expression by incubating all tubes at 30°C with shaking at 220 RPM for a suitable duration.
  2. Transfer 1 mL from each tube to a microcentrifuge tube and measure OD600 of each tube.
  3. Analyze the results using the designated hardware (Smartbox) and software (PestiGuard biosensor platform).
  4. Repeat step 1 -- step 10 two times to get total three replicates for each pesticide.

Generating a Standard Curve

  • Linear Regression: This is a statistical method used to model the relationship between two variables by fitting a linear equation to observed data. [4] When studying GFP fluorescence (in terms of RGB) as a function of pesticide concentration, linear regression can be employed to determine how fluorescence intensity (dependent variable) changes with varying pesticide concentrations (independent variable). The goal is to find the best-fit line:

    RGB = m[P] + b,

    Where [P] is the known pesticide concentrations (in mg/L).

    RGB is the measured GFP fluorescence intensities (intensity values from an RGB analysis).

    m is the slope of the line (rate of change of fluorescence with respect to concentration),

    b is the y-intercept (fluorescence intensity value when pesticide concentration is zero).

  • R-Squared Value: This is calculated to measure the goodness of fit. It represents the proportion of variance in fluorescence explained by the pesticide concentration. Once the regression line is established, the equation can be used to predict pesticide concentrations from new GFP fluorescence measurements. For example, given a fluorescence intensity value, the corresponding pesticide concentration can be calculated by rearranging the regression equation:

    [P] = (RGB - b) / m.

Identifying the concentration range where a linear relationship exists between GFP fluorescence and pesticide concentration is critical in ensuring the accuracy, reliability, and interpretability of the results. This is particularly important when focusing on lower concentrations, as these are often the most biologically or environmentally relevant. For the accuracy in quantification, linear regression assumes a direct proportionality between GFP fluorescence and pesticide concentration. If this relationship is not linear across all concentrations, using the regression equation outside the linear range can lead to inaccurate quantifications. For lower concentrations, fluorescence intensities are often more subtle and closer to background levels. Ensuring linearity in this range allows for sensitive and precise measurements of pesticide presence, which is often critical in detecting trace amounts. This also help avoiding signal saturation at higher concentrations.

Results

For glyphosate

2. Linear Range and R² Values

  • 0–1000 mg/L (Top-Left):
    • Equation: RGB = -0.2271[P] + 247.85
    • R² = 0.3324
    • Observation:
      • Weak linear relationship (low R² value).
      • Nonlinear behavior dominates, likely due to saturation or quenching at high concentrations.
      • Not suitable for quantitative analysis over this wide range.
  • 0–100 mg/L (Top-Right):
    • Equation: RGB = -2.3947[P] + 298.27
    • R² = 0.5127
    • Observation:
      • Moderate linearity, but the R² value suggests significant deviation from linear behavior.
      • Higher concentrations (e.g., near 100 mg/L) likely introduce nonlinear effects.
  • 0–10 mg/L (Bottom-Left):
    • Equation: RGB = -20.475[P] + 365.63
    • R² = 0.9461
    • Observation:
      • Strong linear relationship with high R² value.
      • This range demonstrates good sensitivity and reliable quantification.
  • 0–1 mg/L (Bottom-Right):
    • Equation: RGB = -31.245[P] + 378.76
    • R² = 0.9777
    • Observation:
      • Excellent linearity with the highest R² value.
      • Ideal range for detecting trace concentrations of glyphosate.

3. Nonlinear Behavior at Higher Concentrations

  • Nonlinearity at High Concentrations (e.g., 100 mg/L and above):
    • The graphs for 0–1000 mg/L and 0–100 mg/L show clear nonlinear behavior.
    • Causes:
      • Fluorescence saturation: At high concentrations, GFP fluorescence may reach its maximum intensity, leading to a plateau.
      • Quenching effects: High concentrations of glyphosate might interfere with GFP fluorescence by altering the chemical environment of the chromophore.
      • Instrumental limitations: Detection systems may lose sensitivity or accuracy at very high fluorescence intensities.

4. Importance of Linear Range Identification

  • The 0–10 mg/L and 0–1 mg/L ranges are identified as the most reliable for quantifying glyphosate concentrations, based on their high R² values.
  • These ranges enable:
    • Accurate calibration curves for quantification.
    • Reliable detection of trace concentrations, which is often critical in environmental or agricultural analyses.

Conclusion

  • The 0–1 mg/L range is the most suitable for sensitive and accurate glyphosate detection, with an R² of 0.9777.
  • Nonlinear behavior at higher concentrations (e.g., 100 mg/L and above) limits the utility of these ranges for quantification, likely due to fluorescence saturation or quenching.
  • Focusing on the lower concentration ranges ensures robust and reproducible measurements for practical applications.

Further work on the trial of other pesticides with comparison to glyphosate

Linear Range and R² Values

  • For 0–1000 mg/L,
    • Equation: RGB = -0.2271[P] + 247.85
    • R² = 0.3324
    • Observation:
      • Weak linear relationship (low R² value).
      • Nonlinear behavior dominates, likely due to saturation or quenching at high concentrations.
      • Not suitable for quantitative analysis over this wide range.
  • For 0–100 mg/L,
    • Equation: RGB = -2.3947[P] + 298.27
    • R² = 0.5127
    • Observation:
      • Moderate linearity, but the R² value suggests significant deviation from linear behavior.
      • Higher concentrations (e.g., near 100 mg/L) likely introduce nonlinear effects.
  • For 0–10 mg/L,
    • Equation: RGB = -20.475[P] + 365.63
    • R² = 0.9461
    • Observation:
      • Strong linear relationship with high R² value.
      • This range demonstrates good sensitivity and reliable quantification.
  • For 0–1 mg/L,
    • Equation: RGB = -31.245[P] + 378.76
    • R² = 0.9777
    • Observation:
      • Excellent linearity with the highest R² value.
      • Ideal range for detecting trace concentrations of glyphosate.

The 0–1 mg/L range is the most suitable for sensitive and accurate glyphosate detection, with an R² of 0.9777. Nonlinear behavior at higher concentrations (e.g., 100 mg/L and above) limits the utility of these ranges for quantification, likely due to fluorescence saturation or quenching. Focusing on the lower concentration ranges ensures robust and reproducible measurements for practical applications.





Summary Comparison of Standard Curves

Pesticide Conc. Range(mg/L) R² Value Observations Linear Trend
Glyphosate 0–1000 0.3324 Weak linear relationship; significant nonlinear behavior at higher concentrations. Not suitable for quantification.
0–100 mg/L 0.5127 Moderate linearity; deviations suggest nonlinearity near 100 mg/L. Fairly linear but limited.
0–10 mg/L 0.9461 Strong linear relationship; good sensitivity and reliability. Highly suitable for detection.
0–1 mg/L 0.9777 Excellent linearity; ideal for detecting trace concentrations. Best range for quantification.
Acephate 0–1000 mg/L 0.1892 Weak linearity; poor fit for quantification. Not suitable for quantification.
0–100 mg/L 0.4265 Moderate linearity; still not ideal for reliable quantification. Limited linear relationship.
0–10 mg/L 0.8374 Stronger fit than higher concentrations; more reliable. Fairly suitable for detection.
0–1 mg/L 0.967 Excellent linearity; good for trace detection. Highly suitable for quantification.
Acetamiprid 0–1000 mg/L 0.3099 Weak linearity; ineffective for quantification. Not suitable for quantification.
0–100 mg/L 0.3492 Slightly improved fit, but still low for reliable analysis. Limited linear relationship.
0–10 mg/L 0.8875 Stronger fit; better sensitivity. Fairly suitable for detection.
0–1 mg/L 0.9009 Excellent fit; capable of detecting low pesticide concentrations. Highly suitable for quantification.
Chlorpyrifos 0–1000 mg/L 0.3099 Weak linearity; not suitable for quantification. Not suitable for quantification.
0–100 mg/L 0.5296 Moderate linearity; still limited for accurate measurements. Fairly limited linearity.
0–10 mg/L 0.8687 Strong linearity; better sensitivity than higher ranges. Fairly suitable for detection.
0–1 mg/L 0.9619 Excellent fit; ideal for trace detection. Highly suitable for quantification.
Malathion 0–1000 mg/L 0.3481 Weak linearity; ineffective for quantification. Not suitable for quantification.
0–100 mg/L 0.4954 Moderate linearity; still not reliable for precise measurements. Limited linear relationship.
0–10 mg/L 0.5839 Fair fit; better than higher concentrations but still limited. Moderate suitability for detection.
0–1 mg/L 0.7916 Strong linearity; good for trace detection. Fairly suitable for quantification

Discussion

In general, all pesticides fluorescence data show a decrease in GFP fluorescence intensity (RGB values) with increasing concentrations. However, the patterns differ in terms of their linearity, sensitivity (slope), and behavior at higher concentrations due to differences in their interaction with GFP. Glyphosate performs better over a broader concentration range (0–10 mg/L) due to stronger linearity and having the highest sensitivity in range of 0-1mg/L.

Conclusion

Our modeling work confirms a central challenge in translating DNA aptamers to their corresponding RNA transcripts: the changes in 3D structure. This structural difference is expected to trigger substantial refolding in the RNA aptamer, which likely leads to the inaccessibility of the original DNA-binding pocket while simultaneously forming new binding sites. Given this, our primary computational goal was to determine whether the RNA aptamers could still fold into a functional structure that maintains a stable interaction with the targeted pesticides, and subsequently to predict and rank the binding ability of the RNA aptamers under study.

The computational results provided evidence for the potential functionality of the RNA transcripts. mFold analysis confirmed that all predicted RNA aptamer conformations are thermodynamically stable and feature the stem and loop motifs necessary to facilitate pesticide binding.

Docking estimated the binding affinity ranking (from highest to lowest) for the six targeted pesticides as follows: Acephate, Glyphosate Y, Glyphosate X, Chlorpyrifos, Acetamiprid, and Malathion, which closely matched the results observed in the wet-lab experiments.

Critically, among all the aptamers, Molecular Dynamics simulations indicated that the RNA aptamers for Glyphosate Y demonstrated the most stable binding. This suggests that despite the structural transition from DNA to RNA, the Glyphosate RNA aptamer retains or forms new binding pockets capable of interacting with its target.

The observation in the simulations also shows that the structure of the RNA aptamer changes upon binding, and this may influence the expression of the GFP gene. This computational finding is directly corroborated by our wet-lab result, which showed that glyphosate performed best experimentally. In contrast, the MD performance for Chlorpyrifos, Acephate, Malathion, and Acetamiprid suggests weaker binding. However, it is essential to note that the dry lab modeling is a prediction, and there are inherent limitations to these computational methods, such as the following:

The Limits of Computational Prediction

1. Force-field inaccuracies

Our dry lab work is built upon fundamental computational assumptions. We acknowledge that the force fields used in our MD simulations are an approximation, which means the force fields may not perfectly capture the complex electrostatic and solvent interactions. The dynamic movements are thus not absolute certainties, but rather estimates.

2. Timescale limitation

Furthermore,our MD simulations for the complex are constrained to a 25 ns runtime.It is possible that such duration is insufficient to fully capture the biologically relevant mechanism of aptamer-pesticide binding or conformational changes that may occur over the timescale using nanoseconds, such as microseconds (μs), which is far longer than our current simulation duration, thus potentially missing key events.

Simplified Biological Context

Our models treat the aptamer as an isolated entity in an aqueous solution. However, in reality, the aptamer will be functioning within the crowded cytoplasm of an E. coli cell. Such simplification of the biological environment in our model may neglect the physical presence of other cellular machinery, which could have an effect on the binding effectiveness of the aptamer to the pesticide. Such omission in our simulation for the action of the aptamer may not align with its actual performance in the real-life environment.

Structural and Functional Uncertainties

1. Folding vs. Binding-Induced Change

RNA is an inherently flexible molecule. The conformational change observed in our MD trajectories, particularly at the beginning of the simulation, may represent the aptamer's spontaneous folding into a more stable state within the simulated environment, rather than a direct, ligand-induced conformational shift. Distinguishing between intrinsic folding dynamics and binding-triggered changes remains a subtle challenge.

2. Suboptimal Binding Affinity

The RNA sequences were derived from literature-reported DNA aptamers. As RNA and DNA exhibit structural differences, the transcribed RNA version may not fold into the optimum structure required for maximum binding. Even if the MD results confirm stable binding, this virtual interaction might represent a sub-optimal affinity compared to a genuinely in vitro selected aptamer.

3. Unproven Reporter Gene Silencing

Our dry lab work focuses solely on the aptamer-ligand interaction and the resulting conformational change. We were unable to computationally model the final, crucial step: proving that this conformational change is sufficient to silence the GFP fluorescence gene. This link remains a hypothesis to be validated in the wet lab.

Suggestions for Future Dry Lab Improvements

The previously mentioned limitations indicate that our MD results provide a strong qualitative framework for validating aptamer function, but they still require wet-lab experiments to complementarily confirm the actual effectiveness of our biosensor, while further enhancing our modeling process. Below are some additional possible improvements.

Deepen Computational Precision

To strengthen the certainty of our predictions, our future dry lab work must increase its computational precision. We should move beyond basic stability checks and implement Binding Free Energy Calculations to provide a more rigorous, quantitative, and thermodynamically relevant estimate of binding affinity (ΔG).

Expand the Scope for Comprehensive Validation

Our next step is to expand the scope to model the complete biosensor system. This means building and simulating a structural model of the aptamer-RBS interaction. Validating this silencing mechanism is crucial for demonstrating that the entire design principle is effective.

Furthermore, before any costly wet-lab experiments, we must computationally assess the sensor's real utility by testing its selectivity and specificity by running comparative simulations (docking and MD) with a panel of structurally similar non-target contaminants to check the specificity of aptamers. This will allow us to predict cross-reactivity and prioritize the most selective aptamer candidates.

References

[1] Stephens, C., Goodey, N. M., & Gubler, U. (2025). A beginners guide to SELEX and DNA aptamers. Analytical Biochemistry, 703, 115890.
https://doi.org/10.1016/j.ab.2025.115890

[2] D'Esposito, R. J., Myers, C. A., Chen, A. A., & Sweta Vangaveti. (2022). Challenges with Simulating Modified RNA: Insights into Role and Reciprocity of Experimental and Computational Approaches. Genes (Basel), 13(3), 540–540. https://doi.org/10.3390/genes13030540

[3] Ansari, N., Liu, C., Hédin, F., Hénin, J., Ponder, J. W., Ren, P., Piquemal, J., Lagardère, L., & Hage, K. E. (2025). Targeting RNA with small molecules using state-of-the-art methods provides highly predictive affinities of riboswitch inhibitors. Communications Biology, 8(1).
https://doi.org/10.1038/s42003-025-08809-y

[4] Hegde, S., Akhter, S., Tang, Z., Qi, C., Yu, C., Lewicka, A., Liu, Y., Kushal Koirala, Mikhail Reibarkh, Battaile, K., Cooper, A., Lovell, S., Holmstrom, E., Wang, X., Piccirilli, J., Gao, Q., Miao, Y., & Wang, J. (2025). Mechanistic studies of small molecule ligands selective to RNA single G bulges. Nucleic Acids Research, 53(12). https://doi.org/10.1093/nar/gkaf559

[5] Dickerhoff, J., Warnecke, K. R., Wang, K., Deng, N., & Yang, D. (2021). Evaluating Molecular Docking Software for Small Molecule Binding to G-Quadruplex DNA. International Journal of Molecular Sciences, 22(19), 10801–10801.
https://doi.org/10.3390/ijms221910801

[6] Ruiz-Carmona, S., Alvarez-Garcia, D., Foloppe, N., Garmendia-Doval, A. B., Juhos, S., Schmidtke, P., Barril, X., Hubbard, R. E., & Morley, S. D. (2014). rDock: A Fast, Versatile and Open Source Program for Docking Ligands to Proteins and Nucleic Acids. PLoS Computational Biology, 10(4), e1003571. https://doi.org/10.1371/journal.pcbi.1003571

[7] Trott, O., & Olson, A. J. (2009). AutoDock Vina: Improving the Speed and Accuracy of Docking with a New Scoring function, Efficient optimization, and Multithreading. Journal of Computational Chemistry, 31(2).
https://doi.org/10.1002/jcc.21334

[8] Koes, D. R., Baumgartner, M. P., & Camacho, C. J. (2013). Lessons Learned in Empirical Scoring with smina from the CSAR 2011 Benchmarking Exercise. Journal of Chemical Information and Modeling, 53(8), 1893–1904.
https://doi.org/10.1021/ci300604z

[9] Hartshorn, M. J., Verdonk, M. L., Chessari, G., Brewerton, S. C., Mooij, W. T. M., Mortenson, P. N., & Murray, C. W. (2007). Diverse, High-Quality Test Set for the Validation of Protein−Ligand Docking Performance. Journal of Medicinal Chemistry, 50(4), 726–741.
https://doi.org/10.1021/jm061277y

[10] K. Phopin and T. Tantimongcolwat, "Pesticide Aptasensors—State of the Art and Perspectives," Sensors, vol. 20, no. 23, p. 6809, Nov. 2020, doi: 10.3390/s20236809.

[11] Daniele C. Rodman, "Spectroscopy and RGB-colorimetry for quantification of plant pigment and fruit content in fruit drinks," Lund University, 2018.

[12] "Linear Standard Curves1 - GraphPad."

[13] David Harvey, "Linear Regression and Calibration Curves," Analytical Chemistry 2.1.

Appendix

MD script for the MD run

Aptamer only MD commands: (added ions and built box and solvate via CHARMM-GUI solution builder)

Energy minimization: Assemble the binary input using gmx grompp

gmx grompp -f em.mdp -c step3_input.gro -r step3_input.gro -p topol.top -o em.tpr

gmx mdrun -v -deffnm em

Equilibration - nvt : Assemble the binary input using gmx grompp

gmx grompp -f eq.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr -n index.ndx

gmx mdrun -deffnm nvt

Equilibration - npt : Assemble the binary input using gmx grompp

gmx grompp -f eq.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr -n index.ndx

gmx mdrun -deffnm npt

Note: eq.mdp is generated from CHARMM-GUI that is used for both nvt and npt

MD simulation:

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_10.tpr -n index.ndx

gmx mdrun -deffnm md_0_10 -v

Analysis:

Correcting for Periodicity Effects:

gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_noPBC.xtc -pbc mol -center -n index.ndx

RMSD graph:

gmx rms -s md_0_10.tpr -f md_0_10_noPBC.xtc -o rmsd.xvg -tu ns -n index.ndx

Radius of gyration:

gmx gyrate -s md_0_10.tpr -f md_0_10_noPBC.xtc -o gyrate.xvg -n index.ndx -tu ns

Extraction of aptamer in a stable state for ligand docking:

gmx trjconv -s md_0_10.tpr -f md_0_10_noPBC.xtc -o [aptamer name_extracted].pdb -dump 50000 -n index.ndx

Ligand complex MD commands

Generating topology:

Aptamer topology:

gmx pdb2gmx -f [aptamer name].pdb -o [aptamer name]_processed.gro -ter

Ligand topology:

acpype -i [aptamer ligand].mol2

(extract ligand_GMX.top | ligand_GMX.itp | ligand_GMX.gro | posre_ligand.itp from the folder made)

Building the complex:

1. Create a duplicate of [aptamer name]_processed.gro and rename it to complex.gro

2. Copy the coordinate section of ligand_GMX.gro and paste it into complex.gro

3. Increment the second line of complex.gro by adding the number of atoms in the ligand to the original number of the aptamer.

Building the topology:

1. Copy "Include ligand topology" and "Include ligand parameters" from the ligand_GMX.top and paste it into the topol.top file generated earlier

2. Add "UNL" with "1 mols" at the bottom of the topol.top file

Defining Unit Cell & Adding Solvent:

gmx editconf -f complex.gro -o newbox.gro -c -d 1.2 -bt cubic

gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

Adding ions:

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr

gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

Energy minimization: Assemble the binary input using gmx grompp

gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr

gmx mdrun -v -deffnm em

Thermostats:

gmx make_ndx -f em.gro -o index.ndx

(select: RNA | UNL )

Equilibration - nvt : Assemble the binary input using gmx grompp

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

gmx mdrun -deffnm nvt

Equilibration - npt : Assemble the binary input using gmx grompp

gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr

gmx mdrun -deffnm npt

MD simulation:

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_10.tpr

gmx mdrun -deffnm md_0_10 -v

Analysis:

Correcting for Periodicity Effects:

gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact -n index.ndx

RMSD graph:

gmx rms -s md_0_10.tpr -f md_0_10_center.xtc -o rmsd.xvg -tu ns -n index.ndx

Radius of gyration:

gmx gyrate -s md_0_10.tpr -f md_0_10_center.xtc -o gyrate.xvg -n index.ndx -tu ns

Distance graph:

Extracting frames from the trajectory:

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o conf.gro -sep

Measure the distance of each of these frames between the aptamer and complex

Bash get_distances.sh ; written by Justin A. Lemkul, Ph.D.

Plot a graph for it via Xmgrace

Xmgrace summary_distances.dat