Differential processing of quorum sensing signals through phosphotransfer: structural insights from molecular dynamics simulations

The auto-inducer-mediated virulence gene expression and biofilm formation in Vibrio sp. uses a highly evolved two-component phosphotransfer system, involving a histidine sensor kinase (LuxQ), an Hpt domain protein (LuxU) and a universal response regulator (LuxO), to process the signal. At low and high cell density, the phosphotransfer reaction occurs differently leading to the activation or deactivation, respectively, of the global repressor which in turn regulates the virulence. Here the molecular details of signal processing and signal decay have been studied using structural modelling and molecular dynamics simulation of LuxQ, LuxU and LuxO individual proteins and protein–protein complexes with and without the phosphate group. The stability, conformational flexibility and structural changes associated with phosphotransfer of the individual protein and the protein–protein complexes are compared. The root mean square deviations and the root mean square fluctuations of the phosphorylated and unphosphorylated proteins showed significant differences in these two processes. The principal component analysis points out the remarkable differences in the essential motions of the systems, which depend not only on the phosphorylated complex but also on the key phosphorylation of the individual protein component. This observation is also highlighted by the dynamic cross-correlation matrix (DCCM) analysis where concerted motions are found to differ depending on the state of phosphorylation. Evaluation of the equilibrated structures and their free energy reveals that the reverse transfer of phosphate during signal decay is energetically less favourable.


Introduction
Quorum sensing is a cell-cell communication process used by bacteria to coordinate gene expression in response to population density (Bassler and Losick 2006;Henke and Bassler 2004a, b;Miller and Bassler 2001), which includes production and detection of chemical inducers, known as auto-inducers (Ng and Bassler 2009;Wang and Ma 2014). The auto-inducer-mediated signal recognition and transduction during bacterial quorum sensing is a well-orchestrated molecular process that occurs through a series of phosphotransfer reactions (Capra and Laub, 2012;Salazar and Laub 1 3 2015). Quorum-sensing signal processing and transduction in Vibrio harveyi, a well-studied process, occur through a modified two-component phosphotransfer reaction involving hybrid sensor/kinase (HK), Hpt and a response regulator (RR) protein (Freeman and Bassler 1999;Stewart 2010;Stock et al. 2000). However, a comprehensive structural and conformational analysis of the detailed phosphotransfer mechanism is not completely understood. The integration and interpretation of quorum-sensing signal and their "stop and go" traffic, at the molecular level, directly depend upon the phosphorylation and the associated conformational transition of the respective proteins involved in the pathway (Swem et al. 2008;Tu et al. 2008).
Two-component signalling system holds a crucial role within the microbial community, which helps them to process critical stimulus-response from the environment and allows them to act as a multicellular organism (Bhate et al. 2015;Papenfort and Bassler 2016;Rutherford and Bassler 2012) to survive various environmental stress (Aguilar et al. 2001;Albanesi et al. 2009;Cheung and Hendrickson 2009;Edurado and Groisman 2001;Laub and Goulian 2007;Ulrich and Zulin 2010). The diversity of environmental signals received by these systems is enormous and hence different molecular mechanisms have evolved to respond to the wide range of environmental stimuli. The basic architecture and specificity are also diverse among different kingdoms of life. Thus, the molecular level understanding of the structure and dynamics of two-component signalling systems (TCS) is extremely important to decode the microbial communication (Faloon et al. 2010(Faloon et al. -2013a. Quorum sensing in V. harveyi activates bioluminescence (lux CDABE light genes) at high cell density (Scheme 1) in which two parallel signalling systems are employed to regulate light production. System 1 includes LuxM, an AHL (acyl homoserine lactone) signal synthase, and LuxN, a hybrid sensor/kinase that respond to the AHL signal (Feng et al. 2015;Henke and Bassler 2004a, b). System 2 includes LuxQ that responds to auto-inducer AI-2. AI-2 binds to the periplasmic LuxP protein, and this complex interacts with the hybrid sensor/kinase, LuxQ, and the signal is passed on to LuxU and LuxO (Cámara et al. 2002;Papenfort et al. 2015). At low cell density, LuxN and LuxQ act as kinases, wherein upon auto-phosphorylation they transfer the phosphate present in its active site Asp residue to the His residue of LuxU and subsequently to the Asp residue of the response regulator LuxO, which activates the production of an unidentified repressor inhibiting bioluminescence (Scheme 1) Tu and Bassler 2007;Waters and Bassler 2006). At high cell density, AHL and AI-2 signal molecules interact with their cognate sensors, LuxN and LuxQ, which now act as phosphatases, leading to the loss of phosphate from LuxU and LuxO. LuxO is now inactivated and the unidentified repressor is not expressed; therefore, the repressor is free to activate expression of the bioluminescence genes (Feng et al. 2015;Rutherford et al. 2011;Shao et al. 2013;Tu and Bassler 2007). AI-2 binding to the LuxP, the membrane-bound receptor, brings definite conformational changes in LuxPQ protein complex which promotes the phosphatase activity (Ulrich et al. 2005) and stops the flow of signal from LuxQ to LuxU. Our interest in this system lies in its similarity to the quorum sensing in Vibrio cholerae. Quorum sensing in V. cholerae represses virulence gene expression and biofilm formation at high cell density using three parallel signalling systems; system 2, LuxSPQ in V. cholerae is identical to system 2 in V. harveyi (Cámara et al. 2002;Henke and Bassler 2004a, b;Papenfort et al. 2015).
We have investigated the structural changes accompanying phosphotransfer in the proteins involved in quorum sensing in Vibrio harveyi during forward and reverse signal processing. In this report, we describe the possible mechanism of signal processing and transduction through phosphotransfer and the associated structural changes in the proteins using protein-protein interaction and molecular dynamics (MD) simulations. This study mainly focuses on the auto-inducer 2 (AI-2) pathway of Vibrio harveyi where Scheme 1 Schematic representation of signal processing during quorum sensing via phosphotransfer reaction in Vibrio sp. The phosphotransfer process has been shown. a Low cell density: the signal transfer through phosphorylation occurs and the global repressor is active and the virulence genes are expressed and b high cell density: the reverse phosphotransfer deactivates the global repressor and the expression of virulence genes are stopped the protein complex LuxPQ acts as signal receptor, LuxU as a histidine kinase and LuxO as a response regulator. For the first time, this report generates the complexes formed by protein-protein interactions, which are then subjected to molecular dynamics simulations to understand the underlying mechanism of phosphotransfer and the associated structural changes required for signal processing in both the directions (Scheme 1). Analysis of phosphorylated and unphosphorylated individual proteins and protein-protein complexes of AI-2 pathway showed the differential changes of the major motions in the protein-protein complexes, revealing the variation in processing and decay of quorum sensing signals.

Methods
The structures used for simulation from the PDB were 1y6d, an NMR structure of LuxU from Vibrio harveyi (Ulrich et al. 2005) with and without phosphorylated His58 (shuttles phosphate from LuxQ to LuxO); 3cfy, the signal receiver domain of LuxO from Vibrio parahaemolyticus (Boyaci et al. 2016) (which has 84% sequence identity to the protein from Vibrio harveyi) with and without phosphorylation at Asp66 which is important in phosphotransfer; 2hje, the histidine sensor kinase subunit LuxQ from Vibrio harveyi (Ulrich et al. 2005;Patskovsky et al. 2008) with and without phosphate at Asp107 (involved in phosphotransfer). Only the first model from the NMR structure of LuxU was considered. For each system, the active site residue (Asp or His) was phosphorylated. For phosphohistidine, parameters corresponding to single protonated phosphate group, phosphorylated at NE1, HIS(NE1)-PO2(OH) were taken from http://resea rch.bmh.manch ester .ac.uk/bryce /amber /. The parameters for phosphoaspartate, on the other hand, were separately generated by employing ANTECHAMBER (Wang et al. 2004) on the basis of a representative phosphoaspartate structure downloaded from PDB (http://ligan d-expo.rcsb.org/repor ts/P/PHD/index .html). The complexes LuxUO and LuxUQ have been modelled using docking software GRAMM-X (Tovchigrechko and Vakser 2006). From the literature (derived from the papers citing the monomer structures) (Ulrich et al. 2005;Boyaci et al. 2016;Neiditch et al. 2005), we had selected a list of residues important for binding and these were provided as input to the GRAMM-X server. After docking, the best model ranked by the GRAMM-X procedure, based on its scoring function, was taken for further analysis and simulations.
Simulations were performed for the phosphorylated as well as the unphosphorylated structures in explicit solvent at 300 K for 40 ns employing the Sander module of AMBER 10.0 software package (Case et al. 2005). Periodic boundary conditions and the particle-mesh Ewald method for the treatment of long-range electrostatic interactions were used for the simulation. The protein was solvated using TIP3P box with the FF99SB force field parameters at a 9 Å radius surrounding the protein. The system (protein and the solvent) was neutralized using Na + or Cl − ions depending on the net charge of the system. The molecules were minimized using steepest descent for 500 cycles, followed by conjugate gradient method for 20,000 cycles prior to dynamics run. All systems were heated to 300 K and equilibrated within 40 ps following minimization. Bonds involving hydrogen were constrained with the help of the SHAKE algorithm. The production run at 300 K was carried for 40 ns, employing constant pressure periodic boundary conditions. A non-bonded cut-off distance was set to 12 Å for columbic and electrostatic interactions, and 2 fs integration time step was used. Replicates were run for each system thrice (each comprising of 40 ns simulation time); then RMSD and Rg (Lobanov et al. 2008;Torre et al. 2000) were compared to check for consistencies among the production runs.
The changes in the trajectories were analysed employing quantitative measures, such as RMSD (measured in Å which tells us about the gross structural changes taking place in the structure), root mean square fluctuations (RMSF, the average flexibility of each residue in the protein) and radius of gyration (Rg, the compactness of the protein molecule). Both RMSD and RMSF were calculated using backbone atoms (C, O, N and C α atoms) in the ptraj module of amber. Rg was calculated using only the C α atoms using an in-house PERL program (Lobanov et al. 2008;Torre et al. 2000). The energy of the molecule was estimated by generalized Born solvent accessible surface area, GBSA, employing ffgbsa module of amber (Onufriev et al. 2002). We have also employed the Consurf server (http://consu rf.tau.ac.il/2016/) (Glaser et al. 2003;Landau et al. 2005;Ashkenazy et al. 2016;Celniker et al. 2013;Ashkenazy et al. 2016) to identify regions of functional importance in terms of evolutionary conservation for all three monomers.
Further analysis of the motions was performed using the principal component analysis (PCA) and dynamic cross-correlation maps (DCCM). The identification of residues and regions that move in or out of phase during simulations can be achieved through dynamic cross-correlation maps. The calculation is performed based on the equation given below: The symbol <> indicates an average over the sample period t, t k is the time at that instant ∆x i and ∆x j are vectors that represent the motions of residues i and j. The crosscorrelation coefficient, C ij , varies from − 1 (completely anticorrelated motion) to + 1 (completely correlated motion). In 1 3 this study, the cross-correlation coefficients were computed by considering the C α atoms using the trajectory obtained after 40 ns of simulation employing the ptraj module of Amber (Case et al. 2005), and the plots were generated in Gnuplot (http://www.gnupl ot.info/). Principle component analysis (PCA) is a mathematical technique generally used to find patterns in high-dimensional datasets, such as protein structures or simulation trajectories. It is a technique to reduce multi-dimensional data sets to lower dimensional data. PCA is often used to decipher relationships or patterns in MD simulation trajectories that can help, detect global, correlated motions of the system (the principal components) (Lauria et al. 2009;Yang et al. 2009). PCA is used to extract larger amplitude motions (most prominent changes) observed in MD trajectories through the eigenvectors (principal components, PCs) of the mass-weighted covariance matrix (C) of the atomic positional fluctuations. The dynamics of the subsets of PCS are termed as essential dynamics (ED) and are instrumental while studying large conformational changes. The eigenvectors ( ⃗ v ) and eigenvalues (λ) are calculated from the fluctuation covariance matrix that gives us the principal components of the system: where x i is the coordinate of the protein and the equivalent coordinate for the reference structure is denoted as x 0 i , and the average ensemble is denoted by the braces ˂ ˃. PCA was performed using only C α of the structures from the trajectory of 40 ns employing the Bio3D package of R (Grant et al. 2006). The movies were generated using VMD (Humphrey et al. 1996) and figures using PyMOL (http://www.pymol .org) (The PyMOL Molecular Graphics System, Version 1.8 2019). The calculations have been performed considering only the structural/dynamic effects of the presence or absence of a phosphate group on our domains of interest. Presence of ATP or other co-factors has not been considered for the present study.

Results and discussion
We present the results on the three individual proteins, first unphosphorylated and then phosphorylated. Then we discuss the two complexes, LuxUQ and LuxUO, in both cases considering the unphosphorylated state first, followed by phosphorylation occurring at one or the other active site. Table S1 provides the summary of RMSD values between the first and the last structures from simulations, and the average energy of all structures along with the standard deviation as well as some other details from the analysis performed on the trajectories. The RMSD and Rg show qualitatively similar values for each replicate MD simulation (data not shown).
The protein is composed of four helices and six anti-parallel β-strands (Fig. 1). The structure does not undergo much change except at about 6.5 ns and at about 35 ns (RMSD being 3.9 Å); however, the RMSD stabilizes to 2.3 Å towards the end of 40 ns (Fig. 2a). A stretch of residues, 127-137, comprising parts of two β-strands and the turn connecting them (encircled in red, Fig. 1a) was found to move towards the core of the molecule. A zoomed view of Asp107 and the changes in its conformation during the time period is presented in figure S1a.

On phosphorylation
Not much change was observed on the phosphorylated form of LuxQ (Fig. 1b); RMSD between first and last structures of the simulation is 2.5 Å (Fig. 1c). Similar movement in a loop consisting of active site residues Ser106 and Asp107 was observed in both the phosphorylated and the non-phosphorylated structures ( Figure S1, only Asp107 is shown). RMSD values are quite similar till ~ 30 ns for both the phosphorylated and unphosphorylated structures; the secondary structures remain the same during the simulation. The RMSF profiles for both the structures show a similar trend. Subtle changes were observed in hydrogen bonding patterns on phosphorylation, the hydrogen bonds between Asp107 and Ser109 are maintained in both structures, although the percentage is decreased when Asp107 is phosphorylated, which also forms new hydrogen bonds with Lys201 (Table S2a). Rg of a protein is a measure of its compactness, as mentioned before, and indicates if the molecule is properly folded; the structures from simulation should maintain a relatively steady value of Rg (Lobanov et al. 2008;Torre et al. 2000). We find quite similar Rg values for the unphosphorylated and phosphorylated forms of the protein throughout the time period (18.7 ± 0.1 and 18.6 ± 0.1 Å, respectively). The stability of phosphorylated structure has decreased from the one without phosphate (− 5362 ± 91 vs. − 5606 ± 89 kcal/mol, energy estimated by generalized Born solvent accessible surface area, GBSA (Onufriev et al. 2002).

LuxU with unphosphorylated and phosphorylated His58
Unphosphorylated LuxU is a phosphotransferase protein consisting of a fourhelix bundle (His58 is located in the centre of a helix in the bundle, Fig. 2) (Ulrich et al. 2005). During the time of simulation, large amounts of structural deviations were observed for the whole structure (RMSD between first and last structures of simulation is 4.9 Å (Fig. 2a), affecting some portions in the four-helix bundle that were disrupted. Percentage of helicity decreases from 57 to 38% at the end of simulation. The trail of motion of His58 can be observed in figure S2a.

On phosphorylation
The structural changes were spread all over the phosphorylated structure as seen during the 40 ns of simulation (Fig. 2c, RMSD between first and last structures from simulation is 6.4 Å) and a part of a helix in the bundle was disrupted although overall percentage helicity increased from 54 to 64%. On comparing the energy profile of the two structures, the stability of the phosphorylated form was found to be more than the unphosphorylated one (− 3166 ± 50 vs. − 3008 ± 49 kcal/mol). The reason for additional stability could be due to the presence of positively charged residues (including Lys54 and Lys61) in the active site around the His58, consistent with the binding site for the negatively charged phosphate ( Figure S3). The phosphorylated LuxU undergoes more structural deviations than when unphosphorylated and the change takes place gradually after 2 ns of simulation (Fig. 2c). The RMSF (root mean square fluctuations measure the average fluctuations of the residues in the structure and give an insight into the flexible regions of the protein. The protein with phosphorylated His58 has more flexibility than the non-phosphorylated one (Fig. 2d). A subtle increase in Rg is also observed (14.8 ± 0.2 vs. 15.4 ± 0.3 Å between unphosphorylated and phosphorylated forms), indicating higher looseness of the protein on phosphorylation. The hydrogen bonding partners of His58 change on phosphorylation, the bond with Ser62 although present when non-phosphorylated, the duration or the percentage of occurrence of the hydrogen bonds during the time period of simulation is increased when His58 is phosphorylated. Lys54 and Lys61 that form a part of the positive pocket around His58 have increased occurrence of hydrogen bonds when His58 is phosphorylated (Table S2b). The flipping of the imidazole

Unphosphorylated
The LuxO protein receives the phosphate group from LuxU. It is composed of both α-helices and parallel β-strands, the amphipathic five helices surround the five strands and the active site residue Asp66 is present in the middle of a helix (Fig. 3). The protein is very stable and does not undergo much change during simulation (RMSD between first and last structures of simulation is 1.6 Å, Fig. 3c). No changes were observed for secondary structural elements as well.

On phosphorylation
The phosphorylated structure shows more change than the unphosphorylated one (RMSD between the first and last structures is 2.02 Å, are shown in Fig. 3b and RMSD plot is shown in Fig. 3c). A swaying movement of active site residue Asp66 is observed in both phosphorylated and unphosphorylated structures ( Figure S4). The energetics of the two systems are very similar (− 3533 ± 60 vs − 3478 ± 59 kcal/ mol). Slight changes in the hydrogen bonding patterns of Asp66 take place when phosphorylated. The presence of hydrogen bond between Asp66 and Lys95, although insignificant when non-phosphorylated (0.06%), showed substantial presence when Asp was phosphorylated (36% as in Table S2c). Two other hydrogen bonds involving Asp66 (with Glu70 and Glu62) showed decreased presence upon phosphorylation. The phosphorylated structure undergoes greater structural deviations, as evident from its RMSD and RMSF profiles (Fig. 3c, d).
Among the three systems, the stability of LuxQ is more than that of LuxO that is in turn more than of LuxU. Both LuxQ and LuxO undergo much less structural deviations when compared to LuxU, even in the phosphorylated state, as seen from both RMSD and RMSF profiles (Figs. 1, 2, 3 respectively). The effect of phosphorylation is opposite in case of LuxU as mentioned earlier (stabilizes contacts and decrease in total molecular mechanical energy), for LuxO the difference is negligible on phosphorylation and the phosphorylated structure is less stable in case of LuxQ (Table S1).
Conservation analysis of positions among member proteins from the same family can often reveal the importance of each position for the protein's structure or function. Therefore, the three proteins were also compared in terms of which regions were conserved ( Figure S5) using Consurf server (Glaser et al. 2003;Landau et al. 2005;Ashkenazy et al. 2010;Celniker et al. 2013;Ashkenazy et al. 2016), and we observe that all three proteins have regions of high evolutionary importance. The region containing the active site residues showed more conservation in all three, especially for LuxQ and LuxU. For LuxO, the active site region is slightly less conserved. Surprisingly, LuxU and LuxO have regions of high variability adjacent to the conserved active site region. Thus, not only the residues at the phosphorylation site, but also the region around it are important in all three structures.

LuxUQ, unphosphorylated and with phosphorylation at Asp107 or His58
In the LuxUQ complex, the phosphate group is passed from Asp107 in LuxQ to His58 in LuxU. The structures from the start, middle and end points of simulation are shown in Fig. 4a. Here we found His58 to be quite far apart from Asp107 at start and they remain apart through the time period. This is also reflected in figure S7a where the distance between His58 NE2 and Asp107 OD2 versus time was plotted. The overall structural deviations are within 2-3 Å (Fig. 4c). Compared to LuxQ subunit, the LuxU subunit undergoes more structural changes ( Figure S6), as also observed from the rearrangement of the helices and from the RMSF profile as well (movie SM1, Fig. 4d). These changes do not affect the assembly of the complex. The hydrogen bonds between His58 and Leu53, Lys61 are mainly observed, while Asp107 forms hydrogen bonds mostly with Ser109 and His110 during the time period (Table S3a). The average Rg is 20.7 ± 0.1 Å.

Phosphorylated Asp107
Major changes were observed when Asp107 was phosphorylated, parts from two helices of the LuxU subunit (residues 75-100) moved away from the other two helices near the interface of the complex (Fig. 4b, movie SM2), while the active site residues His58 and Asp107 remain within (~ 7-8 Å towards the end of 40 ns, Figure S7a). This large structural change is clearly reflected in the RMSD and RMSF plots (Fig. 4, movie SM2), where the complex with phosphorylated Asp107 shows larger values compared to both unphosphorylated and phosphorylated His58 structures. This is also evident from the Rg value that increases significantly from the unphosphorylated complex (23.1 ± 0.9 vs 20.7 ± 0.1 Å). Hydrogen bonding patterns changes upon phosphorylation at Asp107, hydrogen bonds with Ser109 (LuxQ), Ser57, Lys61 (near His58, LuxU) predominate during 40 ns (Table S3a). Fig. 4 The protein-protein docked structures of LuxUQ are shown here from the beginning, middle and end points of simulations; LuxU is shown in green and LuxQ in slate blue, and the residues important for phosphorylation are shown as sticks, His (LuxU) in red and Asp (LuxQ) in blue. The non-phosphorylated structure is seen in a, while the complexes with phosphorylated Asp and His are shown in b, c, respectively. In b the four helices in LuxU are numbered sequentially; helix H2 contains the active site His. The backbone d RMSD and e RMSF plots of the LuxUQ complexes are compared. Residues 1-114 belong to LuxU (chain A), and residue numbers 115-324 belong to LuxQ (chain C), the residue numbering being sequential

Phosphorylated His58
His58 is initially at a distance of ~ 10 Å from active site residue Asp107, and remains so for a while and then move further apart towards the end of 20 ns ( Figure S7a). This is brought about by the movement of the LuxU domain which was found to shift away from the LuxQ domain (Fig. 4c, movie SM3). Quantitatively, this change can be observed from the RMSD and RMSF profile (Fig. 4). The RMSD gradually increases to ~ 6 Å in 15 ns, but gradually the structure was stabilised as the LuxU moves away (3-4 Å near 40 ns) from LuxQ. The RMSF plot shows the LuxU domain to be more flexible compared to the LuxQ domain, the former showing more structural changes ( Figure S6). Compared to a Rg value of 20.7 ± 0.1 in the unphosphorylated complex phosphorylation at His58 increases the value to 22.1 ± 0.3 Å. Phosphorylated His58 mainly forms hydrogen bonds with Met34, Ser77 (in LuxQ) and Glu33 (Table S3a).
The total molecular mechanical energy of the three complexes was estimated and it was seen that phosphorylation at His58 increases the stability of the complex (− 8748 ± 131 kcal/mol) when compared to the unphosphorylated structure (− 8632 ± 53 kcal/mol) and that with phosphorylated Asp107 (− 8548 ± 131 kcal/mol).

Principal component analysis and collective modes of motion in the unphosphorylated and phosphorylated LuxUQ
To further our understanding of the major motions seen in the trajectory, we have performed PCA (Lauria et al. 2009;Yang et al. 2009) on the MD trajectories within the time interval 1-40 ns. Thereafter, the principal components (PCs) were sorted according to their contribution to the total fluctuation along the ensemble of conformation. From the proportion of variance plot or screen plot ( Figure S9), one can observe that the first two PCs account 33% and 17% of the variance, respectively, and that the first three PCs cumulatively explain ~ 57.7% of the variance in LuxUQ without phosphorylation. Similarly, for the phosphorylated structures most of the variance was captured within the first three components (e.g. 74, 81 and 83% for LuxUQ-phosphorylated Asp107; and 29, 44 and 54.6%, for phosphorylated His58, respectively, Figure S9). This indicates that the first three principal motions can capture most of the internal motions of our system. The direction of motions of the first three PCs for LuxUQ can be seen in the porcupine plots (Figs. 5, S10).
For the complex consisting of phosphorylated Asp107, the first eigenvector alone accounts for 74% of the motion observed in the 40 ns of the simulation trajectory, whereas the contribution of the first PC for the unphosphorylated structure and the one with phosphorylated His58 was much less (~ 30%). This motion in the phosphorylated Asp107 complex corresponds to an opposing movement of the residues 75-100 in LuxU with respect to the other two helices in the four-helix bundles and the LuxQ subunit (movie SM5), a sort of hinge-like movement that leads to the loosening up of the assembly. The second and third modes correspond to mainly local twisting and turning involving mostly the LuxU subunit. This is in contrast to the motion of the first eigenvector obtained from trajectories of LuxUQ unphosphorylated and with phosphorylated His58 (Fig. 5a, c, movies SM4 and SM6), where the essential modes are influenced by all three PCs. For the unphosphorylated complex, the eigenvector corresponds to the two subunits moving closer towards each other whereas in the phosphorylated His58 structure, the LuxU subunit was seen to move away from LuxQ. PCA confirms that all three differ significantly in their essential motions.

Dynamic cross-correlation analysis for the unphosphorylated and phosphorylated LuxUQ
A two-dimensional dynamical cross-correlation map (DCCM) was calculated (using only C α atoms) (Grant et al. 2006) to reveal the concerted motions of the unphosphorylated and phosphorylated trajectories of LuxUQ. The colour ranges from blue (completely anti-correlated) to red (correlated or moving in the same direction). A diagonal point represents the C α atom of the same residue along both the axes; therefore, the diagonal elements show the maximum correlation. Figure 6 shows the three DCCMs to be quite different. In the unphosphorylated structure (Fig. 6a), most of the residues are either not or weakly correlated, except for the dynamics of the residues 55-60 and 105-110 in LuxU subunit. This is consistent with the fact that from PCA we could see mostly local motions as opposed to concerted global changes (movie SM4). On phosphorylation, the correlated motions have changed significantly, in case of phosphorylation at His58, the intensity of anti-correlation increases especially for residues 115-200 and 220-250 in the LuxQ subunit, while the residues 55-60 and 105-110 of the LuxU subunit show increased intensity of correlated motion (Fig. 6c). However, in the phosphorylated Asp107 complex, the correlated motions have changed drastically, the residues 75-100 move in an anti-correlated fashion with respect to residues 0-50, 101-170 and 220-320 which are in turn strongly correlated to each other (Fig. 6b). We had observed a similar motion from the first eigenvector of PCA, where the residues 75-100 in LuxU move in the opposite direction to the rest of the protein. This indicates that phosphorylation at Asp107 may have a more pronounced effect on the dynamics of the complex compared to when His58 was phosphorylated and the unphosphorylated complex as well.

LuxUO, unphosphorylated and with phosphorylation at His58 or Asp66
In the LuxUO complex, the phosphate group is passed from His58 in LuxU to Asp66 of LuxO. At high cell density, the phosphate is stripped from LuxO thus downregulating the bioluminescence genes. When no phosphate group is present in the LuxUO complex, simulation reveals how His58 and Asp66 come closer and then move apart within the time period (Fig. 7a, S7b and movie SM7). The structural deviations are observed for the whole complex (Fig. 7c, RMSD between the first and the last structures from simulation is 6.6 Å), although the domain corresponding to LuxU shows more change (RMSD increases to 4.6 Å, Figure S8) compared to the LuxO (RMSD, and 2.4 Å) part in the complex. This is also depicted in the RMSF profile where we see the residues of LuxU having higher flexibilities (residues 1-114, Fig. 7d) compared to the residues (115-244) in the LuxO domain. A part of the first helix in the helix bundle of LuxU (residues 40-45) moves towards the helical end of LuxO (residues 126-130), but they move away from each other at the end of 40 ns (Fig. 7a, movie SM7). Overall, the structure remains compact throughout the simulation, Rg is 19.4 ± 0.3 Å. His58 was observed to form hydrogen bonds with Ser62, Asp66 (in LuxO) for the most part of 40 ns, and Asp66 (LuxO) chiefly forms contacts with Gln70 (LuxO) ( Table S3b).

Phosphorylated His58
As seen in the unphosphorylated complex, when His58 was phosphorylated, the active site residues were close to each other (remain within ~ 4.5 Å) initially but after 12 ns they move further apart (Figs. 7b,S7b,movie SM8). A lot of rearrangements of helices are observed especially in the LuxU domain, this subunit undergoes greater structural deviations (RMSD values after 40 ns are 6.7 and 3.3 Å, for LuxU and LuxO, respectively, Figure S8). Rg values were similar to the unphosphorylated structure, 19.6 ± 0.3 vs 19.4 ± 0.3 Å, respectively. His58 primarily forms hydrogen bonds with Ser62, Lys54, Lys61 and Lys95 (LuxO) ( Table S3b).

Phosphorylated Asp66
The active site residues (His58 and Asp66) were quite distant from each other at the start of simulation (~ 7 Å) and move further apart during the time period (~ 14 Å,Figs. 7c,S7b). The two subunits were also seen to move away from each other, affecting the assembly of LuxUO complex (movie SM9). The overall structural changes are comparable in all three LuxUO complexes (Fig. 7d), LuxU was observed to be more flexible (residues 1-144, Fig. 7c) and its RMSD is greater in all three complexes ( Figure S8), especially for the complex with phosphorylation at His58. Rg values (19.6 ± 0.3 Å) are similar to the values in the other two LuxUO complexes, discussed above. The hydrogen bonding partners of Asp66 include Gln70, Lys95 (both in LuxO) and Lys61, Arg71 (of LuxU subunit) ( Table S3b).

Principal component analysis and collective modes of motion in the unphosphorylated and phosphorylated LuxUO
As described for LuxUQ, we have also performed PCA on the simulation trajectories of LuxUO considering both unphosphorylated and phosphorylated systems. In the screen plot ( Figure S11a), one can see that the first principal Fig. 7 The structures of LuxUO are taken from the beginning, middle and end points of simulation; LuxU is shown in green cartoon and LuxO in orange, the residues important for phosphorylation, His (LuxU, red) and Asp (LuxO, blue), are shown as sticks. The non-phosphorylated structure is seen in a, while the complexes with phosphorylated His and Asp residues are shown in b, c, respectively. The d RMSD and e RMSF plots of the LuxUO complexes. Residues 1-114 belong to LuxU (chain A), and residue numbers 115-244 belong to LuxO (chain B) ◂ component (PC1) explains almost two-third of the overall variance (70%), the first three components explain more than 80% of the complete variance. The global motion as seen in PC1 describes a gradual twisting movement about the interface bringing the two subunits and the residues His58 and Asp66 closer (Fig. 8a, movie SM10). For the complex with phosphorylation at His58, PC1 accounts for 54% of the variance ( Figure S13b), here also a twisting motion dominates bringing the two domains together although the direction of the motions differ (Fig. 8b, movie SM11). The marked difference in the LuxUO complex containing phosphorylated Asp66 is that for this system PC1 describes ~ 40% of the variance and the motion defined by it moves the two subunits apart rather than bridging them closer ( Figure S11c, Fig. 8c, movie SM12). The motions described by PC2 and PC3 for the three complexes can be seen in Figure S12.

Dynamic cross-correlation analysis for the unphosphorylated and phosphorylated LuxUO
Additional changes in the direction of motion due to phosphorylation were looked into with the help of dynamic cross-correlations of atomic fluctuations averaged over the simulation period of 40 ns. The results, given in colourcoded form, show correlated and anti-correlated regions of the structure (Fig. 9). For the unphosphorylated complex, we found deeper shades of red and blue, distinguishing regions of high correlation and anti-correlation that corroborate with the fact that mostly global motions are seen for this system. The residues 0-50, 60-70, 130-150 and 220-240 have concerted motions and they move in anti-correlated fashion with the rest of the structure (Fig. 9a). Phosphorylation changes Fig. 8 Porcupine plots showing direction of motion seen in PC1 for LuxUO a without phosphorylation, b with phosphorylated His and c phosphorylated Asp. In the three plots, PC1 accounts for 70, 54 and 39%, respectively, of the variation seen in the trajectories the direction and the cooperativity of motion significantly. In the phosphorylated His58 complex, we find residues 100-150 moving in the same direction as 210-230 (Fig. 9b). A similar trend was observed in case of phosphorylation at Asp66 (Fig. 9c), but here the intensity of correlation and anti-correlation in the map was less pronounced compared to the other two complexes.

Conclusion
To summarize, among the three proteins, LuxU, the one at the centre (Scheme 1), was found to be the most flexible, both in phospho-and unphosphorylated forms. LuxQ was less stable on phosphorylation, which may lead the transfer of the phosphate group to LuxU, which has a number of positive groups surrounding the active site His to stabilize the negatively charged phosphate. LuxO does not undergo much change structurally or energetically on phosphorylation.
The free energy values generated from MD simulated structures of LuxU, LuxQ and LuxO (with and without phosphate) complexes are compared in Scheme 2. The scheme reveals how the free energy changes during complex formation and how they differ when different proteins are phosphorylated. In LuxUQ complex, the structure is more stable when the phosphate group is at His58 of LuxU compared to the other two structures. Similarly, in LuxUO complex, the structure is most stable (among the three) when the phosphate group is at His58, explaining clearly that the formed phosphotransferase is energetically more favourable that the reverse one.
We considered the changes that occur in the two complexes. When there is no phosphate group in the structure, LuxU shows more change as compared to LuxQ, as is the case when the two proteins are isolated. Among the two phosphorylated states (at either Asp107 of LuxQ or His58 of LuxU), when the phosphorylation is at LuxQ, this component undergoes major structural (mainly global) changes and becomes less stable, which may cause the transfer of the phosphate group to LuxU. The complex with phosphorylated LuxU that showed large local changes in LuxU at the initial stage eventually gets stabilized. The energy of this complex is the minimum compared to the unphosphorylated forms, or the complex where LuxQ is phosphorylated. This indicates that the system is geared such that LuxU is phosphorylated. Even in the LuxUO complex, LuxU is more flexible. The major motions as seen in the two complexes (LuxUQ and LuxUO) differ depending on the site of phosphorylation (Asp or His). For example, in LuxUQ, when Asp is phosphorylated, there is more of a global movement effecting the position of two of its helices H3 and H4. When His is phosphorylated, the trajectories are dominated by more local changes. Interestingly, when the same His is phosphorylated in the LuxUO complex, the changes are There are more local changes when Asp is phosphorylated. Thus, considering the flow of phosphate from LuxQ to LuxU and then to LuxO through the participation of two separate complexes LuxUQ and LuxUO, the component that receives the phosphate group exhibits more of local motions on being phosphorylated.
In both Vibrio species, the phosphate is channelled via a single phosphotransfer protein LuxU in three different systems. The centrality of LuxU is also seen in the dynamics of LuxU that is significantly different from LuxQ and LuxO. Of the three proteins, LuxU is the only one, which is stabilized on phosphorylation; this is because of the presence of positive residues (Lys60 and Lys61) near the phosphorylation site. Simulation also reveals LuxU to be more flexible even in the complex assemblies.
Phosphorylation changes the major motions of the complex, as seen from principal components and the dynamic cross-correlation maps. Depending on the site of phosphorylation (His or Asp) the major motions of the complex differ significantly, indicating that the steps of phosphate transfer along the pathway are not interchangeable. In conclusion, the study brings about the important role played by the dynamics in protein structures that may drive a series of reaction steps to proceed along a particular direction leading to the observation of a particular phenotype.