Molecular dynamics simulations reveal the selectivity mechanism of structurally similar agonists to TLR7 and TLR8.

TLR7 and TLR8 are key members of the Toll-like receptor family, playing crucial roles in the signaling pathways of innate immunity, and thus become attractive therapeutic targets of many diseases including infections and cancer. Although TLR7 and TLR8 show a high degree of sequence homology, their biological response to small molecule binding is very different. Aiming to understand the mechanism of selective profiles of small molecule modulators against TLR7 and TLR8, we carried out molecular dynamic simulations on three imidazoquinoline derivatives bound to the receptors separately. They are Resiquimod (R), Hybrid-2 (H), and Gardiquimod (G), selective agonists of TLR7 and TLR8. Our MD trajectories indicated that in the complex of TLR7-R and TLR7-G, the two chains forming the TLR7 dimer tended to remain "open" conformation, while the rest systems maintained in the closed format. The agonists R, H, and G developed conformational deviation mainly on the aliphatic tail. Furthermore, we attempted to quantify the selectivity between TLR7 and TLR8 by binding free energies via MM-GBSA method. It showed that the three selected modulators were more favorable for TLR7 than TLR8, and the ranking from the strongest to the weakest was H, R and G, aligning well with experimental data. In the TLR7, the flexible and hydrophobic aliphatic side chain of H has stronger van der Waals interactions with V381 and F351 but only pick up interaction with one amino acid residue i.e. Y353 of TLR8. Unsurprisingly, the positively charged side chain of G has less favorable interaction with I585 of TLR7 and V573 of TLR8 explaining G is weak agonist of both TLR7 and TLR8. All three imidazoquinoline derivatives can form stable hydrogen bonds with D555 of TLR7 and the corresponding D543 of TLR8. In brief, the set of total 400ns MD studies sheds light on the potential selectivity mechanisms of agonists towards TLR7 and TLR8, indicating the van der Waals interaction as the driving force for the agonists binding, thus provides us insights for designing more potent and selective modulators to cooperate with the hydrophobic nature of the binding pocket.

1. Introduction

Toll-like receptors (TLRs) are a large family of proteins, playing an important part in innate immune system, that recognize structurally conserved molecules, such as single-stranded (ss) or double-stranded (ds) RNAs or DNAs, lipoproteins and lipopolysaccharides derived from microbes and then activate immune cell responses [1]. A typical TLR is a single-spanning receptor consisting of three domains: an extracellular domain (ECD) with variable number of leucine-rich repeat sequences (LRRs) for the recognition of pathogen-associated molecular patterns (PAMPs), a transmembrane domain (TMD), and an intracellular Toll-interleukin 1 receptor (TIR) domain initiating downstream signaling [2]. Until now, thirteen TLRs have been identified, and among which, TLR3, TLR7, TLR8, and TLR9 are located in intracellular membranes because they are sensors of nucleic acids. Specifically, TLR3 recognizes viral dsRNA, and TLR9 senses unmethylated cytosine phosphate guanosine (CpG) containing DNA, whereas TLR7 and TLR8 both are located in endosomal membrane and function as viral ssRNA sensors [3].

Many studies have revealed that the expression levels of TLR7 and TLR8 are altered in some autoimmune diseases, such as arthritis, cancers [4-8], or in antiviral regimes, including corona virus prevention and HIV [9]. Thus, novel drug design and development against TLR7 or TLR8 became very attractive.

In the last few years, many TLR7 or TLR8 agonists with different scaffolds have been developed. These agonists leading to the induction of certain IFNs, cytokines and chemokines can be applied to the treatment of some diseases and can be used as good adjutants of vaccines [10-12].

The drug development of TLR modulators requires a solid understanding of TLR7 and TLR8 activity regulation. TLR7 and TLR8, sharing high degree of sequence homology and three dimensional structure similarity, are both known to serve as endosomal pattern recognition receptors (PRRs) for a number of RNA viruses, such as HIV, coronaviruses, influenza etc [9, 13–16]. However, there still exist many distributional and functional differences between these two closely related proteins. TLR7 is mainly expressed in plasmacytoid dendritic cells (pDC) and B cells [17-19]. But TLR8 is mainly expressed in myeloid dendritic cells (mDC), monocytes, macrophages and neutrophils [17, 20, 21]. TLR7 recognizes guanosine and its derivatives [22]while TLR8 serves as a uridine receptor [23]. Moreover, interferons induced by pDC are the major production of TLR7 signaling [24], whereas TLR8 signaling mainly results in NF-kB pathway activation and subsequent proinflammatory cytokines and chemokines expression [25]. To understand how TLR7 and TLR8 can recognize different ligands, and as consequence to activate different signaling pathways, we design computational simulations to investigate the selectivity mechanisms of small molecule agonists of TLR7 and TLR8.

Molecular dynamics (MD) simulation is a very established computational technique to understand the protein structure-function relationship and guide the drug design. MD simulation also has played an important role in characterizing receptor-ligand interaction [26-28], providing the guidance of structure-based drug design [29-31] and typical conformations for virtual screening [32-37]. Besides, it helps to reveal the novel binding sites which have not been captured by NMR and X-ray crystallographic analysis, for example, cryptic binding sites in HIV-1 integrase [38-41]. Up to date, MD simulations have been successfully applied to many large systems, such as the complete HIV1 capsid with 64 million atoms up to 100ns [42-45]. MD simulation has been also instrumental on understanding the protein folding and function regulation with a simulation time of 10-100us [46-48].

Because of the above mentioned success, not surprisingly, some MD simulations were carried out on TLRs, such as the stability of vaccine and TLRs [49, 50], TLRs model [51, 52] and effect of mutations on TLRs dimer [53]. MD simulation was applied to equilibrate homology model of TLR7, proposing the appropriate TLR7 dimer structure and studying the binding site and residues significant for dimerization [54]. And it was used to explain the difference of interaction mode between agonists and antagonists (including imidazoquinoline and adenine derivatives) for TLR7 [55]. Targeted molecular dynamics (TMD) was employed to study the conformational transition of TLR8 dimerization, illuminating the internal mechanism of relatively aggregate movement of two TLR8 chains [56]. These research results have inspired our research ideas.

Some imidazoquinoline derivatives (Fig 1), i.e. Resiquimod (R) [57, 58], Hybrid-2 (H) [59] and Gardiquimod (G) [60] are agonists of TLR7 and TLR8. They carry the common core of the same parent nucleus, whereas the side chain is oxygen, carbon, nitrogen atom, respectively. However, this one-atom difference in the chemical structure results in magnitude difference in their biological activity on TLR7 and TLR8. In general, they are more potent on TLR7 than TLR8 (Table 1). These two facts aroused our interest in investigating the agonist selectivity for TLR7 and TLR8 by MD simulations.

Fig 1

The chemical structures of the three agonists: Resiquimod (R), Hybrid-2 (H) and Gardiquimod (G).

Table 1

Comparison of predicted binding free energies and experimental EC50 values of R, H and G with TLR7 and TLR8.

Compound	TLR7		TLR8
	Predicted binding free energy (kcal/mol) a	EC50 (nM)	Predicted binding free energy (kcal/mol) a	EC50 (nM)
R	-44.21 ± 0.24	1400 [78]	-39.55 ± 0.20	6400 [78]
H	-50.50 ± 0.25	2.5 [59]	-42.09 ± 0.28	19 [59]
G	-33.65 ± 0.31	2000 [79]	-20.99 ± 0.19	No activation of NF-κB [80]

a The predicted binding free energies were obtained based on 30–50 ns MD simulation trajectory. A total of 201 snapshots evenly extracted from the 30–50 ns MD trajectory of each complex system were used for MM-GBSA calculations.

In this work, initially, eight systems were built, and they are TLR7 (apo), TLR7-R, TLR7-H, TLR7-G, TLR8 (apo), TLR8-R, TLR8-H, and TLR8-G. Thereafter, 50ns MD simulations of each of the eight systems were performed. Subsequently, the binding free energy of each ligand with TLR7, and TLR8 was calculated respectively using the MM-GBSA method. The interactions between the three agonists with TLR7 and TLR8 were analyzed to explain the intrinsic mechanism of the selectivity of the three agonists at the atomic level.

2. Materials and methods

2.1 Research systems

To investigate the selectivity of R, H and G with TLR7 and TLR8, respectively, the eight systems for the study are TLR7 (apo), TLR7-R, TLR7-H, TLR7-G, TLR8 (apo), TLR8-R, TLR8-H and TLR8-G. The crystal structures obtained from the Protein Data Bank (PDB) were used as the templates, and the codes are 5GMH (TLR7-R) [22], 5ZSG (TLR7-G), 3W3N (TLR8-R) [61], 4R6A (TLR8-H). The apo TLR7 was modeled from the X-ray complex of TLR7-R (5GMH) by removing the ligand and same to the apo TLR8 from TLR8-R (3W3N). The TLR7-H system was modeled based on coordinates of TLR7-R by replacing the oxygen on aliphatic chain of R to carbon atom. Similarly, the TLR8-G system was modeled by replacing the oxygen atom on aliphatic chain of R in TLR8-R with a nitrogen atom. The crystal structures were connected for the missing residues via SWISS-MODEL server [62]. The modeled structures were aligned to their crystal structures and were retained with corresponding ligand. The sequence fragment between amino acid residue 434–458 was removed to align with the biological understanding that those motifs in both TLR7 and TLR8 are cleaved prior to the activation process [61, 63]. In fact, these residues are also missing in 3W3N and 4R6A crystal structures. With reference on the recent research that human TLR7 and TLR8 proteins are in acidic endolysosom (pH is around 5) [64, 65], the protonation states of agonists were predicted by Maestro (from Schrödinger V.2019-4) [66] at the pH value of 5 (±0.5). The nitrogen atoms on quinoline ring in all complex systems and nitrogen atoms on aliphatic chain in TLR7-G system and TLR8-G system are protonated.

2.2 Molecular dynamics simulations

All MD simulations were carried out in the isothermal isobaric (NPT) ensemble with periodic boundary condition by the program GROMACS (version 2020.4-GPU package) [67]. The AMBER ff99SB force filed [68] was applied to model proteins and the general amber force filed (GAFF) [69] was assigned to model the small molecule agonists. Each system was placed in a rectangular box of TIP3P explicit water molecules with a minimum distance to the water box wall of 9 Å [70], and counterions (Cl-) were added to neutralize the system. Each simulation system was first subjected to energy minimization using the steepest descent algorithm. Then, a simulation was carried out to heat the system to 300K with the protein fixed using a harmonic restraint. The temperature was kept close to 300K by V-rescale thermostat [71] and the pressure was kept at 1 bar using the Parrinello-Rahman pressure coupling scheme [72]. The LINCS method [73] was used to restrain bond lengths that including hydrogen atoms, allowing an integration step of 2 fs. Long-range electrostatic interaction was calculated using PME with the cutoff 9Å. Finally, based on the relaxed system, the simulation so called the production phase was performed without any constraints for 50ns.

2.3 Binding free energy calculation

The MM-GBSA binding free energy calculations between the protein and agonist were using MMPBSA.py module in AmberTools20 package. The binding free energy (ΔGbinding) between a receptor and a ligand can be estimated using Eqs 1 and 2: ΔGbinding: the relative binding free energy. ΔEMM: the gas phase energy consisting of electrostatic (ΔEelec) and van der Waals (ΔEvdW) terms. ΔEsolv: solvation energy including both the polar solvation energy, ΔEpolar and the nonpolar solvation component, ΔEsurf. The ΔEpolar in above equation is calculated by the GB model [74] and ΔEsurf is estimated by the solvent accessible surface area (SASA). TΔS: the entropy term. It is usually suitable for systems with large conformational changes, thus was ignored in our simulations, aligning with other previous computational studies [75, 76].

3. Result and discussion

3.1 Characterization of conformational changes in the systems

The conformational changes of the eight simulated systems, ie TLR7 (apo), TLR7-R, TLR7-H, TLR7-G, TLR8 (apo), TLR8-R, TLR8-H, TLR8-G systems were first analyzed in terms of the root mean square deviation (RMSD) of protein backbone in each system (Fig 2). The RMSD values of TLR7 (Fig 2A) in the three complex systems change from the initial value (0 ns) to a scope of 2.5–4.0 Å, maintaining in this scope after 30 ns. Among these three, TLR7-G system, represented by the green line has the largest RMSD value close to 4 Å, and therefore the largest deviation comparing with the rest agonists. Throughout most of the simulation, TLR7-R and TLR7-H systems have similar RMSD plots to one another.

Fig 2

The overall RMSD values of TLR7 and TLR8 with and without the agonists.

(A) TLR7 (apo) system, TLR7-R system, TLR7-H system and TLR7-G system. (B) TLR8 (apo) system, TLR8-R system, TLR8-H system and TLR8-G system.

The RMSD values of TLR8 (Fig 2B) in the three complex systems change from the initial value (0 ns) to a scope of 2.0–3.0 Å and then also maintain in this scope after 30 ns. Among these three, TLR8-R, TLR8-H and TLR8-G systems have the similar RMSD plot during most of the simulation time. These reveal the larger conformational change of complex systems presented on TLR7 compared to that of TLR8. It can be observed that the RMSD value in the TLR8 (apo) system is slightly higher than those of the TLR8 complex systems after 30 ns, which indicates that the TLR8 is more flexible in the absence of agonist. The RMSD value in the TLR7 (apo) system is slightly lower than those of the TLR7 complex system after 30 ns, which indicates that the TLR7 is less flexible in the absence of agonist.

In addition, the trajectories displayed in Fig 3 using VMD software [77] to observe conformational changes show that the two chains (i.e. chain A and chain B) in TLR7 of the TLR7-R and TLR7-G systems were gradually open measured by the distances between the anchor points in the first 30 ns, but the TLR7-H system maintained the closed conformation state. In order to describe the closed and open conformational changes, the centroid of R784 backbone atoms in chain A and chain B of TLR7 were defined as anchor points to calculate the time-dependent distance between two chains (Fig 3A). In the same way, the centroid of P773, the corresponding counterpart in TLR8 of R784 in TLR7, backbone atoms in chain A and chain B of TLR8 were chosen as anchor points (Fig 3B). These two pairs of amino acid residues were chosen as anchor points because they are located on the central axis and closest to the bottom of the two chains. In complex systems, the distance between the two R784 gradually increased in TLR7-R and TLR7-G systems, though the distance in TLR7-G is longer than TLR7-R (Fig 3C). The distance differences between the structures during the last 20 ns MD simulations and initial structure are 7.4 ± 0.15 Å and 7.1 ± 0.13 Å in TLR7-R and TLR7-G systems, respectively. However, the distance changes in TLR7-H, TLR8-R, TLR8-H and TLR8-G systems are not evident (Fig 3C and 3D), which are 2.8 ± 0.14 Å, 1.8 ± 0.08 Å, 3.0 ± 0.12 Å and 2.3 ± 0.05 Å, respectively. In apo systems, the distance changes in TLR7 and TLR8 systems are also not evident (Fig 3C and 3D), which are 3.6 ± 0.15 Å and 0.6 ± 0.08 Å, respectively.

Fig 3

The distance changes of chain A and chain B for TLR7 and TLR8.

(A) The R784 of each chain is defined as anchor point on TLR7. (B) The P773, counterpart of R784 of TLR7, of each chain is defined as anchor point on TLR8. (C) The time dependent distance between chainA-R784 and chainB-R784. (D) The time dependent distance between chainA-P773 and chainB-P773. The Δ_distance refers to the difference of the chain A and chain B distance between the structures during the last 20 ns MD simulations and initial structure. The initial distances (0 ns) are 15.1 Å, 12.4 Å, 16.5 Å and 15.3 Å in TLR7 (apo), TLR7-R, TLR7-H, TLR7-G systems, respectively. The initial distances (0ns) are 11.9 Å, 8.3 Å, 11.9 Å and 14.1 Å in TLR8 (apo), TLR8-R, TLR8-H, TLR8-G systems, respectively. The average distances during the last 20 ns are 18.7 ± 0.15 Å, 19.8 ± 0.15 Å, 19.3 ± 0.14 Å and 22.4 ± 0.13 Å in TLR7 (apo), TLR7-R, TLR7-H, TLR7-G systems, respectively. The average distances during the last 20 ns are 11.3 ± 0.08 Å, 10.1 ± 0.08 Å, 14.9 ± 0.12 Å and 11.8 ± 0.05 Å in TLR8 (apo), TLR8-R, TLR8-H, TLR8-G systems, respectively.

3.2 Conformational changes of pocket residues and agonists

The amino acid residues within 6 Å of the agonists were defined as the pocket residues in this work. The RMSD values of backbone atoms (Fig 4A and 4B) and heavy atoms (Fig 4C and 4D) of pocket residues were calculated respectively. The RMSD values of backbone atoms in the four TLR7 systems change from the initial value (0 ns) to a scope of 1.0–2.0 Å, and 1.5–2.5 Å for heavy atoms, maintaining in this range after 30 ns. In these four systems, the RMSD values of pockets in the TLR7 (apo) and TLR7-H systems are the lowest and the RMSD values of residues in TLR7-R and TLR7-G systems are nearly the same. It reveals that the pocket of TLR7 is less flexible in the TLR7 (apo) and TLR7-H systems. The RMSD values of backbone atoms and heavy atoms in the four TLR8 systems change from the initial value (0 ns) to a scope of 1.2–2.5 Å and 1.5–3.0 Å, maintaining in this scope after 30 ns. In these four systems, the RMSD value of pocket residues in the TLR8-H system is the lowest of the four, while the RMSD values of pocket residues in the TLR8 (apo), TLR8-R and TLR8-G systems are nearly the same. It reveals that the pocket of TLR8 is less flexible in the TLR8-H system compared to TLR8 (apo), TLR8-R and TLR8-G systems. The above analysis suggests that the pockets residues of TLR8 were more flexible than that of TLR7, though the overall conformation change of TLR8 is smaller than that of TLR7 (Fig 3).

Fig 4

The RMSD values of the pocket residues in each system defined by 6 Å away from the agonists.

(A) RMSD of backbone atoms of TLR7 pocket. (B) RMSD of backbone atoms of TLR8 pocket. (C) RMSD of heavy atoms of TLR7 pocket. (D) RMSD of heavy atoms of TLR8 pocket.

We also examined the conformational changes of agonists. The heavy atoms’ RMSD of agonists was calculated in the six TLR-agonist systems (Fig 5B and 5E) with respect to the initial conformation. The RMSD of agonists in TLR7 and TLR8 complex systems increase from the initial value (0 ns) to a scope of 1.0–3.0 Å and maintain in this scope after 30 ns. However, three agonists show significant differences in fluctuation range. The fluctuations are largest in TLR7-G and TLR8-G systems, and the fluctuations are smallest in TLR7-R and TLR8-R. The conformations of the three agonists were superimposed on the imidazoquinoline ring, the initial frame (0 ns) and a frame near the late stage of the simulation (40 ns) were selected for comparison. In the initial frame, all atoms of the agonists of TLR7 (Fig 5A) and TLR8 (Fig 5D) complex systems are superimposed nicely. The alignment at t = 40 ns of snapshot was chosen to represent every conformation in the simulation 30–50 ns, the side chain orientations of the ligands in TLR7 (Fig 5C) and TLR8 (Fig 5F) are apparently different.

Fig 5

The conformational change of the agonists.

Schematic diagram of three TLR7 agonists superimposed on the imidazoquinoline ring at 0 ns (A) and 40 ns (C). Schematic diagram of three TLR8 agonists superimposed on the imidazoquinoline ring at 0 ns (D) and 40 ns (F). RMSD values of agonists in TLR7 complex systems (B) and TLR8 complex systems (E).

3.3 MM-GBSA binding free energy

Aiming to quantify the selectivity profile between the agonists and TLR7 and TLR8, the last 20 ns of simulation data was applied for MM-GBSA binding free energy analysis. The predicted binding free energies of R, H and G binding to TLR7 and TLR8 are summarized in Table 1. The binding free energies of R, H and G binding to TLR7 are -44.21 kcal/mol, -50.50 kcal/mol and -33.65 kcal/mol, respectively. The binding free energies of R, H and G binding to TLR8 are -39.55 kcal/mol, -42.09 kcal/mol and -20.99 kcal/mol. The results indicate that H shows the strongest binding affinity, while G shows the lowest. Additionally, the binding affinity of each of the agonists with TLR7 is stronger than that of TLR8, consistent with previously published experimental data of R, H and G [59, 78–80].

To better understand the binding profile, key binding contributors were also analyzed and summarized in Table 2 for TLR7, and in Table 3 for TLR8. Specifically, van der Waals interaction makes a major contribution to binding free energy. In TLR7-H and TLR7-R systems, there are stronger van der Waals interactions compared to the TLR7-G system. In TLR8-H and TLR8-R systems, there are also stronger van der Waals interactions compared to the TLR8-G system. The van der Waals interactions of TLR7 complex systems are overall stronger than that of TLR8 complex systems. The energy is unfavorable to binding is the total of electrostatic interaction and polar solvation energy. TLR7-H and TLR8-H systems have the strongest binding affinity as compared to that of corresponding complex systems, because of lower van der Waals energies and contribution from the electrostatic interaction and polar solvation energy. The above indicates that the pockets of TLR7 and TLR8 are hydrophobic pockets. Furthermore, the energy decomposition values of binding free energy were calculated to get insight into the binding mode of the agonists with TLR7 and TLR8. Items that contributed less than -1.0 kcal/mol to the binding free energy are listed in S1-S6 Tables in S1 File. In S1 and S2 Figs in S1 File, the values of these main residues are displayed in a bar graph.

Table 2

Results of MM-GBSA method of R, H and G to TLR7.

Energetic contributions	TLR7-H	TLR7-R	TLR7-G
ΔE vdW	-45.72 ± 0.20	-45.24 ± 0.23	-34.60 ± 0.23
ΔE elec	398.07 ± 0.84	393.98 ± 0.63	788.50 ± 1.45
ΔE polar	-396.95 ± 0.79	-387.33 ± 0.57	-783.04 ± 1.35
ΔE surf	-5.90 ± 0.01	-5.63 ± 0.01	-4.52 ± 0.02
ΔE pol,ele	1.12	6.65	5.46
ΔG binding	-50.50 ± 0.25	-44.21 ± 0.24	-33.66 ± 0.31

bΔ Epol,ele is ΔEpolar plus ΔEelec. All units are kcal/mol.

Table 3

Results of MM-GBSA method of R, H and G to TLR8.

Energetic contributions	TLR8-H	TLR8-R	TLR8-G
ΔE vdW	-41.41 ± 0.22	-38.32 ± 0.19	-27.90 ± 0.19
ΔE elec	74.98 ± 0.76	27.72 ± 0.68	157.83 ± 1.26
ΔE polar	-70.10 ± 0.76	-23.78 ± 0.63	-146.77 ± 1.21
ΔE surf	-5.56 ± 0.02	-5.16 ± 0.02	-4.14 ± 0.02
ΔE pol,ele	4.88	3.94	11.06
ΔG binding	-42.09 ± 0.28	-39.55 ± 0.20	-20.99 ± 0.19

cΔ Epol,ele is ΔEpolar plus ΔEelec. All units are kcal/mol.

Fig 6A–6C depicts the interaction plot of H, R and G with TLR7 pocket residues. There are electrostatic interactions, H-bond interactions, π-π interactions, and C-H-π interactions between D555, T586, F408 and L557 with agonists. F351, Y356, V381 and I585 form the hydrophobic interaction regions. Compared with H and R, G forms a stronger H-bond interaction with Thr586 of TLR7 through a nitrogen atom. However, the movement of the side chain weakens the van der Waals interactions between G and F351, V381 and I585. Compared with R, the carbon atom on the side chain of H is more flexible and hydrophobic than the oxygen atom on the side chain, leading H to be more adapted to pocket environment. H forms stronger binding free energies with D555, V381 and F351 (-6.768, -1.450 and -3.503 kcal/mol) than R (-4.542, -0.964 and -2.849 kcal/mol).

Fig 6

The interactions between the agonists R, H and G against TLR7 and TLR8.

TLR7-H system (A), TLR7-R system (B), TLR7-G system (C), TLR8-H system (D), TLR8-R system (E), TLR8-G system (F). The electrostatic interactions are shown in magentas dashes, the H-bond interactions are shown in yellow dashes, the π-π interactions are shown in green dashes, and the C-H-π interactions are shown in red dashes. TLR7 and TLR8 are shown in white cartoon. Agonists and representative residues are shown in cyan and orange stick.

Fig 6D–6F depicts the interaction plot of H, R and G with of TLR8 pocket residues. There are electrostatic interactions, H-bond interactions and π-π interactions between D543, T574 and F405 with H and R. Y348, Y353, V378 and V573 form the hydrophobic interaction regions. Compared with H and R, G forms two electrostatic interactions with D545 and D543. However, this interaction weakens the H-bond interaction with T574 and Van der Waals interaction with Y348 and V573. Zhu et al proposed that G is an agonist of TLR7 but not of human TLR8, which in accordance with our modeling. Using NF-κB reporter assay to measure the activation of human TLR7 and human TLR8, they found that G only activated human TLR7, but not human TLR8 in Cos-7 cells and 293T cells [80]. Compared with R, the carbon atom on the side chain of H is more flexible and hydrophobic than the oxygen atom on the side chain, leading H to be more adapted to pocket environment. It forms stronger electrostatic interaction with D543 and stronger van der Waals interactions with Y353.

Agonists R, H and G form C-H-π interaction with L557 of TLR7, which has been also reported in other study for R [55]. However, in TLR8 the corresponding residue of L557 is D545, which could not form strong interactions with agonists due to its shorter side chain.

3.4 Hydrogen bond interaction and van der Waals interaction between TLR7/8 and agonists

To further explore the interactions between three agonists and the receptors, the occupancy of hydrogen bonds with more than 10% occupancy between the agonists and residues atoms were analyzed to discard the extremely weak hydrogen bond interaction. The important atoms on the three agonists are shown in Fig 7C. The details of the hydrogen bond between agonists and TLR7 are shown in Fig 7A; S7-S9 Tables in S1 File. In TLR7-H, TLR7-R and TLR7-G systems, the stable hydrogen bonds formed between D555 and N and N1 atoms of the agonists were maintained during the last 20 ns of the conformations. In TLR7-H and TLR7-R systems, the occupancy of hydrogen bonds between T586 and N1 atoms of the agonists (TLR7-H: 34.88%, TLR7-R: 43.83%) are similar. Nevertheless, the occupancy of the hydrogen bond between N1 and T586 of TLR7-G was very low (TLR7-G: 12.99%). Since N4 on G was protonated, a hydrogen bond formed between T586 of TLR7-G and N4, resulting in a change in the orientation of the side chain and hence the loss of the hydrogen bond between T586 of TLR7-G and N2.

Fig 7

Analysis of hydrogen bond interactions between three agonists and TLR7 and TLR8. Occupancy of hydrogen bonds between agonists and TLR7.

(A) and TLR8 (B). (C) The position of the important atoms on the agonists H, R and G.

The details of the hydrogen bonds between three agonists and TLR8 are shown in Fig 7B; S10-S12 in Tables in S1 File. In TLR8-H, TLR8-R and TLR8-G systems, the stable hydrogen bonds formed between D543 and N and N1 atoms of the agonists were maintained during the last 20 ns of the conformations. In TLR8-H and TLR8-R systems, the hydrogen bonds form between T574 of TLR8 and N1 (TLR8-H: 30.68%, TLR8-R: 25.29%). Since N4 on G was protonated, hydrogen bond formed between G572 of TLR8-G and N4, which changed the orientation of the side chain.

The occupancy of residues less than 5 Å away from the agonists were analyzed in detail. To analyze the difference in data, all occupancy less than 20% or close to 100% in all six systems were ignored. The TLR7-R (Fig 8C) and TLR8-R (Fig 8D) systems are chosen to present the position relationship between residues and agonists. As shown in Fig 8A, the occupancy of residues N265, F349, E352, L353, G354, G379, Y380, T406, N407, F466, Y579 and H587 are highest in TLR7-H, lower in TLR7-R, and lowest in TLR7-G. Among them, residues N265, F349, E352, L353, Q379, Y380, T406, N407 and F466 are located around the aliphatic tail. As shown in Fig 8B, the occupancy of residues F261, N262, F346, G351, G376, Y377, I403 and T574 are highest in TLR8-H, lower in TLR8-R and lowest in TLR8-G. Among them, residues F261, N262, F346, G376, Y377 and I403 are located around the aliphatic tail. The above indicates that the occupancy of residues around three aliphatic tails are quite different, resulting in the van der Waals interaction between three agonists and receptors to be ranked H, R, then G, from highest to lowest. Comparing Fig 8A and 8B, the residues F466 and E583 in TLR7 complex systems were located in the 5 Å range around the agonists for a certain period of time, while the occupancy of corresponding residues P463 and A571 in TLR8 complex system is zero. Further, the occupancy of residues Y264, N265, L353 and H587 in three TLR7 complex systems are higher than residues F261, N262, L350, and H575 in three TLR8 complex systems, respectively. In addition, the occupancy of residues E352, L556 and T586 in two or one TLR7 complex systems is higher than I349, F544 and T574 in TLR8 complex systems. For example, the occupancy of residues E352 in TLR7-H complex system and L556 in TLR7-R and TLR7-G complex systems and T586 in TLR7-G system is higher than I349, F544 and T574 in TLR8 corresponding systems. This explains why the van der Waals interaction and binding affinity between agonists and TLR7 is stronger than TLR8. The corresponding counterpart of the residue of TLR7 in Fig 8A is the residue of TLR8 in Fig 8B.

Fig 8

Occupancy of residues less than 5 Å away from the three agonists.

TLR7 complex systems (A) and TLR8 complex systems (B). Red labels of horizontal coordinate represent pockets residues around side chain of agonists. The residues on x-axis of (A) are in alignment with that of (B). Position relationship between these residues and agonist in TLR7-R system (C) and TLR8-R system (D). The residues are shown in surface and stick. Red surface represents negative charge and blue surface represents positive charge. Agonists are shown in cyan stick.

Conclusion

The aim of this research was to explore the intrinsic mechanisms underlying the selectivity of R, H and G for TLR7 and TLR8 at atomic level. MD simulations and MM-GBSA method were used to model the overall conformational changes and calculate the binding free energies between three agonists and the receptors TLR7 and TLR8. Trajectory analysis showed that TLR7-R and TLR7-G systems formed open conformations during the simulation, however, other systems kept in closed conformations. The pocket residues in TLR7 are conformationally less flexible than those in TLR8, suggesting tight binding in TLR7. This is confirmed by the predicted binding free energies via MMM-GBSA method. Moreover, the calculated binding free energies indicated that the three agonists are more sensitive for TLR7 than TLR8, and the rank of the binding free energy values are in agreement with the experimental EC50 values in the cellular assay. In brief, in the last 20 ns of the complex systems, the flexible and hydrophobic aliphatic side chain of H forms van der Waals interactions with V381 and F351 of TLR7 and Y353 of TLR8. The side chain nitrogen of G is positively charged in an acidic environment, leading to less favorable interactions with I585 of TLR7 and V573 of TLR8. Stable hydrogen bonds were formed between agonists and D555 of TLR7 and D543 of TLR8. The occupancy of residues around less than 5 Å away from three agonists is quite different, which account for the deviation of van der Waals interaction between agonists and receptors. An atomic difference on the aliphatic tail of each agonist results in the occupancy of residues and the change of van der Waals interaction. Thus, MD simulations provide explanation of differences in interaction modes of three agonists binding with TLR7 and TLR8 at the atomic level, paving the way for further design of more effective TLR7 and TLR8 modulators.

Energy decomposition of binding free energy of agonists to TLR7 and TLR8.

Details of hydrogen bonds between agonists and TLR7 and TLR8. The binding free energies between each residue of TLR7 and TLR8 and agonists.

(DOCX)

Click here for additional data file.

16 Dec 2021

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Anand Gaurav

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf.

2. Thank you for stating the following in the Competing Interests section:

[The authors have declared that no competing interests exist.]

We note that one or more of the authors are employed by a commercial company: Shanghai ChemPartner Co., Ltd., 576 Libing Road, Shanghai.

a) Please provide an amended Funding Statement declaring this commercial affiliation, as well as a statement regarding the Role of Funders in your study. If the funding organization did not play a role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript and only provided financial support in the form of authors' salaries and/or research materials, please review your statements relating to the author contributions, and ensure you have specifically and accurately indicated the role(s) that these authors had in your study. You can update author roles in the Author Contributions section of the online submission form.

Please also include the following statement within your amended Funding Statement.

“The funder provided support in the form of salaries for authors [insert relevant initials], but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of these authors are articulated in the ‘author contributions’ section.”

If your commercial affiliation did play a role in your study, please state and explain this role within your updated Funding Statement.

b) Please also provide an updated Competing Interests Statement declaring this commercial affiliation along with any other relevant declarations relating to employment, consultancy, patents, products in development, or marketed products, etc.

Within your Competing Interests Statement, please confirm that this commercial affiliation does not alter your adherence to all PLOS ONE policies on sharing data and materials by including the following statement: "This does not alter our adherence to  PLOS ONE policies on sharing data and materials.” (as detailed online in our guide for authors http://journals.plos.org/plosone/s/competing-interests) . If this adherence statement is not accurate and  there are restrictions on sharing of data and/or materials, please state these. Please note that we cannot proceed with consideration of your article until this information has been declared.

Please include both an updated Funding Statement and Competing Interests Statement in your cover letter. We will change the online submission form on your behalf

3. In your Data Availability statement, you have not specified where the minimal data set underlying the results described in your manuscript can be found. PLOS defines a study's minimal data set as the underlying data used to reach the conclusions drawn in the manuscript and any additional data required to replicate the reported study findings in their entirety. All PLOS journals require that the minimal data set be made fully available. For more information about our data policy, please see http://journals.plos.org/plosone/s/data-availability.

Upon re-submitting your revised manuscript, please upload your study’s minimal underlying data set as either Supporting Information files or to a stable, public repository and include the relevant URLs, DOIs, or accession numbers within your revised cover letter. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories. Any potentially identifying patient information must be fully anonymized.

Important: If there are ethical or legal restrictions to sharing your data publicly, please explain these restrictions in detail. Please see our guidelines for more information on what we consider unacceptable restrictions to publicly sharing data: http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions. Note that it is not acceptable for the authors to be the sole named individuals responsible for ensuring data access.

We will update your Data Availability statement to reflect the information you provide in your cover letter.

4. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: http://journals.plos.org/plosone/s/supporting-information.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The article entitled "Molecular dynamics simulations reveal the selectivity mechanism of structurally similar agonists to TLR7 and TLR8". Through molecular dynamics simulations, the authors aim to understand the mechanism of selective profiles of small molecule modulators against TLR7 and TLR8. The article falls within the scope of PLOS ONE, and may be considered for publication subject to rectification of the following issues:

1. Introduction is lengthy, kindly reduce the MD simulations part in the introduction section.

2. Page 14, line 295: However, the movement of the side chain weakens the van der Waals interactions between G and Phe351, Val381 and Ile585.

line 298: It forms stronger electrostatic interactions with Asp555 and stronger van der Waals interaction with Val381 and Phe351.

These two statements are contradictory. In the first statement, you mentioned the movement of the side chain weakens the van der Waals interactions between G and Phe351 and Val381. This means van der Waals interactions are weak in this complex.

In the second statement, you mentioned molecule-H of G forms stronger van der Waals interaction with Val381 and Phe351. The van der Waals interactions in TLR7-G complex are weak or strong interactions? Kindly re-write this part.

3. In Fig6C, the residues Val381 and Ile585 do not appear in the figure, kindly include these two residues in the figure if they are important for the interactions as mentioned in page 14, line 296 and 299.

4. In Fig6F, the residues Thr574 and Val573 do not appear in the figure, kindly include these two residues in the figure as you mentioned there are H-bond interaction with Thr574 and Van der Waals interaction with Val573, page 14, line 304 and 305.

5. Page 14, line 302, you have mentioned interactions with Val378 of TLR8 but in Fig6D-E Val378 do not appear in any figure. Kindly include the residue in the figures.

6. There are many errors in English that should be improved. kindly use the same tense in your manuscript.

7. Errors need to be corrected as follow:

1. page 4, line 91: Up to now, MD simulations have been successfully applied to many large systems....million atoms.

Change to “Up to date, MD simulations have been successfully applied to many large systems....million atoms.”

2. Page 5, line 110: In general, they are more potent on TLR7 than those on TLR8 (Table 1).

Change to "In general, they are more potent on TLR7 than TLR8 (Table 1).

3. Page 7, line 160: The binding free energy (ΔGbinding) between a receptor and a ligand and can be estimated using Equation 1:

Change to "The binding free energy (ΔGbinding) between a receptor and a ligand can be estimated using Equation 1:"

4. Page 7, line 166: ΔGbinding: the relative binding free energy. Kindly move this sentence to the beginning of the paragraph (line 163) before ΔEMM : the gas phase...

5. Page 8, line 187: It can be observed that the RMSD value in the TLR8 (apo).....that the TLR8 is more flexible in the presence of agonist. Is it presence or absence, kindly check as higher RMSD is associated with protein flexibility.

6. Page 10, line 239: On the other hand of the trajectory analysis, we also examine the conformational change of agonists.

Change to "On the other hand of the trajectory analysis, we also examine the conformational changes of the agonists."

7. Page 11, line 241: The RMSD of agonists in TLR7 and TLR8 complex systems increase from the initial value (0 ns) to a scope of 1.0-3.0 Å and all reach maintain in this scope after 30 ns.

Change to "The RMSD of the agonists in TLR7 and TLR8 complex systems increase from the initial value (0 ns) to a scope of 1.0-3.0 Å and maintain in this scope after 30 ns.

8. Page 11, line 246: The conformations of the three agonists were superimposed on the imidazoquinoline ring, and the initial frame (0 ns) and a frame near the late stage of the simulation (40 ns) were selected for comparison.

Change to "The conformations of the three agonists were superimposed on the imidazoquinoline ring, the initial frame (0 ns) and a frame....comparison."

9. Page 11, line 248: The alignment at t= 40 ns of snapshot was chose to represent every conformation in the simulation 30-50 ns, the side chain orientations of the ligands in TLR7 (Fig 5C) and this in TLR8 (Fig 5F) are apparently different.

Change to "The alignment at t= 40 ns of snapshot was chosen to represent every conformation in the simulation 30-50 ns, the side chain orientations of the ligands in TLR7 (Fig 5C) and TLR8 (Fig 5F) are apparently different.

10. Page 15, line 324: Agonists and representative residues are shown in orange and cyan stick.

Change to "Agonists and representative residues are shown in cyan and orange stick."

11. Page 15, line 333: In TLR7-H, TLR7-R and TLR7-G systems, the stable hydrogen bonds formed in Asp555 and N and N1 atoms on agonists were maintained during the last 20 ns of the conformations.

Change to "In TLR7-H, TLR7-R and TLR7-G systems, the stable hydrogen bonds formed between Asp555 and N and N1 atoms of the agonists were maintained during the last 20 ns of the conformations."

12. Page 15, line 334: In TLR7-H and TLR7-R systems, the occupancy of hydrogen bonds between Thr586 and N1 atoms on agonists (TLR7-H: 34.88%, TLR7-R: 43.83%) are similar.

Change to "In TLR7-H and TLR7-R systems, the occupancy of hydrogen bonds between Thr586 and N1 atoms of the agonists (TLR7-H: 34.88%, TLR7-R: 43.83%) were similar."

13. Page 15, line 336: Nevertheless, the occupancy of the hydrogen bond between N1 and Thr586 of TLR7-G is very low (TLR7-G: 12.99%).

Change to "Nevertheless, the occupancy of the hydrogen bond between N1 and Thr586 of TLR7-G was very low (TLR7-G: 12.99%)."

14. Page 15, line 338: Since N4 on G was protonated, a hydrogen bond formed between Thr586 of TLR7-G and N4, resulting in a change the orientation of the side chain and the loss of the hydrogen bond between Thr586 of TLR7-G and N2.

Change to "Since N4 on G was protonated, a hydrogen bond formed between Thr586 of TLR7-G and N4, resulting in a change in the orientation of the side chain and loss of the hydrogen bond between Thr586 of TLR7-G and N2."

15. Page 15, line 342: The details of the hydrogen bonds... formed in Asp543 and N and N1 atoms on agonists were maintained during the last 20 ns of the conformations.

Change to "The details of the hydrogen bonds... formed between Asp543 and N and N1 atoms of the agonists were maintained during the last 20 ns of the conformations.

16. Page 16, line 344: The occupancy of hydrogen bond between N1 and Thr574 of TLR8-G is zero. Do you mean no hydrogen bond formed between Thr574 and N1 of G in TLR8-G complex?

17. Page 16, line 359: are located around the aliphatic tail. And as shown in Fig 8B, the....TLR8-G.

Change to "are located around the aliphatic tail. As shown in Fig 8B, the....TLR8-G."

18. Page 16, line 367: complex system is zero. And the occupancy...respectively.

Change to "complex system is zero. Further, the occupancy...respectively.

19. Page 17, line 388: MD simulations and MM-GBSA method were used to model the overall conformational changes and calculate the binding free energies between three agonists and the TLR7 and TLR8.

Change to "MD simulations and MM-GBSA method were used to model the overall conformational changes and calculate the binding free energies between three agonists and the receptors TLR7 and TLR8."

20. Page 17, line 389: Trajectory analysis showed that TLR7-R and TLR7-G systems form more “open” conformations during the simulation, however, other systems kept in closed conformations.

Change to "Trajectory analysis showed that TLR7-R and TLR7-G systems formed open conformations during the simulation, however, other systems kept in closed conformations.

21. Page 17, line 392-393: Plus, the calculated binding free energies indicated that three agonists are more sensitive for TLR7 than TLR8, and the rank of the binding free energy values are in agreement with the experimental EC50 values in the cellular assay.

Change to "Moreover, the calculated binding free energies indicated that the three agonists are more sensitive for TLR7 than TLR8, and the rank of the binding free energy values are in agreement with the experimental EC50 values in the cellular assay.

22. Page 18, line 397: The side chain nitrogen of G is positively charged in an acidic environment, leading to its much less favorite interaction with Ile585 of TLR7 and Val573 of TLR8.

Change to "The side chain nitrogen of G is positively charged in an acidic environment, leading to less favorite interaction with Ile585 of TLR7 and Val573 of TLR8.

Reviewer #2: The present work describes routine but useful simulations to understand the selectivity of agonists towards TLR7 and TLR8. The study was guided by binding free energy analysis, hydrogen bond analyses, decomposition of free energy. All the data have been presented either in main or supplementary text. The results obtained were almost aligned with the experimental results and helpful in determining the selectivity mechanisms. However the manuscript still require few explanations and language corrections.

Suggestions for clarification or modifications:

Line 13: replace “Highly” with “high”

Line 25: Replace “experiment” with “experimental”

Line 28: Replace “favor” with “favorable”

Line 29: Replace “agonist in” with “agonist of”

Line 29: Replace “imidazoquinolines” with “imidazoquinoline derivatives”

Line 31: Replace “selective” with “selectivity”

Line 32: Insights for “designing” more potent

Line 36: Molecular dynamics simulations

Line 50: both are located

Line 57: diseases and can be used as good adjutants of vaccines

Line 62: “many distributional and functional”

Line 71: “selectivity mechanism of small molecule agonists of TLR7 and TLR8d TLR8”

Line 80: pH

Line 89: replace “structural” to “structure”

Line 100: replace “against” with “for”

Line 111: replace “first” with initially

Line 112: remove “respectively”. Replace second with “thereafter”

Line 116: replace “selectivity mechanism” with “selectivity”

Line 128-29: The TLR7-H system was modeled based on coordinates of TLR7-R by replacing the oxygen to carbon atom. (Mention the aliphatic chain)

Section 2.2  change heading to “Molecular dynamics simulations”

Line 148: replace Cl-1 to Cl-

Line 149: “steepest descent algorithm”

Line 154: 50ns is not a long-time simulation. Just mention 50 ns simulation.

Section 2.3 equations must be included for “∆EMM : the gas phase energy consisting of electrostatic (∆Eelec) and van der Waals (∆EvdW) terms. ∆Gsolv : solvation energy including both the polar solvation energy, ∆Gpolar and the nonpolar solvation component, ∆Gsurf”

Line 192: “confirmational”

Line 197: capitalize “In”

Line 201:  remove “are both”.

Line 202: “The distance changes in TLR7-H, TLR8-R, TLR8-H and TLR8-G systems are not evident (Fig 3C and Fig 3D). In apo systems, the distance changes in TLR7 and TLR8 systems are also not evident (Fig 3C and Fig 3D).”

This point should be made clearer by showing the distance (angstrom) in Fig 3A and 3B. or by superimposing TLR7 and TLR8 with highlighted anchor residue distances.

Line 217: replace scope with “range”

Line 235: replace “On the other hand of the” by “In another”

& change to “changes”

Line 245: replace chose to “chosen”

Line 246: remove “this is”

Line 256: TLR7 and TLR8

Line 260: while G shows the lowest

Line 270: energy is unfavorable to binding include the values here.

Table 1. Comparison of predicted binding free energies and experimental EC50 values of R, H and G with TLR7 and TLR8

SEM (standard error of means) should be reported for all the BFE values e.g -44.21 ± ?? in all the tables.

Table 2. replace “Item” to “Energetic contributions”  All energetic contributions should have uniform representation e.g. ∆Evdw, ∆Eele……

All units are in kcal/mol.

Line 293: “It forms stronger electrostatic interactions with Asp555 and stronger van der Waals interaction with Val381 and Phe351” should be supported with energy values

Line 306-09: these are repeated lines. Remove.

Line 310: “The residue of TLR8 alignment to Leu557 of TLR7 is Asp545, which might result in a weaker interaction with agonists.” rewrite to make it more understandable.

Line 327: Similar codes (either 3 letter or single) for amino acids should be used in text and figures.

Line 332: side chain and hence the

Caption Fig 7.  “the important atoms”

Line 364: two or one TLR7 please specify

Line 382: remove “the” before TLR7 and TLR8.

Line 384: “pockets” to pocket

Line 392: less favorable interactions

Line 396: each agonist

Line 394: quiet to “quite”

Need to address these questions?

**********

How have you treated long range electrostatic?

Any specific reason for using only MM-GBSA (Not MM-PBSA) protocol for binding free energy calculations.

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

10 Feb 2022

Responses to reviewer#1

The article entitled "Molecular dynamics simulations reveal the selectivity mechanism of structurally similar agonists to TLR7 and TLR8". Through molecular dynamics simulations, the authors aim to understand the mechanism of selective profiles of small molecule modulators against TLR7 and TLR8. The article falls within the scope of PLOS ONE, and may be considered for publication subject to rectification of the following issues:

We greatly appreciate the reviewer’s valuable and constructive suggestions. We have incorporated all the suggestions of the reviewer that helped to improve this manuscript, and highlighted the revisions in red. We hope that you will find the revised manuscript suitable for publication in PLOS ONE.

1. Introduction is lengthy, kindly reduce the MD simulations part in the introduction section.

Response: Many thanks for the constructive suggestion. The MD simulations part has been reduced and revised paragraph in section 1 is following:

Molecular dynamics (MD) simulation is a very established computational technique to understand the protein structure-function relationship and guide the drug design. In the late 70s, the first case of MD is Bovine pancreatic trypsin inhibitor [26]. Especially in recent years, benefit from the development of modern hardware and force fields, MD simulation has made huge progress. Protein conformation dynamics that is hardly available by current experimental techniques can be observed using MD simulations [27-34]. MD simulation also has played an important role in characterizing receptor-ligand interaction [35-37], providing the guidance of structure-based drug design [38-40] and typical conformations for virtual screening [41-46]. Besides, it helps to reveal the novel binding sites which have not been captured by NMR and X-ray crystallographic analysis, for example, cryptic binding sites in HIV-1 integrase [47-50]. Up to date, MD simulations have been successfully applied to many large systems, such as the complete HIV1 capsid with 64 million atoms up to 100ns [51-54]. MD simulation has been also instrumental on understanding the protein folding and function regulation with a simulation time of 10-100us [55-57].

2. Page 14, line 295: However, the movement of the side chain weakens the van der Waals interactions between G and Phe351, Val381 and Ile585.

line 298: It forms stronger electrostatic interactions with Asp555 and stronger van der Waals interaction with Val381 and Phe351.

These two statements are contradictory. In the first statement, you mentioned the movement of the side chain weakens the van der Waals interactions between G and Phe351 and Val381. This means van der Waals interactions are weak in this complex.

In the second statement, you mentioned molecule-H of G forms stronger van der Waals interaction with Val381 and Phe351. The van der Waals interactions in TLR7-G complex are weak or strong interactions? Kindly re-write this part.

Response: We are thankful to the reviewer’s suggestions. The sentence “However, the movement of the side chain weakens the van der Waals interactions between G and Phe351, Val381 and Ile585.” describes the molecule G. However, the sentence “It forms stronger electrostatic interactions with Asp555 and stronger van der Waals interaction with Val381 and Phe351.” is following the sentence “The carbon atom on the side chain is more flexible and hydrophobic than the oxygen atom on the side chain, leading the molecule-H to be more adapted to pocket environment.”. Thus, the latter sentence describes the molecule H. To make this discussion clear, we revised these sentences as following:

However, the movement of the side chain weakens the van der Waals interactions between G and F351, V381 and I585. Compared with R, the carbon atom on the side chain of H is more flexible and hydrophobic than the oxygen atom on the side chain, leading H to be more adapted to pocket environment. H forms stronger electrostatic interactions with D555 and stronger van der Waals interaction with V381 and F351.

3. In Fig6C, the residues Val381 and Ile585 do not appear in the figure, kindly include these two residues in the figure if they are important for the interactions as mentioned in page 14, line 296 and 299.

Response: We are thankful to the reviewer’s suggestions. Figure 6 marks the important residues that have the binding energy of less than -1kcal/mol with the agonists. However, the movement of the side chain weakens the van der Waals interactions between G and Phe351, Val381 and Ile585. The binding energy between Val381 and Ile585 and G is more than -1kcal/mol (Val381: -0.190 kcal/mol, Ile585: -0.096 kcal/mol). Thus, residues Val381 and Ile585 are not shown in Figure 6C.

4. In Fig6F, the residues Thr574 and Val573 do not appear in the figure, kindly include these two residues in the figure as you mentioned there are H-bond interaction with Thr574 and Van der Waals interaction with Val573, page 14, line 304 and 305.

Response: Many thanks for the constructive suggestion. Figure 6 marks the important residues that have the binding energy of less than -1kcal/mol with the agonist. However, this interaction weakens the H-bond interaction with Thr574 and Van der Waals interaction with Tyr348 and Val573. The binding energy between Thr574 and Val573 and G is more than -1kcal/mol (Thr574: 0.089 kcal/mol, Val573: -0.211 kcal/mol). Thus, residues Thr574 and Val573 are not shown in Figure 6F.

5. Page 14, line 302, you have mentioned interactions with Val378 of TLR8 but in Fig6D-E Val378 do not appear in any figure. Kindly include the residue in the figures.

Response: Many thanks for the constructive suggestion. Val378 of TLR8 forms good van der Waals interaction with H (-1.372 kcal/mol), R (-1.204 kcal/mol) and G (-1.037kcal/mol). However, the binding free energies between Val378 and H, R and G are -0.790 kcal/mol, -0.911 kcal/mol and -0.622kcal/mol. These values are more than -1 kcal/mol, thus they are not marked in the Figure 6.

6. There are many errors in English that should be improved. Kindly use the same tense in your manuscript. Errors need to be corrected as follow:

Response: We greatly appreciate the reviewer’s valuable suggestions. We have double-checked our manuscript and revised it.

(1) Page 4, line 91: Up to now, MD simulations have been successfully applied to many large systems....million atoms.

Change to “Up to date, MD simulations have been successfully applied to many large systems....million atoms.”

(2) Page 5, line 110: In general, they are more potent on TLR7 than those on TLR8 (Table 1).

Change to "In general, they are more potent on TLR7 than TLR8 (Table 1).

(3) Page 7, line 160: The binding free energy (ΔGbinding) between a receptor and a ligand and can be estimated using Equation 1:

Change to "The binding free energy (ΔGbinding) between a receptor and a ligand can be estimated using Equation 1:"

(4) Page 7, line 166: ΔGbinding: the relative binding free energy. Kindly move this sentence to the beginning of the paragraph (line 163) before ΔEMM : the gas phase...

(5) Page 8, line 187: It can be observed that the RMSD value in the TLR8 (apo).....that the TLR8 is more flexible in the presence of agonist. Is it presence or absence, kindly check as higher RMSD is associated with protein flexibility.

(6) Page 10, line 239: On the other hand of the trajectory analysis, we also examine the conformational change of agonists.

Change to "On the other hand of the trajectory analysis, we also examine the conformational changes of the agonists."

(7) Page 11, line 241: The RMSD of agonists in TLR7 and TLR8 complex systems increase from the initial value (0 ns) to a scope of 1.0-3.0 Å and all reach maintain in this scope after 30 ns.

Change to "The RMSD of the agonists in TLR7 and TLR8 complex systems increase from the initial value (0 ns) to a scope of 1.0-3.0 Å and maintain in this scope after 30 ns.

(8) Page 11, line 246: The conformations of the three agonists were superimposed on the imidazoquinoline ring, and the initial frame (0 ns) and a frame near the late stage of the simulation (40 ns) were selected for comparison.

Change to "The conformations of the three agonists were superimposed on the imidazoquinoline ring, the initial frame (0 ns) and a frame....comparison."

(9) Page 11, line 248: The alignment at t= 40 ns of snapshot was chose to represent every conformation in the simulation 30-50 ns, the side chain orientations of the ligands in TLR7 (Fig 5C) and this in TLR8 (Fig 5F) are apparently different.

Change to "The alignment at t= 40 ns of snapshot was chosen to represent every conformation in the simulation 30-50 ns, the side chain orientations of the ligands in TLR7 (Fig 5C) and TLR8 (Fig 5F) are apparently different.

(10) Page 15, line 324: Agonists and representative residues are shown in orange and cyan stick.

Change to "Agonists and representative residues are shown in cyan and orange stick."

(11) Page 15, line 333: In TLR7-H, TLR7-R and TLR7-G systems, the stable hydrogen bonds formed in Asp555 and N and N1 atoms on agonists were maintained during the last 20 ns of the conformations.

Change to "In TLR7-H, TLR7-R and TLR7-G systems, the stable hydrogen bonds formed between Asp555 and N and N1 atoms of the agonists were maintained during the last 20 ns of the conformations."

(12) Page 15, line 334: In TLR7-H and TLR7-R systems, the occupancy of hydrogen bonds between Thr586 and N1 atoms on agonists (TLR7-H: 34.88%, TLR7-R: 43.83%) are similar.

Change to "In TLR7-H and TLR7-R systems, the occupancy of hydrogen bonds between Thr586 and N1 atoms of the agonists (TLR7-H: 34.88%, TLR7-R: 43.83%) were similar."

(13) Page 15, line 336: Nevertheless, the occupancy of the hydrogen bond between N1 and Thr586 of TLR7-G is very low (TLR7-G: 12.99%).

Change to "Nevertheless, the occupancy of the hydrogen bond between N1 and Thr586 of TLR7-G was very low (TLR7-G: 12.99%)."

(14) Page 15, line 338: Since N4 on G was protonated, a hydrogen bond formed between Thr586 of TLR7-G and N4, resulting in a change the orientation of the side chain and the loss of the hydrogen bond between Thr586 of TLR7-G and N2.

Change to "Since N4 on G was protonated, a hydrogen bond formed between Thr586 of TLR7-G and N4, resulting in a change in the orientation of the side chain and loss of the hydrogen bond between Thr586 of TLR7-G and N2."

(15) Page 15, line 342: The details of the hydrogen bonds... formed in Asp543 and N and N1 atoms on agonists were maintained during the last 20 ns of the conformations.

Change to "The details of the hydrogen bonds... formed between Asp543 and N and N1 atoms of the agonists were maintained during the last 20 ns of the conformations.

(16) Page 16, line 344: The occupancy of hydrogen bond between N1 and Thr574 of TLR8-G is zero. Do you mean no hydrogen bond formed between Thr574 and N1 of G in TLR8-G complex?

Response: Yes. There is no hydrogen bond formed between Thr574 and N1 of G in TLR8-G complex system.

(17) Page 16, line 359: are located around the aliphatic tail. And as shown in Fig 8B, the....TLR8-G.

Change to "are located around the aliphatic tail. As shown in Fig 8B, the....TLR8-G."

(18) Page 16, line 367: complex system is zero. And the occupancy...respectively.

Change to "complex system is zero. Further, the occupancy...respectively.

(19) Page 17, line 388: MD simulations and MM-GBSA method were used to model the overall conformational changes and calculate the binding free energies between three agonists and the TLR7 and TLR8.

Change to "MD simulations and MM-GBSA method were used to model the overall conformational changes and calculate the binding free energies between three agonists and the receptors TLR7 and TLR8."

(20) Page 17, line 389: Trajectory analysis showed that TLR7-R and TLR7-G systems form more “open” conformations during the simulation, however, other systems kept in closed conformations.

Change to "Trajectory analysis showed that TLR7-R and TLR7-G systems formed open conformations during the simulation, however, other systems kept in closed conformations.

(21) Page 17, line 392-393: Plus, the calculated binding free energies indicated that three agonists are more sensitive for TLR7 than TLR8, and the rank of the binding free energy values are in agreement with the experimental EC50 values in the cellular assay.

Change to "Moreover, the calculated binding free energies indicated that the three agonists are more sensitive for TLR7 than TLR8, and the rank of the binding free energy values are in agreement with the experimental EC50 values in the cellular assay.

(22) Page 18, line 397: The side chain nitrogen of G is positively charged in an acidic environment, leading to its much less favorite interaction with Ile585 of TLR7 and Val573 of TLR8.

Change to "The side chain nitrogen of G is positively charged in an acidic environment, leading to less favorite interaction with Ile585 of TLR7 and Val573 of TLR8.

Reviewer #2: The present work describes routine but useful simulations to understand the selectivity of agonists towards TLR7 and TLR8. The study was guided by binding free energy analysis, hydrogen bond analyses, decomposition of free energy. All the data have been presented either in main or supplementary text. The results obtained were almost aligned with the experimental results and helpful in determining the selectivity mechanisms. However the manuscript still requires few explanations and language corrections.

We are thankful to the reviewer’s positive and valuable suggestions. We have incorporated all the suggestions of the reviewer that helped to improve the manuscript. The revisions were highlighted in red in the revised manuscript. We hope that you will find the revised manuscript suitable for the publication in PLOS ONE.

1. Suggestions for clarification or modifications:

Response: We greatly appreciate the reviewer’s valuable suggestions. We have double-checked our manuscript and revised it.

(1) Line 13: replace “Highly” with “high”

(2) Line 25: Replace “experiment” with “experimental”

(3) Line 28: Replace “favor” with “favorable”

(4) Line 29: Replace “agonist in” with “agonist of”

(5) Line 29: Replace “imidazoquinolines” with “imidazoquinoline derivatives”

(6) Line 31: Replace “selective” with “selectivity”

(7) Line 32: Insights for “designing” more potent

(8) Line 36: Molecular dynamics simulations

(9) Line 50: both are located

(10) Line 57: diseases and can be used as good adjutants of vaccines

(11) Line 62: “many distributional and functional”

(12) Line 71: “selectivity mechanism of small molecule agonists of TLR7 and TLR8d TLR8”

(13) Line 80: pH

(14) Line 89: replace “structural” to “structure”

(15) Line 100: replace “against” with “for”

(16) Line 111: replace “first” with initially

(17) Line 112: remove “respectively”. Replace second with “thereafter”

(18) Line 116: replace “selectivity mechanism” with “selectivity”

(19) Line 128-29: The TLR7-H system was modeled based on coordinates of TLR7-R by replacing the oxygen to carbon atom. (Mention the aliphatic chain)

(20) Section 2.2 change heading to “Molecular dynamics simulations”

(21) Line 148: replace Cl-1 to Cl-

(22) Line 149: “steepest descent algorithm”

(23) Line 154: 50ns is not a long-time simulation. Just mention 50 ns simulation.

(24) Section 2.3 equations must be included for “∆EMM : the gas phase energy consisting of electrostatic (∆Eelec) and van der Waals (∆EvdW) terms. ∆Gsolv : solvation energy including both the polar solvation energy, ∆Gpolar and the nonpolar solvation component, ∆Gsurf”

(25) Line 192: “confirmational”

(26) Line 197: capitalize “In”

(27) Line 201: remove “are both”.

(28) Line 202: “The distance changes in TLR7-H, TLR8-R, TLR8-H and TLR8-G systems are not evident (Fig 3C and Fig 3D). In apo systems, the distance changes in TLR7 and TLR8 systems are also not evident (Fig 3C and Fig 3D).”

This point should be made clearer by showing the distance (angstrom) in Fig 3A and 3B. or by superimposing TLR7 and TLR8 with highlighted anchor residue distances.

Response: We are thankful to the reviewer for this valuable suggestion. We showed the distance in Fig 3 C and Fig 3D. We added the values of initial distance and average distance during the last 20 ns for six complex systems and two apo systems in section 3.1 as following:

The distance differences between the structures during the last 20 ns MD simulations and initial structure are 7.4 ± 0.15 Å and 7.1 ± 0.13 Å in TLR7-R and TLR7-G systems, respectively. However, the distance changes in TLR7-H, TLR8-R, TLR8-H and TLR8-G systems are not evident (Figs 3C and 3D), which are 2.8 ± 0.14 Å, 1.8 ± 0.08 Å, 3.0 ± 0.12 Å and 2.3 ± 0.05 Å, respectively. In apo systems, the distance changes in TLR7 and TLR8 systems are also not evident (Figs 3C and 3D), which are 3.6 ± 0.15 Å and 0.6 ± 0.08 Å, respectively.

Fig 3. The distance changes of chain A and chain B for TLR7 and TLR8. (A) The R784 of each chain is defined as anchor point on TLR7. (B) The P773, counterpart of R784 of TLR7, of each chain is defined as anchor point on TLR8. (C) The time dependent distance between chainA-R784 and chainB-R784. (D) The time dependent distance between chainA-P773 and chainB-P773. The Δ_distance refers to the difference of the chain A and chain B distance between the structures during the last 20 ns MD simulations and initial structure. The initial distances (0 ns) are 15.1 Å, 12.4 Å, 16.5 Å and 15.3 Å in TLR7 (apo), TLR7-R, TLR7-H, TLR7-G systems, respectively. The initial distances (0ns) are 11.9 Å, 8.3 Å, 11.9 Å and 14.1 Å in TLR8 (apo), TLR8-R, TLR8-H, TLR8-G systems, respectively. The average distances during the last 20 ns are 18.7 ± 0.15 Å, 19.8 ± 0.15 Å, 19.3 ± 0.14 Å and 22.4 ± 0.13 Å in TLR7 (apo), TLR7-R, TLR7-H, TLR7-G systems, respectively. The average distances during the last 20 ns are 11.3 ± 0.08 Å, 10.1 ± 0.08 Å, 14.9 ± 0.12 Å and 11.8 ± 0.05 Å in TLR8 (apo), TLR8-R, TLR8-H, TLR8-G systems, respectively.

(29) Line 217: replace scope with “range”

(30) Line 235: replace “On the other hand of the” by “In another” & change to “changes”

(31) Line 245: replace chose to “chosen”

(32) Line 246: remove “this is”

(33) Line 256: TLR7 and TLR8

(34) Line 260: while G shows the lowest

(35) Line 270: energy is unfavorable to binding include the values here.

(36) Table 1. Comparison of predicted binding free energies and experimental EC50 values of R, H and G with TLR7 and TLR8

(37) SEM (standard error of means) should be reported for all the BFE values e.g -44.21 ± ?? in all the tables.

(38) Table 2. replace “Item” to “Energetic contributions” All energetic contributions should have uniform representation e.g. ∆Evdw, ∆Eele……

All units are in kcal/mol.

(39) Line 293: “It forms stronger electrostatic interactions with Asp555 and stronger van der Waals interaction with Val381 and Phe351” should be supported with energy values.

Response: Many thanks for the insightful comments. As suggested, we have added the energy values to consolidate this sentence, which is also copied as following:

Compared with R, the carbon atom on the side chain of H is more flexible and hydrophobic than the oxygen atom on the side chain, leading H to be more adapted to pocket environment. H forms stronger binding free energies with D555, V381 and F351 (-6.768, -1.450 and -3.503 kcal/mol) than R (-4.542, -0.964 and -2.849 kcal/mol).

(40) Line 306-09: these are repeated lines. Remove.

(41) Line 310: “The residue of TLR8 alignment to Leu557 of TLR7 is Asp545, which might result in a weaker interaction with agonists.” rewrite to make it more understandable.

Response: We are thankful to the reviewer for this valuable suggestion. The sentence has been rewritten in section 3.3 as following:

Agonists R, H and G form C-H-π interaction with L557 of TLR7, which has been also reported in other study for R [64]. However, in TLR8 the corresponding residue of L557 is D545, which could not form strong interactions with agonists due to its shorter side chain.

(42) Line 327: Similar codes (either 3 letter or single) for amino acids should be used in text and figures.

(43) Line 332: side chain and hence the

(44) caption Fig 7. “the important atoms”

(45) Line 364: two or one TLR7 please specify

Response: We are thankful to the reviewer for this valuable suggestion. The sentence has been rewritten as following:

In addition, the occupancy of residues E352, L556 and T586 in two or one TLR7 complex systems is higher than I349, F544 and T574 in TLR8 complex systems. For example, the occupancy of residues E352 in TLR7-H complex system and L556 in TLR7-R and TLR7-G complex systems and T586 in TLR7-G system is higher than I349, F544 and T574 in TLR8 corresponding systems.

(46) Line 382: remove “the” before TLR7 and TLR8.

(47) Line 384: “pockets” to pocket

(48) Line 392: less favorable interactions

(49) Line 396: each agonist

(50) Line 394: quiet to “quite”

2. Need to address these questions?

(1) How have you treated long range electrostatic?

Response: We are thankful to the reviewer for this valuable suggestion. Long-range electrostatic interaction was calculated using PME with the cutoff 9Å. And this sentence is added to the section 2.2 as following:

Long-range electrostatic interaction was calculated using PME with the cutoff 9Å.

(2) Any specific reason for using only MM-GBSA (Not MM-PBSA) protocol for binding free energy calculations.

Response: According to the researches of Hou et al1, the accuracy of MM/GBSA is comparable to MM/PBSA, even in the test of the dataset including 1864 crystal structures MM/GBSA is better than MM/PBSA. However, it is known that MMPBSA is more time expensive than MM/GBSA. Therefore, MM/GBSA was applied to predict the protein-ligand binding free energy in this study.

[1] Sun H, Li Y, Tian S, Xu L, Hou T. Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set. Phys Chem Chem Phys. 2014 Aug 21;16(31):16719-29. doi: 10.1039/c4cp01388c.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

11 Mar 2022

3 Apr 2022

Responses to reviewer#1

The authors have addressed all the comments, however, there are additional comments that need to be addressed by authors before the acceptance of the manuscript.

Many thanks for taking your time to review our responses. We have incorporated all the suggestions of the reviewer that helped to improve this manuscript, and highlighted the revisions in red.

1. Page 4, line 77. Kindly remove this part from the paragraph, "In the late 70s, the first case of MD is Bovine pancreatic trypsin inhibitor [26]. Especially in recent years, benefit from the development of modern hardware and force fields, MD simulation has made huge progress. Protein conformation dynamics that is hardly available by current experimental techniques can be observed using MD simulations [27-34]." As it is extra and does not add scientific value to your manuscript.

So the updated paragraph will be "Molecular dynamics (MD) simulation is a very established computational technique to understand the protein structure-function relationship and guide the drug design. MD simulation has played an important role in characterizing receptor-ligand interaction .... with a simulation time of 10-100us [55-57].

Response: Many thanks for the constructive suggestion. This paragraph has been reduced. The revised paragraph in section 1 is following:

Molecular dynamics (MD) simulation is a very established computational technique to understand the protein structure-function relationship and guide the drug design. MD simulation also has played an important role in characterizing receptor-ligand interaction [26-28], providing the guidance of structure-based drug design [29-31] and typical conformations for virtual screening [32-37]. Besides, it helps to reveal the novel binding sites which have not been captured by NMR and X-ray crystallographic analysis, for example, cryptic binding sites in HIV-1 integrase [38-41]. Up to date, MD simulations have been successfully applied to many large systems, such as the complete HIV1 capsid with 64 million atoms up to 100ns [42-45]. MD simulation has been also instrumental on understanding the protein folding and function regulation with a simulation time of 10-100us [46-48].

2. Page 11, line 246. "In another trajectory analysis, we also examine the conformational changes of agonists. The heavy atom RMSD of agonists was calculated in the six TLR-agonist systems (Figs 5B and 5E) with respect to the initial conformation."

Kindly change it to " We also examine the conformational changes of agonists. The heavy atoms' RMSD of agonists was calculated in the six TLR-agonist systems (Figs 5B and 5E) with respect to the initial conformation."

Response: We are thankful to the reviewer’s positive suggestion. The revised paragraph in section 3.2 is following:

We also examined the conformational changes of agonists. The heavy atoms' RMSD of agonists was calculated in the six TLR-agonist systems (Figs 5B and 5E) with respect to the initial conformation.

3. Page 16, line 349. "The occupancy of hydrogen bond between N1 and T574 of TLR8-G is zero."

Kindly remove this sentence since there is no hydrogen bond between N1 and T574 of TLR8-G.

Response: We are thankful to the reviewer’s valuable suggestions. This sentence has been removed.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

12 Apr 2022

Molecular dynamics simulations reveal the selectivity mechanism of structurally similar agonists to TLR7 and TLR8

PONE-D-21-35895R2

Dear Dr. Jinxia Nancy Deng ,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Anand Gaurav

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Dear Author,

The revised manuscript addresses all the comments of the reviewers and the manuscript can be accepted in the current form.

Regards

Anand Guarav, PhD

Reviewers' comments:

14 Apr 2022

PONE-D-21-35895R2

Molecular dynamics simulations reveal the selectivity mechanism of structurally similar agonists to TLR7 and TLR8

Dear Dr. Deng:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Anand Gaurav

Academic Editor

PLOS ONE