In molecular model, we focused on the structural prediction and functional analysis of key proteins involved in the gene circuits of engineered bacteria. The first task involved predicting the 3D structure of aconitate isomerase (TbrA) using AlphaFold 2.0, followed by molecular dynamics simulations at various temperatures to assess substrate affinity. Molecular docking was performed using Rosetta and MOE to evaluate interactions between cis-aconitate and the enzyme.

Additionally, we predicted the HrpRS-Phrpl complex structure using AlphaFold 3 and molecular docking strategies, followed by binding free energy calculations to estimate the dissociation constant (K_d). This work aimed to understand the interaction between the HrpRS protein complex and the Phrpl promoter.

Lastly, to improve the sensitivity of the NahR protein to salicylic acid (SA), we performed saturation mutagenesis around the binding site, generating 5000 mutant sequences. After several rounds of screening based on energy calculations, we identified 22 mutant sequences with enhanced binding affinity to SA.

Why we built the model?

Since the crystal structure of the chosen TbrA has not been resolved, we aim to predict its structure using bioinformatics methods to guide subsequent modeling efforts.

How will this model be implemented?

First, we obtained the tbrA gene sequence encoding TbrA from Bacillus thuringiensis YBT-1520 and translated its coding sequence(CDS) into its corresponding amino acid sequence. We initially attempted to reconstruct the 3D structure of TbrA through homology modeling based on this amino acid sequence. However, after performing a BLAST search in the PDB database, the highest sequence similarity we found was only 41.8%. Therefore, we concluded that there is a lack of suitable templates among the resolved protein structures. Ultimately, we chose to predict the 3D structure of TbrA using AlphaFold 2.0 [1] and downloaded the best-predicted conformation.

Conclusion

From the Ramachandran plot, it can be observed that the majority of amino acid residues in the TbrA structure predicted by AlphaFold are within the allowed regions. Only four residues fall outside of these regions. However, since molecular dynamics simulations will be conducted later, we still consider this to be a reliable prediction.

Why we built the model?

In natural environments, the temperature of shallow soils exhibits significant variation due to diurnal cycles and seasonal fluctuations. These temperature differences may influence enzyme activity in practical applications, thereby affecting the overall effectiveness of their use.

What question we want to answer?

How will this model be implemented?

Ⅰ.Molecular dynamics simulations

First we used AmberTools[2] to conduct molecular dynamics simulations at 273K, 288K, and 301K (RT) to obtain the structures closest to their natural states. Besides, we employed MOE for molecular docking and analyzed the affinity between the substrate and enzyme.

Before running the molecular dynamics simulations, we completed the following steps:

  1. Selected the AMBER14SB force field.
  2. Chose the TIP3P model.
  3. Added an 8Å water box for protein solvation.
  4. Minimized the energy of the dynamical system.

We then performed molecular dynamics simulations at temperatures(1500ps) of 273K, 288K, and 301K (RT).

The MD trajectories were run on the bioinformatics server provided by Professor Wang Xiaocong, and we generated root mean square deviation (RMSD) graphs for the trajectories at different temperatures.

Ⅱ.Molecular docking

As the mechanism of aconitase isomerase catalyzing the isomerization of CAA to trans-Aconitic acid(TAA) has not been resolved, we proposed two docking strategies:

  1.Blind Docking: Aiming for optimal binding sites by increasing the number of structures screened.
  2.Site Prediction: Using the Site Finder module in MOE[3] to predict potential binding sites on TbrA, selecting the site with the highest PLB score as the docking site.

Based on the above two strategies, we performed molecular docking using MOE, downloading the ionized form of CAA from PubChem for the docking process.

Table 1. Docking score of MOE

Conclusion

During the molecular docking process, we observed that whether through blind docking or docking after predicting the binding site, the ligand consistently bound to the same cavity. Since the overall performance of docking after predicting the binding site was better, we chose the results from this approach.

After docking, we evaluated the top five structures based on the docking results and found that all ligands were generally located within the same cavity. We believe this cavity to be a potential active site of the isomerase.

Through structural alignment, we discovered that temperature had little impact on the structure around the binding site. However, the structural changes still affected the specific binding sites of CAA and TbrA. After analyzing the interactions between the ligand and receptor in the docking results, we found that at the temperatures we selected, TbrA could form relatively stable interactions with CAA. However, we were unable to identify specific key residues involved in the process of CAA isomerization to TAA.

Why we built the model?

In our circuit, signals from engineered bacteria responding to SA upregulation in plants are transmitted through a cascade involving the two-component system HrpR, HrpS, and the P_hrpl. We performed detailed modeling analyses of this process in the bacteria lysis model. However, despite an extensive literature review, we found no relevant data on the binding interaction between the HrpRS protein complex and the P_hrpl.

What question we want to answer?

How will this model be implemented?

Ⅰ.Preparation of HrpRS-P_hrpl complex structures

We used the online server of AlphaFol 3[4] to predict the structure of the HrpR, HrpS, and P_hrpl complex.

Ⅱ.Calculation of binding free energy

We use the Interface Analyzer module in Rosetta[5] to calculate several metrics, including the binding free energy (ΔG), for the HrpRS complex with P_hrpl based on the structure predicted by AlphaFold 3.

For the reaction:

$$\text R + L \stackrel{k_{\text {on}}}{\underset{k_{\text {off}}}{\rightleftharpoons}} RL $$

when the reaction reaches a steady state, we can derive the fyzaerlowing kinetic equation:

$$k_{on}[R][L]=k_{off}[RL]$$

The binding constant(K_d) is given by:

$$K_{d}=\frac{k_{off}}{k_{on}}$$

There is a quantitative relationship between the dissociation constant and the molar Gibbs free energy (ΔG, binding free energy):

$$ΔG=RTlnK_{d}$$

Conclusion

Through calculations using the Interface Analyzer, we obtained the binding free energy (ΔG) between HrpRS and P_hrpl, which is -36.602 kJ/mol. Based on this value, the dissociation constant (K_d) was calculated to be 3.78 × 10^-7 M , providing a reasonable parameter for use in the bacteria lysis model.(Click here learn more about this model)

Why we built the model?

In expert interviews conducted by the HP group, Professor Shunping Yan from the College of Life Sciences and Technology at Huazhong Agricultural University mentioned that "the concentration of SA varies across different plant tissues". Most available information on SA concentration is related to leaves, as the synthesis pathway of SA in plants is closely linked to chloroplasts.

Therefore, we hypothesize that the concentration of SA in plant roots is relatively low. To enhance NahR's sensitivity to SA in order to address potential variations in SA concentration, we performed saturation mutagenesis on residues near the docking site based on NahR's affinity for SA.

What question we want to answer?

We wanted to screen out mutants with higher affinity for SA through virtual mutation, hoping to improve the sensitivity of NahR to SA through this method.

How will this model be implemented?

Ⅰ.Generation of NahR-SA Complex

Before performing saturation mutagenesis, we first need to identify the ligand binding site and prepare a conformational library of SA isomers,the steps are as follows.

1.We initially used MOE to conduct preliminary docking of SA to identify its binding site with NahR.

2.We utilized the BCL::Conf web server[6] to generate a conformational library of SA isomers, obtaining all possible structures of SA for molecular docking.

3.The ligands were placed into the binding pocket of the protein.

Ⅱ.Rosetta Design

1.Generation of mutant structures

After obtaining the complex structure, we selected residues near the docking site as mutation sites for saturation mutagenesis. The mutant structures were generated using RosettaLigand[7][8] based on RosettaScript [9], in order to obtain 5000 mutants along with their energy and other parameter information.

2.Preliminary screening of mutants

Using the total scores from the energy information generated by RosettaLigand, we performed a preliminary screening based on the van der Waals interactions between the protein and the ligand, as well as among amino acid side chains. This process yielded 3681 sequences.

3.Refinement of mutant energy calculations and selection

In this step, we employed the Interface Analyzer to conduct a more detailed calculation of interactions between the mutants and SA from the preliminary screening results. Based on the calculated results, we filtered the mutants using metrics such as packstat and delta unsatHbonds, progressively identifying those with improved binding affinity for SA. After nine rounds of selection, we ultimately obtained 22 mutant sequences with enhanced binding capabilities for SA.

Conclusion

Based on the interactions between the ligand and the protein, we screened 22 NahR mutants with enhanced binding affinity to SA from the 5000 mutants generated by RosettaLigand.

Future work & cooperation

Initially, we planned to measure the binding ability of NahR mutants to SA through Surface Plasmon Resonance (SPR) and verify its function as a biosensor. However, unfortunately, due to some problems in the functional characterization of wild-type SA biosensors, we did not have enough time to verify it through experiments.

Click here to download the mutant sequences.