
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:
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.
The following workflow was implemented to validate the RNA aptamer's function and binding mechanism, moving from structure prediction to dynamic simulation.
The initial phase converts the aptamer's linear sequence into a functional three-dimensional model, which is essential for predicting binding behavior.
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.
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.
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.
This phase predicts the initial binding event and provides an empirical score for aptamer selection.
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.
This is the final, most rigorous validation step, which models the dynamic behavior of the entire system over time.
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.
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.
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.
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.
![]() |
![]() |
| 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.
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.
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:
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:
We performed MD runs on the aptamer only and obtained RMSD graphs to analyze the stability of the aptamers.
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.
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.
![]() |
![]() |
| Figure 2a. RMSD graph of Glyphosate Y aptamer | Figure 2b. RMSD graph of Glyphosate X aptamer |
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.
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.
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.
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.
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.
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.
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.
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.
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.
After we completed testing the glyphosate aptamers (Glyphosate X and Y), we proceeded to test the other aptamers in the study.
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.
![]() |
![]() |
| 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 |
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.
![]() |
![]() |
| 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 |
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.
However, the docking results indicate that these aptamers exhibit a weaker binding ability compared to the glyphosate aptamers we have tested.
Lastly, we performed the MD runs on the Aptamer-Ligand Complexes, and the results are as follows:
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
The model can be expressed as:
RGB = R + G + B
Where:
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.
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).
[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.
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.
| 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 |
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.
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:
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.
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.
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.
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.
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.
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.
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.
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).
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.
[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.
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
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