Conformational Dynamics of thiM Riboswitch To Understand the Gene

Subscriber access provided by UNIV OF DURHAM

Computational Biochemistry

Conformational Dynamics of thiM Riboswitch to understand Gene Regulation Mechanism using Markov State Modeling and Residual Fluctuation Network Approach Manish kesherwani, Kutumbarao N. H. V., and Devadasan Velmurugan J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00155 • Publication Date (Web): 25 Jun 2018 Downloaded from http://pubs.acs.org on June 26, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Conformational Dynamics of thiM Riboswitch to understand Gene Regulation Mechanism using Markov State Modeling and Residual Fluctuation Network Approach Manish Kesherwani, Kutumbarao N. H. V. and Devadasan Velmurugan* Centre for Advanced Study in Crystallography and Biophysics, University of Madras, Guindy Campus, Chennai-600025 India *[email protected]

Correspondence Address: Prof. D. Velmurugan UGC-BSR Faculty CAS in Crystallography & Biophysics University of Madras Chennai-600025 Email Id: [email protected]

1 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract: Thiamine pyrophosphate (TPP) riboswitch is a cis-regulatory element in the non-coding region of messenger RNA (mRNA). The aptamer domain of TPP riboswitch detects the high abundance of coenzyme thiamine pyrophosphate (TPP) and modulates the gene expression for thiamine synthetic gene. The mechanistic understanding in recognition of TPP in aptamer domain and ligand-induced compactness for folding of expression platform are most important to design novel modulators. To understand the dynamic behavior of TPP riboswitch upon TPP binding, molecular dynamics simulations were performed for 400 ns in both apo and TPP bound forms of thiM riboswitch from E.Coli and analyzed in terms of eRMSD based Markov state modeling and residual fluctuation network. Markov state models show good correlations in transition probability among metastable states from simulated trajectory and generated models. Structural compactness in TPP bound form is observed which is correlated with SAXS experiment. The importance of junction of P4 and P5 is evident during dynamics which correlates with FRET analysis. The dynamic nature of two sensor forearms is due to flexible P1 helix, which is its intrinsic property. Transient state in TPP bound form was observed in Markov state model along with stable states. We believe that this transient state is responsible to assist the in- and out- flux of ligand molecule by creating a solvent channel around the junction region of P4 and P5 and such a structure was anticipated in FRET analysis. The dynamic nature of riboswitch depends on the interaction between residues on distal loops L3 and L5/P3 and junction P4 and P5, J3/2 which stabilize the J2/4. It helps in transfer of allosteric information between J2/4 and P3/L5 tertiary docking region through the active site residues. Understanding such information flow will benefit in highlighting crucial residues in highly dynamic and kinetic systems. We report here the residues and segments in riboswitch which play vital role in providing stability and this can be exploited in designing inhibitors to regulate the functioning of riboswitches.

Keywords: TPP riboswitch, Markov modelling, Residual fluctuation network, eRMSD, Bacterial metabolism regulation


Page 2 of 37

Page 3 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction: Riboswitches are structured elements commonly found in the non-coding region of messenger RNAs of many bacteria, fungi, and plants. Globally, riboswitch consists of two distinct functional domains. The aptamer domain is responsible for cellular metabolites binding at an elevated concentration to stabilize the global fold which promotes folding of expression platform domain, thus transduces small molecule binding into a genetic regulatory signal. The allosteric communication between aptamer and expression platform domain is established due to the formation of the secondary structural switch (such as terminator /antitermination stemloop or ribosome binding site (RBS)) which regulates gene expression by transcription termination, translation inhibition, alternative splicing and RNA stability1–3. There are 24 classes of riboswitches known till now and three-dimensional structures of several aptamerligand complexes of riboswitches such as purines, thiamine pyrophosphate (TPP), S-adenosyl methionine (SAM), flavin mononucleotide (FMN), glycine, lysine and c-di-AMP riboswitch are reported2. Functions of riboswitches are mainly classified based on their charges of sensor ligands, riboswitch size, shape and mode of regulation (transcription vs. translation) as well as structural and binding properties of the aptamer. In terms of structural organization and ligand binding, riboswitches are classified into two types: Type I riboswitches recognize cognate ligands at junction of several helices or at pseudoknots such as purines or glms riboswitches, respectively and induce local rearrangement in binding site while Type II riboswitch uses distinct peripheral structural elements (other than junctions) to induce global rearrangement in binding site in response to cell's metabolic state and genetic outcome4. The riboswitch-mediated regulation originates from the alternative pairing of sequences from aptamer domain and expression platform domain and is highly depends on the concentration of the free ligand. Recent studies has described that the riboswitch is responsive to the cellular conditions using a self-regulated feedback loop mechanism5. The refolding of riboswitch has been well studied and has provided the linkage between the unbiased mappings of conformational space to the kinetics derived from the ligand-induced release of the ON state. Thus, importance of kinetics and thermodynamic properties in metastable states and the role of ligand in inducing the fold and on/off switching in the riboswitch, play a major role in understanding the transcription regulation. T. The structural features that assist in the on and off mechanism of guanine sensing xpt-pbuX are also studied based on riboswitch folding6. These studies also explained a metastable conformation, the antiterminator state in a transcriptional riboswitch, and the influence of ligand binding coupled 3 ACS Paragon Plus Environment


with successive refolding. This has highlighted that the regulation of riboswitch is kinetically driven and folding of RNA has to compete with the rate of transcribing polymerase. TPP is an active form of vitamin B1 (thiamine) and considered as an indispensable coenzyme in many cytosolic and mitochondrial enzymatic reactions. TPP riboswitch belonging to type II classification of riboswitches recognizes TPP in aptamer domain and acts as the gene regulatory element in mRNA to regulate gene expression of thiamine synthesis. TPP riboswitch has been found in all three kingdoms of life. In eukaryotic plant Arabidopsis thaliana, TPP riboswitch is discovered in 3' untranslated intron region of the THIC operon. It modulates an alternate gene splicing event through the thermodynamic regulation mechanism and its kinetic of folding depends upon temperature7,8. In prokaryotic E.Coli bacteria, there are two types of TPP riboswitches identified in 5' untranslated mRNA of thiM and thiC operons and named as thiM and thiC riboswitch, respectively9. The thiM riboswitch regulates translation of thiM operons whereas genetic regulation of thiC riboswitch occurs both at the level of translation and transcription. In E.Coli, riboswitch-mediated gene regulation is of kinetic nature, thiM riboswitch mediates the gene regulation at low TPP concentration (0.2-1.5 µM), whereas at high TPP concentration (~1.5-15 µM), thiC riboswitch mediates the regulation7. Despite their difference in functional mechanism, crystal structures of TPP riboswitch from E.Coli (thiM) and Arabidopsis thaliana (THIC) showed similar structural arrangements such as they contain two parallel helical arms (P2/J3-2/P3/L3 and J2-4/P4/P5/L5) connected to helix (P1) by the three-way junction and form h-shape fold (Figure 1). They also share a common binding conformation where TPP binds in an extended conformation in two main helical sensor arms; first sensor helix arm (P2/J3-2/P3/L3) forms an intercalation pocket for the pyrimidine moiety of TPP. The other sensor helix arm (J2-4/P4/P5/L5) offers a water lined binding pocket for pyrophosphate moiety of TPP which engages divalent metal ion and stabilizes tertiary contacts in J2-4/P2/P4 regions for modulation of gene expression10. The long tertiary interaction between the P3 stem and L5 loop formed region in these two helical sensor arms showed the necessity to stabilize overall fold of riboswitch. The main difference in both thiM and THIC riboswitches is the extended docking region of P3 stem and L5 loop of the two helical sensor arms in THIC riboswitch which forms a stable P3 stem and leads to slower helix arm dynamics. In thiM, the docking region of P3 stem and L5 loop is smaller and show higher dynamic behaviour on helix arm in TPP binding.


Page 4 of 37

Page 5 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Crystal Structure of thiM Riboswitch from E.Coli (PDB ID: 2GDI) Several MD simulations studies have been performed for oligonucleotides, tetraloops, ribozymes, and riboswitches to investigate the stability of intermolecular interactions, structural dynamics such as conformational rearrangement due to ligand binding and inherited flexibility in the active site as well as the physiological concentration of ions11,12. MD simulations in aptamer of add A and pbuE A riboswitch have been performed in ligand bound and free form to analyze conformational changes associated with ligand binding and associated expression platform for understanding allosteric switching mechanism. It has also highlighted the importance of relative stability in helices to determine the order of riboswitch folding and response to expression platform13–16. Dynamic cross-correlation and network analysis from MD simulations of tetrahydrofolate sensing riboswitch and preQ1 riboswitch have shown metabolite binding and an allosteric switching mechanism for gene regulation17,18. Recently, Markov state modeling has become a complementary method to reduce the complexity of biomolecular dynamics into simple kinetics description to establish thermodynamic and kinetics of slow processes in protein and nucleic acid. Markov Models are generated using Markov chain on the discrete partition of sub-space of configuration space and allow to analyze the longest-living (metastable) set of structures and transition rate between them19,20. Recently, Markov state modeling has been used for protein folding pathways analysis of peptides such as PinWW and Fip35 from MD simulations trajectory21,22. Markov state model from single molecule Foster resonance energy transfer (smFRET) data has also been used in several riboswitches to correlate observed fluorescence spectrum to structural dynamics for 5 ACS Paragon Plus Environment


understanding differences in metabolite bound and unbound forms23,24. Recently, RNA folding pathways of tetraloops and hairpin were predicted based on Markov state model of threedimensional fragments which are collected from high-resolution crystal structures of RNAs as well as several long dynamics from small simulation trajectories25,26. The resistance to antibiotics has become a major concern for treating bacterial infection caused by E. Coli and methicillin-resistant S. aureus (MRSA) bacteria. Riboswitches have been identified in microorganism but they are not yet detected in human. These riboswitches endeavor as specific potential drug target to design novel antibiotics. Crystal structures of TPP bound aptamer domain of thiM riboswitch provide only static event of ligand binding10, whereas the dynamic behavior of TPP ligand binding and the conformational changes in aptamer domain which drives folding of adjacent expression platform domain is lacking. Biophysical characterization using FRET analysis of ligand free and bound form of aptamer domain points out the importance of dynamic properties of various sensor arms which play important role in folding of riboswitch and further facilitating gene regulation, but they lack in structural explanation23,27. Here, we had carried out the molecular dynamics simulations of TPP bound and unbound thiM riboswitch in an explicit solvent for 400 ns timescales. The simulated trajectory is used for kinetic clustering on top autocovariance modes obtained from TICA component analysis on the eRMSD metric. Markov state modeling and flux analysis were carried out for identifying metastable states and their transition flux path. Further, intermolecular interactions, dynamical network analysis, nucleic acid parameters and free energy calculation were analyzed in metastable states for understanding allosteric interactions and conformational changes associated with metabolite binding.

Material & Methods: Molecular Dynamics Simulations of TPP bound and apo form of thiM riboswitch: The three-dimensional structures of apo and TPP bound form of thiM riboswitch are parameterized using nucleic acid force field, ff99bsc0 with the OL3 correction for each RNA residues28,29. Forcefield parameter of TPP ligand was generated by optimizing initial geometry at HF/6-31G* level of theory in Gaussian 03 software30. The derived partial atomic charges were fitted using RESP program available in Antechamber module of Amber 12 and bonded and nonbonded parameters were adopted from General Amber Force Field (GAFF) parameter sets

31,32

. The apo form of thiM riboswitch was obtained by deletion of TPP from a cocrystal

complex of TPP bound form (PDB Id: 2GDI). The generation of an apo form of adenine-


Page 6 of 37

Page 7 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sensing riboswitch was also obtained by deleting the ligand from the ligand-bound form in recent MD simulations of ade A riboswitch33. Initial structures of these two forms of TPP riboswitch were solvated in TIP3P water box with a buffer distance of 10 Å and neutralized with counter ions of Na+ or Cl-

34

. The SHAKE algorithm was used to constrain hydrogen

linked atoms in RNA. Systems were used for energy minimization of 5000 cycles using steepest descent followed by conjugate gradient energy minimization techniques. Systems were heated up for 500 ps in NVT ensemble for equilibrating solvent molecule and counterions around riboswitch. Further, systems were equilibrated in NPT ensemble protocol for 2 ns with 2 fs integration time step to bring temperature and pressure at 300K and 1 atm respectively, by Langevin thermostat and barostat35. Particle Mesh Ewald (PME) was used to treat long-range electrostatic interactions and the 10Å nonbonded cutoff was used to treat pairwise interaction. Finally, simulation production was performed in periodic boundary condition for 400 ns with 1 fs integration time to the equilibrated systems. MD simulations were carried out by using CUDA version of PMEMD in GPU cores of NVIDIA Tesla C2075 server built in CAS in Crystallography & Biophysics, University of Madras, Chennai, India. Simulation trajectories were analyzed for structural deviation in terms of eRMSD on overall as well as a residual basis fluctuations. Intermolecular contacts such as hydrogen bonding (D-H...A) and hydrophobic contacts (contacts around 7Å of ligand) were also analyzed between TPP and active site of the riboswitch. To understand the structural features of riboswitch, various nucleic acid parameters, namely, helical rise, propeller, roll, hydrogen bonds in Watson -Crick base pair etc. were analyzed using CPPtraj module of Amber12. For structural visualization, Chimera v1.11.2rc program was used.36 Molecular dynamics simulations of TPP bound in L3 loop mutant TPP riboswitch (PDB Id: 3K0J) and recently identified TPP modulator (inhibitor) (PDB ID: 4NYA) were also carried out for 100 ns.

Markov state modeling of TPP bound and apo form of thiM riboswitch: Markov state model for TPP bound and apo form of thiM riboswitch was built on eRMSD metrics using the pyEMMA software37. eRMSD was calculated for each nucleotide in riboswitch (residual-wise) using baRNA software38. Time-lagged independent component analysis (tICA) was used to find low-dimensional subspace of residual eRMSD data of simulated trajectory. Projection of top two TICA components was subsequently used for conformational clustering with interspace distance cutoff of 0.9 Å using regular space clustering algorithm. Maximum Likelihood Estimation (MLE) algorithm was used for 7 ACS Paragon Plus Environment


Page 8 of 37

generating Bayesian Markov model by estimating transition rates on microstate cluster which are further lumped into macrostate by the PCCA+ algorithm. The built Markov state models were 10 fold cross-validated using Chapman-Kolmogorov test and best generated Markov model were used for calculating flux between metastable states37. To differentiate metastable states and to understand the conformational transition, 1000 sample conformations from each metastable state in apo and TPP bound form are retrieved and analysed using fluctuation correlation network and binding interaction study.

Network model construction of each metastables states in TPP bound and apo form: Network model was constructed for each metastable state from TPP bound and free form. The network contains nodes which are defined by N1 or N7 atom of nucleobase in each nucleotide residue and edges are formed only with those nodes which have persistent intermolecular interactions around 7 Å in 75% trajectory. Edge distance is weighted with strong correlation (above 40%) on eRMSD of each node in sampled conformations of metastable states. The network was analyzed for network centrality such as an eigenvector 39. The network was also used for node partition into the community using Girvan Newman algorithm from a script developed by Luthey-Schulten Lab

40,41

. VMD v1.9.3 program was

used for visualising the residual fluctuation network 42.

Binding free energy calculation of TPP bound and unbound thiM riboswitch: Binding free energies are calculated to predict binding affinity of TPP ligand to thiM TPP riboswitch using MMPBSA approach by MMPBSA.py script available in Amber1243. Here 1000 samples from each metstable states were used to calculate the enthalpy of binding free energy for each molecular species in the formation of TPP bound riboswitch complex. Thus enthalpy includes gas phase energy (electrostatic and van dar waal) from molecular mechanics forcefield terms and solvation phase energy was calculated using PoissonBoltzmann Solvent accessible (PBSA) solver for polar solvation energy and solvent accessible surface area approximation for nonpolar solvation energetics44. ∆Gbinding = ∆Hbinding - T∆Sbinding

(i)

∆Hbinding = ∆Ecomplex - ∆Ereceptor – ∆Eligand

(ii)

∆ETot = ∆Egas + ∆Esolv

(iii)

∆Egas = ∆Eele + ∆Evdw

(iv)


Page 9 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Binding free energies for TPP bound in mutant thiM riboswitch and inhibitor-bound form are also calculated using above approach from the last 20 ns simulated trajectory and compared with wild-type of TPP bound form. Residual decomposition of binding free energy is also carried out for these complexes. Free energies of apo and TPP bound thiM riboswitch were also calculated using the above formulae, except that we remove difference (∆) and calculate for particular species.

Results & Discussion: Molecular dynamics simulations of TPP bound and unbound form of thiM riboswitch: Molecular dynamics simulations are hypothetical single molecule experiments and provide an investigation of unlimited temporal and spatial resolution of biomolecules at the atomic level. Here, thermal sampling of biomolecules using single atomistic conformation is performed which mimics thermal fluctuations of real molecules. Molecular dynamics simulations of TPP bound and unbound form (apo form) of thiM riboswitch are carried out for 400 ns timescale in explicit water to understand dynamical behavior of structural and thermodynamic properties. Structural deviation of apo and TPP bound form of thiM riboswitch with initial crystal structure has been analyzed in terms of eRMSD (Figure 2a). eRMSD is a nucleobase centric deviation metric for discriminating structural and kinetically distinct conformation and has been reported as a better parameter in analysing RNA dynamics when compared to standard RMSD metric38. It measures the deviation in relative position and orientation between two nucleobases and relates the rearrangement in base stacking and base pairing interactions in RNA structure38. Here, eRMSD was clearly able to distinguish distinct conformational states in 400 ns simulation timescale of apo and TPP bound complex, where higher degree of conformational restrictions have been observed in TPP bound form compared to apo form (Figure 2a). The standard backbone RMSD is unable to differentiate the conformations (Figure S1).



Figure 2. (a) Structural deviation of apo and TPP bound thiM riboswitch with respect to eRMSD and simulation time (ns) , (b) Projection of simulated trajectory (kinetic map) into top two tICA (IC) components and (c) Relaxation timescale of riboswitch system using Bayesian Markov model (The error bars are at 95% confidence intervals estimated using reversible transition matrix) Markov state modeling of simulated TPP bound and apo form of thiM riboswitch: In Markov state modeling (MSM), conformational dynamics of biomolecule describes a kinetic network of transitions between metastable states and enables the identification of relevant metastable states and transition rates between them. Here, MSM models are constructed for TPP bound and apo form of simulated trajectories based on residue wise 10 ACS Paragon Plus Environment

Page 10 of 37

Page 11 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


eRMSD distance metrics. For capturing the best timescale dynamics, time-lagged independent component analysis (tICA) is used to linearly transform the high dimensional space of residual eRMSD into low dimensional space reflecting the slowest conformations motions. For performing tICA, optimum lag time was identified based on tICA coordinates obtained at different lag time and has a reduction ability to minimize dimensions and capturing cumulative kinetic variance of more than 75%. We used 10 ns lagtime for tICA analysis and reduced feature dimensions of about 7 from a total of 77 which showed 90% relevant kinetic information for both apo and TPP bound form. Top two tICA components of apo form and TPP bound form showed a cumulative variance of 56.6 and 51%, respectively. Projection of top two tICA components (ID1 & ID2) on the simulated trajectory of riboswitch gives a glimpse of conformational states that exist in folding/binding landscape by computing free energy (Fi=-ln zi) on the histogram of top two tICA components, where zi is bin counts (Figure 2b). In our tICA based free energy landscape, apo form shows four possible conformation states whereas TPP bound form shows five conformation states. Conformations mapped in top two tICA components are clustered into 516 microstate cluster for both apo and TPP bound form with interspace distance cutoff of 0.9 Å using regular space clustering protocol (Figure 2b). Implied timescales are computed based on estimating transition probability matrix on these clustered microstates and generate Bayesian Markov state model at different lagtimes of 1 step to 2000 steps (1step=1ps) with error bar (using bootstrapping approach). Implied timescale computed at these lag time shows a clear gap between the slowest and next two slowest timescale processes in apo form, which suggests the possibility of four conformational states in the simulated apo form. But in TPP bound form, four slowest timescale process were observed, indicating five conformational clusters (Figure 2c). Further, based on implied timescales at different lagtimes, lagtime corresponding to constant implied timescales such as 500 ps and 1000 ps are considered as markovian timescale for apo and TPP bound form, and are used to estimate an final Bayesian Markov model by Maximum Likelihood Estimation (MLE) algorithm.Relaxation timescale differences calculated on these Bayesian Markov model also show three gradual steep relaxation processes (four metastable states) in apo form and four gradual steep relaxation processes (five metastable states) in TPP bound form (Figure 3). Further, based on these relaxation timescales, the PCCA+ algorithm was used to group microstates into macrostate for apo and TPP bound which yields four and five states, respectively (Figure 3). Markov state model was 10 fold cross-validated with 95% confidence 11 ACS Paragon Plus Environment


level by conducting the Chapman-Kolmogorov test (Figure S2). Markov models in both apo and TPP bound form of riboswitch correlate well in transition probability among metastable states in estimation from simulation trajectory and calculation from the model (Figure S2).

Figure 3. Distribution of transition probability of metastables on top two tICA component of apo and TPP bound thiM riboswitch


Page 12 of 37

Page 13 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To understand the transition path for observed metastable states in MSM models of both apo and TPP bound form, flux analysis was performed. Flux analysis describes pathway of conformational transition in low populated macrostates to high populated macrostates based on forward committor probability on reversible Markov state process. Flux analysis of the metastable states in apo form is presented in Figure 4 where the disk represents discrete states (with area proportional to the stationary probability of state) and the arrow represents transition flux (µs) between different states. Based on the stationary probability of each metastable state in apo form, state 1 depicts rare state or low populated state and state 4 is highly populated and stable state whereas states 2 and 3 are intermediate states. The optimal pathway of conformational transition of low populated (state 1) to high populated state (state 4) is in the evolving of transition probability of metastable states such as State 1 --> 2 --> 3 -->4. Thus, flux analysis in apo form suggests that state 4 is more probable state whereas other states are either in unfolded conformation in apo form or transient form of TPP bound. On the other hand, flux analysis of MSM derived five conformational states in TPP bound form showed that macrostates 4 and 5 have a high stationary probability which suggests most probable states, whereas macrostates 2 and 3 are in the intermediate states. State 1 is the least probability of occurrence in TPP bound form. Figure 4 describes the path transition between two probable metastable states, state 4 to state 5. The more probable transition path is covered through state 3 whereas other paths through states 1 and 2 have high transition rate (flux). Thus transition path suggests that state 3 acts as a transitional state, where state 4 transforms to state 3 which in turn transforms to state 5 with equal transition rate (~60 ns). There is also direct path between state 4 and 5 but their transition time is too high (~9000 ns) and riboswitch will rarely visit this path to interconvert due to suboptimal energetic cost. The kinITC experiment has also shown that in thiM riboswitch, TPP binding and their gene regulation function exist in two-step kinetic mechanism. In in the first step TPP binds and in the next step TPP binding triggered the RNA folding to form stable conformation7. Thus, based on TPP triggered conformational changes, two types of states are categorized in TPP bound form; strongly bound and weakly bound states. The transition flux analysis in TPP bound form also showed that states 4 and 5 are two prominent states which exist in TPP bound form and are more stable and can interchange with a small perturbation in energy level through the transient state (state 3).



Page 14 of 37

Figure 4. Flux path analysis in metastable states of Apo form and TPP bound form (The circle has an area proportional to equilibrium probability and arrow indicate the transition flux

(in

µs) for distinct transitions between different metastable states) Analysis of metastables states in apo and TPP bound thiM riboswitch: The metastable states observed from Markov state modeling in apo and TPP bound forms are used to calculate free energy in Poisson Boltzmann (PB) implicit solvent using MMPBSA approach. Here, free energy depicts only enthalpy energetic contribution and entropy is not considered. It was also reported that TPP binding in thiM riboswitch is an enthalpy-driven process and the entropic effect is less effective in free energy calculation45. In apo form, free energy in metastable states shows that state 4 is also highly favorable and stable compared to the other states (Table 1). Similarly, in TPP bound form, free energy of state 4 is the most favorable and states 3 and 5 are intermediate favorable states compared to others states. State 1 and 2 have least favorable energy and suggest that these states are unfolded or less prominent states in TPP binding. Thus, free energy complement the stationary probability observed in flux analysis. As state 4 of TPP bound and state 4 of apo form are more favorable, radius of gyration calculation in these states shows that TPP bound form (37.76+/-0.26 Å) has more structural compactness compared to apo form (Table 1). Reported SAXS analysis in TPP bound and unbound form of riboswitch also suggested considerable compactness upon metabolite binding. Radius of gyration (RGYR) in states 3 and 5 of TPP bound form is higher


Page 15 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


compared to other states and it suggests that these states have slight expanded conformation compared to the stable state 4. Table 1. Free energy and Radius of gyration in metastable states in apo and TPP bound thiM riboswitch (Bold value indicate more favorable energy correspond to metastable state of apo and TPP bound form)

RGYR

Free Energy

Binding

(Å)

(kcal/mol)

free energy (kcal/mol)

Apo

State 1

form State 2

State 3

State 4

TPP

State 1

Bound form

State 2

State 3

State 4

State 5

37.53

-15694.91

-

(+/-0.28)

(+/-38.62)

37.43

-15704.55

(+/-0.32)

(+/-36.80)

38.24

-15714.14

(+/-0.42)

(+/-39.42)

38.00

-15709.54

(+/-0.32)

(+/-37.14)

37.66

-13304.11

-18.06

(+/-0.29)

(+/-37.94)

(+/-13.09)

37.83

-13307.83

-21.51

(+/-0.29)

(+/-35.37)

(+/-12.94)

38.59

-13312.02

-39.15

(+/-0.18)

(+/-37.39)

(+/-8.78)

37.76

-13321.00

-37.47

(+/-0.26)

(+/-37.47)

(+/-8.89)

38.57

-13311.70

-35.53

(+/-0.16)

(+/-36.99)

(+/-10.50)

-

-

-

Combined force spectroscopy and single molecule Foster Resonance Energy Transfer (smFRET) studies of TPP bound form had shown that upon ligand binding, P3 and L5 arms undergo rapid conformational changes and this long-range tertiary interaction stabilizes the global fold of aptamer domain and provides the specificity for TPP binding in weakly bound state F' and strongly bound state F''27. Distance distribution of nucleotides between P3-L5 arm 15 ACS Paragon Plus Environment


was analyzed to understand the binding interaction between P3 and L5 (Figure 5A). Results showed that all four metastable states in apo form showed single peak suggesting no conformational changes in P3-L5 contacts whereas in TPP bound form, two peaks were observed (Figure 5A). States 3, 4 and 5 have strong ability to form a tertiary interaction in nucleotide pair C24-A69 in P3-L5 region compared to states 1 and 2. Considering the overall docking region of P3 and L5, state 4 has compact binding in docking region comprising residue pair C21-G72, C22-U71, C23-A70, C24-A69 compared to states 3 and 5, due to the higher density with strong residual contact pair of C21-G72 which is at the proximal of the binding site (Figure 5A). Conformational changes in G72 were also observed in the co-crystal structure of the inhibitor-bound form of thiM complex (PDB ID:4NYA) where the inhibitor does not have direct contact with pyrophosphate binding pocket but inhibitor-induced conformational changes in G72 are observed to regulate the gene regulatory activity of riboswitch46. From the distance distribution analysis, we can conclude that state 4 is more stable state whereas states 3 and 5 are intermediate states (open state) and states 1 and 2 are weakly bound ones. Thus the compactness of P3-L5 arm increases the binding affinity of TPP and the flexibility in L5-P3 arms are shown to be essential for proper orientation of the helix arms, which prime the arms to dock and the aptamer to adopt the tight binding conformation. It was also reported that the dynamic nature of P1 stem is essential for switching and ligand recognition mechanism and it has also been reported that P1 is thermodynamically a weak stem compared to other P2-P5 stems. The conformational changes in the distal region such as binding pocket and P3/L5 tertiary interaction (which is formed by nucleotide A69 of L5 which stacks between neighboring residue A70 and C24 of sensor arm P3) stabilize the P1 stem. In state 4 of TPP bound form, distance density plot for A12-U87 shows weak contact (19 Å) compared to other states and other residual pairs in P1 stem show almost equal density distribution contacts in all states (Figure 5B). In apo form, residue pair A12-U87 has strong contact (15 Å), thus observation of weak contact in state 4 of TPP bound form shows the dynamic nature of P1 stem which is important for information transfer to binding region and P3/L5 distal tertiary interaction region (Figure 5B). In case of state 3, strong contact in A12U87 as well as more density distribution in another residual pair of the P1 stem were observed suggesting transition state between apo and TPP bound form. Analysis of eRMSD based residual fluctuation in metastables states of TPP bound form shows that P4 and P5 stems have more restricted residual fluctuations compared to apo form (Figure S3). Among states 3, 4 and 5 of TPP bound form, state 4 has more restricted fluctuation in the P1 stem, L3, and L5 loop compared to states 3 and 5. It was also noted that fluctuation in the P1 stem is similar in both 16 ACS Paragon Plus Environment

Page 16 of 37

Page 17 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


states 3 and 5 and conformational restriction in L3 loop is more in state 3 compared to state 5. P3 stem and adjacent loop L3 have also shown to stabilize distal tertiary interactions between P3 and L5 loop in previously reported mutation, ITC, and crystal structure studies of L3 mutant TPP bound form45.


Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

Figure 5A. Distance distribution of residues between P3 helix and L5 loop in apo and TPP bound form

18

ACS Paragon Plus Environment

Page 18 of 37

Page 19 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5B. Distance distribution in P1 helix of apo and TPP bound form



Further, comparative analysis of states 3 and 4 in TPP bound form and state 4 in apo form suggests that more conformational restriction exists in TPP binding region such as J3/2 region and P4-P5 junction (G60,-U62, G76-U79) in addition to P1 stem and L3 loop in state 4 of TPP bound form compared to apo form (Figure S4). Residues in P3 and L5 are also more restricted in state 4 of TPP compared to state 3. Eigenvalue calculation in correlated residual fluctuation network defined by the node connected with other important nodes also suggests that in TPP bound form (State 3 and 4), strong node connectivity in the P3 stem and the L5 loop was observed compared to apo form (Figure S4). Less structural deviation in L3 loop in TPP bound form compared to apo form is due to high eigenvector centrality in this region (Figure S4). Betweenness centrality in these networks suggests that junction of P4 and P5 stem have many nodes where shortest paths are passed in the overall network (Figure S4). It suggests that residues from P4 and P5 stems are more important in either RNA folding or TPP binding.

Dynamical residual fluctuation Network analysis in metastable states of apo and TPP bound form: To understand P1 stem dynamics and the allosteric effect in P3/L5 docking region, dynamic contact map of nucleobases in each metastable state of apo and TPP bound riboswitch was filtered with residual cross-correlation of eRMSD and the residual fluctuation network was constructed. Network community was further analyzed using Girvan-Newman algorithm on residual fluctuation network and mapped with secondary structural elements of the riboswitch. Intercommunity connection strength is measured in terms of the sum of betweenness of edges connecting communities to indicate the importance of edge for intercommunication. Nodes belonging to the same community are highly interconnected and can communicate through multiple routes in a network with a small difference in distance (energy). Intercommunity connections between communities are mostly formed with a small number of node connections and indicate a bottleneck for information transfer in the network. Residual fluctuation networks are partitioned into seven major communities, C1 contains P1 stem, C2 contains P2, J2/4, P4, C3 contain junction region of P4 and P5 stems, C4 community contains P5 stem, C5 contains P3 and L5 loop, C6 contains P2 stem and J3/2 junction and C7 contains P3 stem. A detailed description regarding the composition of each community in each metastable state is a given in Table S1. Bottleneck residual coupling in


Page 20 of 37

Page 21 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


intercommunity communications are provided in Table S2 for each metastable states of apo and TPP bound form. Correlated residual fluctuation network in all the four metastables states of apo form has the same node assignment as well as interconnectivity in communities except in state 2, where P2 stem and J3/2 junction forming C6 community are broken into J3/2 forming community NN and P2 stem forming community and specifies the important role of G21A43 in proper folding of aptamer (Figure 6A). The difference in the strength of intercommunication between communities in residual fluctuation network of metastable states was also observed. Intercommunication strength of communities in state 3 and 4 shows almost the same except the strength between C1 community (formed by the P1 stem) and C2 community (formed by J2/4, P2, and P4 ) (Figure 6A). The reason is the difference in bottleneck residual coupling pair G51-A85 in state 3 and G51-A84 in state 4 (Table S2). Residual coupling of G51-A84 was also identified in TPP bound form. Thus, due to strong residual coupling formed by this residue pair in state 4, this state can be considered as a more probable state in apo form. In state 1, community C2 (J2/4 junction, P2 and P4 stem) and community C3 (junction of P4 and P5 stem) have direct connection which lead to change in communication strength observed in the overall community network and makes state 1 as an inactive state for apo form or TPP competent binding form (Figure 6A). In state 4 of apo form, there is no communication between C2 and C3 and it suggests that there is no information flow from J2/4, P2, P4 (C2) to P3/L5 (C5) through the junction of P4/P5 (C3). Previously reported FRET analysis on TPP bound and apo form revealed the importance of P3/L5 interaction, the formation of the stable J2/4 junction and dynamic nature of P1 stem in stabilizing the aptamer-TPP complex which then promotes folding of expression platform by forming an alternative structure in the P1 stem and regulate gene expression23,27. The above transition probability based Markov state model generated for apo form shows state 4 to be a more stable state and state 1 as less stable or TPP binding competent form. The comparative network analyses in these metastable states also support the same. Bottleneck forming residues for communication in a community are almost the same in all four metastable states with slight sequential change (Table S2). Thus, communication strength changes in community of nucleotide fluctuation network leading to a change in conformational shift in thermodynamic stability of apo form.



Figure 6A. Classification of residual communities (Numbered as 1,2,3....) in dynamical network of metastable states in Apo form. Brown color code residues represent bottleneck residues which occur in the majority of shortest path connecting nodes between different communities. For better visualization of the network, intercommunication between communities is replotted in macrolevel nodes where nodes represent the total number of 22 ACS Paragon Plus Environment

Page 22 of 37

Page 23 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


residues belonging to a particular community and mapped with secondary structural elements. The intercommunication strength, represented by the width of an edge in this network, is scaled and labeled in values in terms of the number of a shortest path passing through communities.

In TPP bound form, intercommunication network between communities is almost the same in states 1, 2 and 4 whereas states 3 and 5 have slight changes. In TPP bound form, transition probability in metastable states and free energy suggests states 3, 4 and 5 to be more prominent and states 1 and 2 to be the least. Residual correlated fluctuation network in TPP bound form shows that state 4 has more interconnected nodes in community C1 (P1 stem) and C5 (P3 and L5) compared to states 3 and 5 (Figure 6B). In addition, state 4 also shows strong communication between community C5 (P3 and L5) and community C6 (J3/2 region and P2 stem) by residual pair G40-A43, where stacking interactions are observed with HMP ring of TPP. Community C5 (P3 and L5 region) also has strong intercommunication with C7 (P3 stem). There is also strong information flow between P1 stem and J2/4 junction by G51-A84 residual pair (Figure 6B and Table S2). Due to strong intracommunity in P3 and L5 region and P1 stem, state 4 shows more stable community and strong communication strength in both sensor arms compared to other states (including state 3 and 5).



Figure 6B. Dynamical network of residual communities in metastable states in TPP bound form


Page 24 of 37

Page 25 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Comparison of dynamical fluctuation network in states 3, 4 and 5 shows that there is a difference in the connectivity in contacts between the junction of P4 and P5 stem (C3) and P5 stem (C4) which lead to the formation of new NN community in states 3 and 5 by P4 stem between C3 (junction of P4 and P5) and C4 (P5) communities. It suggests that junction of P4 and P5 have conformational changes in states 3 and 5 whereas state 4 does not have or has less conformational changes and it does not lead to change the residual connectivity in P4 and P5 forearms. In state 4, there is a direct communication between C3 and C4 by A61-C77. In state 3, indirect information flows by residual pair U59-A80 of C3-NN and U64-A75 of NNC4 communities (Figure 6B and Table S2). It was also reported that dynamics in junction of P4 and P5 stem (which is at the base of binding pocket) provide the signal to propagate in P1 stem and P3/L5 forearm and leads the aptamer sensor arms to relatively open Y-shaped conformation. Due to conformational changes in the junction of P4 and P5 stem of states 3 and 5, allosteric information transfer between J2/4 junction region and P3/L5 forearm is disrupted or longer path is covered by a correlated residual pair of C2-C3, C3-NN, NN-C4 and C4-C5 communities compared to state 4 (C2-C3, C3-C4, and C4-C5) (Figure 6B). The shortest path is also calculated to understand the stability between J2/4 (A84 stabilizes J2/4) and L5 forearm (A69 stabilizes P3/L5 docking region). In the shortest path, shorter distance between these nodes indicates larger correlations of nodes along the path and the greater allosteric signal between J2/4 to P3/L5. Thus state 4 has a high allosteric signal between J2/4 (A84) to P3/L5 (A69) region which is governed by most of the known conserved active site residues, basically G17, G18, G19, G21, A43 and A70 and most of them are also identified as bottleneck residues in communication in community analysis. States 3 and 5 have weak allosteric signal due to more distance and their covering path between A84 and A69 is different and it is almost similar to covering path in apo form (state 4) (Table S3). The allosteric signal propensity in state 4 was also significant in shortest path between bottle neck residue observed in community communication (A84-U71) (Table S3). It was also observed that state 3 has a more intracommunity connection in C6 (J3/2 and P2 stem) where HMP ring of TPP interacts with J3/2 compared to state 4. Residual fluctuation network in apo form has also shown that J3/2 and P2 stem forming C6 community has more interconnected subnetwork, thus states 3 and 5 are intermediate states and state 4 is a strongly bound state. Previously, OH footprinting, chemical probing and crystallographic experiment for TPP binding with thiM riboswitch have also suggested that initial TPP binding occurs through strong interactions between pyrimidine ring of TPP with the conserved loop sequence (J3/2) between P2 and P3. Further, ligand binding is stabilized through the formation of a G60 25 ACS Paragon Plus Environment


contact in pyrophosphate binding pocket and protection of C14 and C26 in P3 by the closing of the apical loop L5. Thus, the experimental report of TPP binding complements our dynamical fluctuation network analysis and transition probability such that state 3 is a transition state where TPP binds with less proximal P3-L5 interaction whereas state 4 forms more stable interaction and is considered as a strong binding state. Based on information transfer in these networks, we can suggest the importance of junction regions, namely, J2/4 junction, junction of P4 and P5 stem and J3/2 junction in TPP binding whereas in the apo form, information transfer in this region is less. Binding interaction analysis in TPP binding: Binding free energy calculation in TPP bound form of riboswitch shows that state 3 (39.15+/-8.78) has more thermodynamically favorable binding than state 4 (-37.47+/-8.89) and state 5 (-35.53+/-10.50). States 1 and 2 have the least binding free energy in TPP binding (Table 1). The crystal structure of TPP bound riboswitch showed hydrogen bond interaction between G19 and G40 with HMP ring and G78 with terminal phosphate of TPP. Hydrogen bonding interaction statistics in all five metastable states obtained in Markov model of TPP bound form (including states 3, 4 and 5 which are energetically favorable) showed G19 to have an equally favorable tendency to form hydrogen bond interaction with pyrimidine moiety of ligand (Table 2). In state 1 and 2, the hydrogen bond formed between G78 of junction P4/5 and pyrophosphate moiety of TPP has less occupancy while states 3, 4 and 5 have high occupancy (Table 2). In a similar way, analysis of non bonded interactions between HMP ring of TPP and J3/2 region (also called as T loop like turn) of TPP bound complex showed that states 3, 4 and 5 have almost equal contacts occupancy whereas states 1 and 2 showed less occupancy (Figure S5). On the other hand, among states 3,4 and 5, state 3 has more non-bonded contact occupancy in A56, C57 of P4 stem and C74, G76, C77 and A80 of the P5 stem which plays role in interaction with pyrophosphate group of TPP. In spite of these differences in occupancy of non-bonded contacts with TPP binding, in all metastable states, coordination of Mg2+ ions is consistent with pyrophosphate group of TPP and junction of P4 and P5 stem (G60, G78, U59, C77, A61) (Figure S6 and Table S4). The buried surface area was identified in residues U39, G42, A43, of J3/2, G51, A56, C57, C73 of J2/4 and A80 of the junction of P4 and P5 in TPP bound form compared to apo form and suggests the importance in folding/ligand binding (Figure S7). The buried surface area in these residues is higher in states 3 and 5 compared to state 4 which directly correlates with the ability to form contacts. Thus binding analysis showed that transient state (state 3) has more 26 ACS Paragon Plus Environment

Page 26 of 37

Page 27 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


favorable binding conformation. Stable state 4 is more favorable compared to state 5 and it has slightly less interaction in the junction of P4/P5 compared to transient state (state 3). States 1 and 2 in TPP bound are weakly bound or unstable states. Binding interaction of TPP with thiM riboswitch reveals that pyrophosphate moiety is more important to form a more stable complex. Table 2. Hydrogen bond interaction statistics in metastable state of TPP binding in thiM riboswitch, Inhibitor bound and mutant TPP bound complex

Hydrogen bond interaction (occupancy in percentage %) TPP bound form (Wild Type)

TPP bound form (Mutant Type) Inhibitor bound form

State 1 G19(0.87), G40(0.52), G78(0.19)

State 2

State 3

State 4

G19(0.88), G19(0.90), G19(0.82), G40(0.45), G40(0.37), G40(0.28), G78(0.26) G78(0.50) G78(0.41) G19(0.93),G78(0.37), G40(0.32)

State 5 G19(0.90), G40(0.30), G78(0.52)

U74(0.98), G40(0.41)

To further validate the importance of P3-L5 docking interactions in TPP binding and their induced fit effect in P4/P5 junction of pyrophosphate sensor domain, binding of TPP in the best three metastable states (state 3, 4 and 5) are compared with the simulations of mutant TPP bound structure which does not have any perturbation in TPP binding core region, but L3 loop is mutated 45. Structural deviation using eRMSD suggests that mutant form has more fluctuation compared to a wild form of TPP bound riboswitch (Figure S8). The binding free energy of TPP binding in mutant complex is less favourable when compared to wild-type (Table 3), this difference in binding free energy can also be observed in decomposed binding free energy in mutant complex, where over stabilization of A61, G73, C77 in P5 helix had incurred due to the binding of TPP (Figure 7). Such effect can also be observed in a decrement of hydrophobic contact occupancy of A61, G72, C73 and C74 residues of the P5 helix (Figure 7). No significant difference in hydrogen bonding and hydrophobic interaction occupancy was observed between the mutant and wild complex. From the simulations and the network clustering analysis, we anticipate that the difference in helix arm orientation of P3 27 ACS Paragon Plus Environment


Page 28 of 37

and P5 which brings P3 to close adjacent loop L5 plays a crucial role in stabilization of overall TPP riboswitch. In addition, from the analysis of known inhibitor bound TPP riboswitch, it is observed that the inhibitor has a stronger binding affinity in term of binding free energy and more restricted structural deviation (eRMSD) compared to that of wild-type and mutant TPP bound complex (Figure S8). This can be interpreted from the observation that there is an indirect stabilization in A56, C57, G60, A61 and C74 of P5 helix in pyrophosphate binding cavity and strong hydrophobic contact occupancy and their favorable binding energetic with G40, G42 and A43 in aminopyrimidine binding pocket in pyrimidine sensor arm. Inhibitor also forms strong hydrogen bonding interaction with U74 in addition to G40 where this interaction is absent in TPP bound. Table 3. Binding free energy calculation of kinetic clusters in TPP bound thiM riboswitch

VDWAALS EEL EPB ENPolar DELTA Ggas DELTA Gsolv DELTA TOTAL

TPP bound @State 3 -27.59 (+/- 4.93) 700.73 (+/-24.31) -708.13 (+/-21.76) -4.16 (+/-0.07) 673.14 (+/-22.99) -712.30 (+/-21.77) -39.15 (+/-8.78)

TPP bound @State 4 -26.30 (+/-4.75) 678.96 (+/-24.06) -686.03 (+/-18.94) -4.09 (+/-0.07) 652.66 (+/-22.98) -690.13 (+/-18.93) -37.47 (+/-8.89)

Mutant TPP bound form -21.41 (+/-4.94) 587.93 (+/-20.34) -598.91 (+/-23.73) -4.07 (+/-0.11) 566.51 (+/-18.42) -602.98 (+/-23.71) -36.46 (+/-13.39)


Inhibitor bound form -29.83 (+/-2.83) -1417.85 (+/-13.76) 1402.08 (+/-11.92) -2.28 (+/-0.04) -1447.69 (+/-13.54) 1399.79 (+/-11.92) -47.89 (+/-5.68)

Page 29 of 37

Hydrophobic contact occupancy

1.2 1

0.8 0.6 0.4 0.2 0 RG RU RG RU RG RA RG RA RA RA RC RG RA RG RC RC RG RC RA 19 20 21 39 40 41 42 43 45 56 57 60 61 72 73 74 76 77 80 TPP bound@state 3 TPP bound@state 4 Mutants@TPP bounds Inhibitor bound form 8 6

∆GTot (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4 2 0 -2

RG RU RG RU RG RA RG RA RA RA RC RG RA RG RC RC RG RC RA 19 20 21 39 40 41 42 43 45 56 57 60 61 72 73 74 76 77 80

-4 -6

TPP bound@state 3 Mutants@TPP bounds

TPP bound@state 4 Inhibitor bound form

Figure 7. Hydrophobic contact occupancy and residual binding free energy decomposition of active site in inhibitor bound, TPP binding in wild and mutant form

Conclusion: Resistance to antibiotics has become a major concern for treating bacterial infection and identification of TPP riboswitches in microorganism endeavor them as specific potential drug target to design novel antibiotics. In this study, conformational dynamics of TPP bound and apo form of thiM (TPP) riboswitch have been described through integrated MD simulations, Markov state modeling and residual fluctuation network on nucleobase 29 ACS Paragon Plus Environment


centric deviation (eRMSD) metric. Transition probability based flux analysis in Markov state model shows the one most probable state (state 4) in apo form whereas in the TPP bound, high stationary probability was observed in state 4 and 5 and state 3 is intermediate. Identification of weak and strong TPP bound metastable states was also reported in ITC experiments. Further, network community analysis in TPP bound form shows that state 4 has strong intercommunication between communities formed by P3/ L5 docking region and J3/2 and P2 stem with residual pair G40-A43. The disturbance in residual contact in junction of P4 and P5 stems in states 3 and 5 leads to decrease in information flow between P1, J2/4 to P3/L5 docking region and forms intermediate states in TPP binding. Shortest path in residual network analysis also suggests that state 4 has strong allosteric information flow between A84 (J2/4) to A69 (P3/L5) through the path of most of the known active site residues G17, G18, G19, G21, A43 and A70 compared to states 3 and 5. The formation of tertiary interaction in P3-L5 docking region due to residual contact C21-G72 and C24-A69 and dynamic nature in terminal end of P1 stem leads to conformational switching in different metastable states in TPP binding. Binding energetic analysis using MMPBSA approach suggests that in TPP bound form, intermediate state (state 3) is more energetically favorable. Although binding of TPP is consistent in all the observed five states in Markov model, state 3 has slight higher occupancy in hydrogen bond (G78) and non bonded contacts (A56,C57, C74,G76, C77 and A80). The buried surface area in ligand interacting regions such as U39, G42, A43 of J3/2 G51, A56, C57, C73 of J2/4 and A80 of the junction of P4 and P5 was also higher in states 3 and 5. It was known that degree of metasbolite-induced aptamer switching in response to gene regulation of riboswitch is likely variably tuned. In TPP bound form, intermediate state 3 showed more favourable but information flows between junction J2/4 and P3/L5 docking region is lesser due to extended P1 stem. State 4 has almost similar binding affinity and binding interactions in binding pocket and higher information flows between junction J2/4 and P3/L5 docking region and stable P1 stem was observed. It is pointed out that intrinsic instability in P1 stem in TPP bound form contribute allosteric effect on the dynamics of junction P4 and P5 stem and P3/L5 docking region. Studies also suggest that these events play major role in riboswitch mediated regulatory response during translation initiation. Identifying such crucial residues and also mapping segments of domains where information flow is observed to stabilize the riboswitch (reported in the discussion) are vital for identification of inhibitors to target the regulation of TPP riboswitch. Thus, Markov state 30 ACS Paragon Plus Environment

Page 30 of 37

Page 31 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


based analysis points out the mechanistic understanding in the recognition of TPP in aptamer domain and ligand-induced conformational changes in aptamer domain for conformational switching which are important for gene regulation mediated by kinetic riboswitch. Supporting Information: The Supporting Information is available free of charge on the ACS Publications website at DOI: Markov state model validations and metastable states analysis based on residual fluctuations and intercommunity communication in the residual community network. Author contributions: MK and NHV have performed the molecular dynamics simulation and trajectory analysis of TPP riboswitch. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. DV has supervised the analysis, framing the research and suggestions in writing the manuscript. Funding source: MK is a recipient of Junior research fellowship in DST-PURSE from Department of Science and Technology, Government of India ACKNOWLEDGMENT: Authors thanks for utilising GPU server facility available in CAS in

Crystallography and Biophysics, University of Madras, Chennai, India. References: (1)

Breaker, R. R. Riboswitches and the RNA World. Cold Spring Harb. Perspect. Biol. 2012, 4 .

(2)

Breaker, R. R. Prospects for Riboswitch Discovery and Analysis. Mol. Cell 2011, 43 , 867–879.

(3)

Garst, A. D.; Edwards, A. L.; Batey, R. T. Riboswitches: Structures and Mechanisms. Cold Spring Harb. Perspect. Biol. 2011, 3 .

(4)

Montange, R. K.; Batey, R. T. Riboswitches: Emerging Themes in RNA Structure and Function. Annu. Rev. Biophys. 2008, 37, 117–133.

(5)

Helmling, C.; Klötzner, D.-P.; Sochor, F.; Mooney, R. A.; Wacker, A.; Landick, R.; Fürtig, B.; Heckel, A.; Schwalbe, H. Life Times of Metastable States Guide Regulatory Signaling in Transcriptional Riboswitches. Nat. Commun. 2018, 9 , 944953.

(6)

Steinert, H.; Sochor, F.; Wacker, A.; Buck, J.; Helmling, C.; Hiller, F.; Keyhani, S.; Noeske, J.; Grimm, S.; Rudolph, M. M.; Keller, H.; Mooney, R. A.; Landick, R.; Suess, B.; Fürtig, B.; Wöhnert, J.; Schwalbe, H. Pausing Guides RNA Folding to Populate Transiently Stable RNA Structures for Riboswitch-Based Transcription Regulation. Elife 2017, 6, 21297-21315. 31 ACS Paragon Plus Environment


(7)

Guedich, S.; Puffer-Enders, B.; Baltzinger, M.; Hoffmann, G.; Da Veiga, C.; Jossinet, F.; Thore, S.; Bec, G.; Ennifar, E.; Burnouf, D.; Dumas, P. Quantitative and Predictive Model of Kinetic Regulation by E. Coli TPP Riboswitches. RNA Biol. 2016, 13 , 373– 390.

(8)

Ontiveros-Palacios, N.; Smith, A. M.; Grundy, F. J.; Soberon, M.; Henkin, T. M.; Miranda-Ríos, J. Molecular Basis of Gene Regulation by the THI-Box Riboswitch. Mol. Microbiol. 2008, 67 , 793–803.

(9)

Winkler, W.; Nahvi, A.; Breaker, R. R. Thiamine Derivatives Bind Messenger RNAs Directly to Regulate Bacterial Gene Expression. Nature 2002, 419 , 952–956.

(10)

Serganov, A.; Polonskaia, A.; Phan, A. T.; Breaker, R. R.; Patel, D. J. Structural Basis for Gene Regulation by a Thiamine Pyrophosphate-Sensing Riboswitch. Nature 2006, 441 , 1167–1171.

(11)

Henriksen, N. M.; Roe, D. R.; Cheatham, T. E. Reliable Oligonucleotide Conformational Ensemble Generation in Explicit Solvent for Force Field Assessment Using Reservoir Replica Exchange Molecular Dynamics Simulations. J. Phys. Chem. B 2013, 117 , 4014–4027.

(12)

Sgrignani, J.; Magistrato, A. The Structural Role of Mg2+ Ions in a Class I RNA Polymerase Ribozyme: A Molecular Simulation Study. J. Phys. Chem. B 2012, 116 , 2259–2268.

(13)

Di Palma, F.; Colizzi, F.; Bussi, G. Ligand-Induced Stabilization of the Aptamer Terminal Helix in the Add Adenine Riboswitch. RNA 2013, 19 , 1517–1524.

(14)

Sharma, M.; Bulusu, G.; Mitra, A. MD Simulations of Ligand-Bound and Ligand-Free Aptamer: Molecular Level Insights into the Binding and Switching Mechanism of the Add A-Riboswitch. RNA 2009, 15, 1673–1692.

(15)

Lin, J.-C.; Hyeon, C.; Thirumalai, D. Sequence-Dependent Folding Landscapes of Adenine Riboswitch Aptamers. Phys. Chem. Chem. Phys. 2014, 16 , 6376–6382.

(16)

Banáš, P.; Sklenovský, P.; Wedekind, J. E.; Šponer, J.; Otyepka, M. Molecular Mechanism of PreQ1 Riboswitch Action: A Molecular Dynamics Study. J. Phys. Chem. B 2012, 116 , 12721–12734.

(17)

Wang, W.; Jiang, C.; Zhang, J.; Ye, W.; Luo, R.; Chen, H.-F. Dynamics Correlation Network for Allosteric Switching of PreQ1 Riboswitch. Sci. Rep. 2016, 6, 3100531015.

(18)

Zhang, J.-M.; Jiang, C.; Ye, W.; Luo, R.; Chen, H.-F. Allosteric Pathways in Tetrahydrofolate Sensing Riboswitch with Dynamics Correlation Network. Mol. Biosyst. 2016, 13 , 156–164.

(19)

Chodera, J. D.; Noé, F. Markov State Models of Biomolecular Conformational Dynamics. Curr. Opin. Struct. Biol. 2014, 25, 135–144.

(20)

Trendelkamp-Schroer, B.; Wu, H.; Paul, F.; Noé, F. Estimation and Uncertainty of Reversible Markov Models. J. Chem. Phys. 2015, 143 , 174101-174121.

(21)

Sarich, M.; Banisch, R.; Hartmann, C.; Schütte, C. Markov State Models for Rare 32 ACS Paragon Plus Environment

Page 32 of 37

Page 33 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Events in Molecular Dynamics. Entropy 2013, 16 , 258–286. (22)

Berezovska, G.; Prada-Gracia, D.; Rao, F. Consensus for the Fip35 Folding Mechanism? J. Chem. Phys. 2013, 139 , 35102-35109.

(23)

Haller, A.; Altman, R. B.; Soulière, M. F.; Blanchard, S. C.; Micura, R. Folding and Ligand Recognition of the TPP Riboswitch Aptamer at Single-Molecule Resolution. Proc. Natl. Acad. Sci. U. S. A. 2013, 110 , 4188–4193.

(24)

Soulière, M. F.; Altman, R. B.; Schwarz, V.; Haller, A.; Blanchard, S. C.; Micura, R. Tuning a Riboswitch Response through Structural Extension of a Pseudoknot. Proc. Natl. Acad. Sci. U. S. A. 2013, 110 , 3256-64.

(25)

Bottaro, S.; Gil-Ley, A.; Bussi, G. RNA Folding Pathways in Stop Motion. Nucleic Acids Res. 2016, 44 (, 5883–5891.

(26)

Huang, X.; Yao, Y.; Bowman, G. R.; Sun, J.; Guibas, L. J.; Carlsson, G.; Pande, V. S. Constructing Multi-Resolution Markov State Models (MSMs) to Elucidate RNA Hairpin Folding Mechanisms. Pac. Symp. Biocomput. 2010, 228–239.

(27)

Duesterberg, V. K.; Fischer-Hwang, I. T.; Perez, C. F.; Hogan, D. W.; Block, S. M. Observation of Long-Range Tertiary Interactions during Ligand Binding by the TPP Riboswitch Aptamer. Elife 2015, 4,12362-12379.

(28)

Zgarbová, M.; Otyepka, M.; Sponer, J.; Mládek, A.; Banáš, P.; Cheatham, T. E.; Jurečka, P. Refinement of the Cornell et Al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput. 2011, 7 , 2886–2902.

(29)

Banáš, P.; Hollas, D.; Zgarbová, M.; Jurečka, P.; Orozco, M.; Cheatham, T. E.; Šponer, J.; Otyepka, M. Performance of Molecular Mechanics Force Fields for RNA Simulations: Stability of UUCG and GNRA Hairpins. J. Chem. Theory Comput. 2010, 6 , 3836–3849.

(30)

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; R, J. A. GAUSSIAN03. In: B.05 R, Editor. Pittsburgh, PA: Pople, Gaussian,Inc.; 2003.

(31)

Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. a; Case, D. a. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25 , 1157–1174.

(32)

Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26 , 1668–1688.

(33)

Allnér, O.; Nilsson, L.; Villa, A. Loop-Loop Interaction in an Adenine-Sensing Riboswitch: A Molecular Dynamics Study. RNA 2013, 19 , 916–926.

(34)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 , 926-935. 33 ACS Paragon Plus Environment


(35)

Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 , 3684-3690.

(36)

Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera - A Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25 , 1605–1612.

(37)

Scherer, M. K.; Trendelkamp-Schroer, B.; Paul, F.; Pérez-Hernández, G.; Hoffmann, M.; Plattner, N.; Wehmeyer, C.; Prinz, J.-H.; Noé, F. PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models. J. Chem. Theory Comput. 2015, 11, 5525–5542.

(38)

Bottaro, S.; Di Palma, F.; Bussi, G. The Role of Nucleobase Interactions in RNA Structure and Dynamics. Nucleic Acids Res. 2014, 42 , 13306–13314.

(39)

Seeley, J. R. The Net of Reciprocal Influence; a Problem in Treating Sociometric Data. Can. J. Psychol. Rev. Can. Psychol. 1949, 3 , 234–240.

(40)

Girvan, M.; Newman, M. E. J. Community Structure in Social and Biological Networks. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 , 7821–7826.

(41)

Sethi, A.; Eargle, J.; Black, A. A.; Luthey-Schulten, Z. Dynamical Networks in TRNA:Protein Complexes. Proc. Natl. Acad. Sci. U. S. A. 2009, 106 , 6620–6625.

(42)

Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14 , 33–38.

(43)

Miller, B. R.; McGee, T. D.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E. MMPBSA.Py : An Efficient Program for End-State Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 3314–3321.

(44)

Shivakumar, D.; Williams, J.; Wu, Y.; Damm, W.; Shelley, J.; Sherman, W. Prediction of Absolute Solvation Free Energies Using Molecular Dynamics Free Energy Perturbation and the Opls Force Field. J. Chem. Theory Comput. 2010, 6 , 1509–1519.

(45)

Kulshina, N.; Edwards, T. E.; Ferre-D’Amare, A. R. Thermodynamic Analysis of Ligand Binding and Ligand Binding-Induced Tertiary Structure Formation by the Thiamine Pyrophosphate Riboswitch. RNA 2010, 16 , 186–196.

(46)

Warner, K. D.; Homan, P.; Weeks, K. M.; Smith, A. G.; Abell, C.; Ferré-D’Amaré, A. R. Validating Fragment-Based Drug Discovery for Biological RNAs: Lead Fragments Bind and Remodel the TPP Riboswitch Specifically. Chem. Biol. 2014, 21 , 591–595.


Page 34 of 37

Page 35 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60





Page 36 of 37

Page 37 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


59x44mm (300 x 300 DPI)

ACS Paragon Plus Environment

Conformational Dynamics of thiM Riboswitch To Understand the Gene

Recommend Documents