After obtaining all the simulation results the most important step to take is to analyze the data, so to give appropriate interpretations of what is happening during the evolution of the system. Choosing the right analysis tool for the kind of information we want to extract is of prime importance. Here we will briefly illustrate the methods used in the work.
RMSD and RMSF
The Root Mean Square Deviation (RMSD) is defined as,
RMSD=
√∑Natoms
i=1 (ri(t1)−ri(t2))2 Natoms
(3.1) where the Cα-backbone is used to compare two conformations of the same protein in terms of structural similarity. In the case of an MD simulation, a reference structure has to be chosen to which all other snapshots are compared. In our case, the initial structure.
The Root mean square fluctuation (RMSF) indicates the flexibility of a molecule in general, is
i.e. the deviation of the atomifrom his mean position at time-step T. Again the Cα-backbone structure is used as reference.
In general the RMSD gives information on deviations from a conformation or position, while the RMSF indicate an overall increased structural flexibility.
Secondary structure assignment
Secondary structure can be determined by the characteristic values of backbone’s Φ and Ψ angles in combination with hydrogen bond energies, as done by the algorithm STRIDE (Structural identifi-cation) [81]. The STRIDE energy function contains a hydrogen-bond term containing a Lennard-Jones-like 8-6 distance-dependent potential and two angular factors reflecting the planarity of the hydrogen bond geometry. It uses, as training data, a verified set of secondary structural elements deposited by crystallographers in the protein data bank and can give accurate predictions in the 60%
of cases.
Results of the secondary structure evolution of the trajectories are presented by residues along the time-line with the color code reflecting the adopted secondary structure.
Dynamic cross correlation
To better describe the dynamic of our systems, we analyzed resulting trajectories by mean of the dy-namic cross correlation (DCC). This method has been extensively applied to quantify the amount of correlation between atom movements. DCC between two atoms (i, j) is defined as follow
DCC(i,j) = √ ⟨Δri(t)·Δrj(t)⟩ anticorrelated motions, to a value of +1, strongly correlated motions. The coordinate vectorri(t) can be taken along the x, y, z or any combination of these coordinates.
Hydrogen bonds and salt bridges
Hydrogen bonds and salt bridges are among the most important stabilizing interactions in protein and knowing how many there are can give important information on the system stability. Hydrogen bonds can be formed between two groups that are donor and acceptor when at the defined distance
and angle, for protein usually between 2.5 and 3.2 Å and 180, giving interaction energies of around 1 Kcal/mol. The distance cutoff used in the analysis is 3 Å and for the angle the cutoff is 20 . Salt bridges are stronger interactions than hydrogen bonds, ranging around 4 Kcal/mol. They occur between amino acid side-chains with opposite positive or negative charges, namely, Glu- or Asp-vs. Arg+ or Lys+, or charged molecules present in solution, such as acids or ions. A salt bridge is generally considered to exist when the centers of charge are 4 Å or less apart. In our analysis we used a cutoff of 4 Å.
Results 4
4.1 DrBph
4.1.1 Dynamic Properties of the Phytochrome’s dimer Photosensory Domain¹ Eukaryote phytochromes are postulated to be obligate dimers[4–6] to ensure complete photocon-version and proper signal transduction[83–85]. Moreover some intermediate conformers can be attributed to dimer’s mixed states in which, for example, one chain is in the Pfr state and the other chain in the Pr state. [14,86–88].
Most computational studies focus on the chromophore and its binding pocket[89–91], but little at-tention was devoted to the influence of a dimer versus a monomer on the overall dynamic properties.
In this section we address this lack simulating DrBph in the Pr state by means of classical atomistic molecular dynamics employing two different force fields, AMBER ff14SB and CHARMM27. We also tested two different approaches, namely, cMD and aMD, to explore a greater conformational space. Special attention is given to the distinct dynamic behaviour and correlation of the DrBphP monomer and dimer forms.
¹The content of this section was published, in a more comprehensive form, in 2020 in the Journal of Phys-ical Chemistry B under the title ”Dynamic Properties of the Photosensory Domain of Deinococcus Radiodurans Bacteriophytochrome”[82]. Figures and tables are taken from there.
Preparation
The initial models for DrBphP monomer and dimer were built by homology modeling with SWISS-MODEL taking as template the crystal structure in the Pr state (2,74 Å, PDB entry 4q0j)[14]. The crystal waters were added to the resulting structure by aligning it to the original structure. Titratable side chains were protonated by estimating their pKa values withKarlsberg 2+program[92] at pH 7.0 and assigning their relative protonation. Finally the protein was solvated in a rectangular box of TIP3P water with Åfrom the edges and the charge was neutralized with sodium ions.
Simulations set up
The simulations with the CHARMM27 force field used the biliverdin (BV) chromophore parame-ters [93] with some corrections on the dihedral parameparame-ters[82]. All simulations were done on GPU using a 2 fs time step under periodic boundary conditions with the particle-mesh-Ewald method with a Å cutoff for the van der Waals interaction and hydrogens were constrained with the SHAKE algorithm. The system was minimized for 40 ps, then heated for 60 ps (NVT) to 300 K using Langevin dynamics and pre-equilibrated for 60 ps (NPT) using the Langevin Piston Nose-Hoover method [94] at 300 K. The 500 ns cMD and aMD (dual boost) production simulations at 300 K (NPT) were performed for both DrBphP monomer and dimer, starting from the pre-equilibrated structure.
For the aMD simulation with NAMD(v2.12), we used the following boosting parameters (units are kcal/mol)Edihed = ⟨Vdihed⟩+2Nres, αdihed = 2
5NresandEtot = ⟨Vtot⟩+0,14Natoms, αtot = 0,14Natoms, for monomer, andEdihed = ⟨Vdihed⟩+1.2Nres, αdihed = 1.2
5 NresandEtot = ⟨Vtot⟩+ 0,10Natoms, αtot =0,10Natoms, for dimer. Mean values of the dihedral and total potential energy were obtained from a 10 ns production run. We used slightly lower coefficients than the suggested ones (for example in ref. [95]) in order to avoid destabilization of the system.
For the simulations with the AMBER ff14SB force field [96], performed with AMBER(v16), with BV parameters[82]. Also these simulations were performed on GPUs with a 2 fs time step under periodic boundary conditions with the particle-mesh-Ewald method for electrostatic interactions, a cutoff of 12 Å for the van der Waals interaction and hydrogens constrained with the SHAKE algo-rithm. The setup was minimized for 40 ps, heated from 0 to 100 K (NVT) in 60 ps and then from 100 to 300 K (NPT) and finally pre-equilibrated at 300 K (NPT) for 120ps. Starting from the last structure a 500 ns cMD and GaMD (dual boost) simulations were performed at 300 K (NVT) for the monomer and dimer. For the GaMD simulation we set the boost potentialE = Vmaxand the upper limit for the standard deviationσ0 =5 kcal/mol for the monomer andσ0 =3 kcal/mol for the dimer. The values forVmax,Vmin,VavgandσV were automatically obtained from an initial 2 ns cMD simulation.
Results
To better evaluate the results, the trajectories were evaluated in terms of RMSD of heavy atoms with respect to the crystal structure and their RMSF values, secondary structure evolution plots and dy-namic cross correlation diagrams.
The RMSD and RMSF values obtained from the cMD and aMD simulations using AMBER ffSB14 and CHARMM27 force fields are comparable. The average RMSD values computed from the cMD and aMD simulations are about 4 Å and 6 Å (Table 4.1.1), respectively, indicating that the sampled dynamics do not strongly depend on the chosen computational approach. The higher RMSD values for the accelerated MD simulations (Figure 4.1.1) are related to the increased con-formational space exploration, as illustrated in the Ramachandran plots (Figure 4.1.2).
Table 4.1.1: Average RMSD values, in Å and standard deviation (in parentheses), of Cαpositions for monomer and dimer structures of DrBphP computed over the last 100 ns MD simulations using CHARMM27 and Amber ffSB14 force fields with respect to the crystal structure.
RMSD avg. last 100 ns Monomer Dimer CHARMM27 (cMD) 3.55 (0.66) 3.79 (0.38) AMBER ffSB14 (cMD) 4.15 (0.48) 4.13 (0.29) CHARMM27 (aMD) 3.80 (0.66) 4.65 (0.32) AMBER ffSB14 (GaMD) 5.10 (0.46) 6.17 (0.26)
Table 4.1.2: Mean RMSF values for the Cα-backbone, after aligning with respect to the GAF domain, during 500 ns simulation, cMD. Units are Å.
Cα-RMSF Monomer Dimer
CHARMM27 2.073 2.166
AMBER ff14SB 1.953 1.653
Figure 4.1.1: Evolution of RMSD values for Cαof DrBphP with respect to initial crystal conforma-tion (PDB entry 4q0j), throughout the MD simulaconforma-tions computed with CHARMM27 (purple) and AM-BER ffSB14 (green) force fields using a) monomer form simulated with cMD, b) monomer form sim-ulated with aMD, c) dimer form simsim-ulated with cMD and d) dimer form simsim-ulated with aMD. Points collected every 0.1 ns (transparent lines) and every 5 ns (solid lines).The first 100 ns (in the gray box) are regarded as initial equilibration of the system and therefore not considered in the analysis.
(a) (b)
(c) (d)
Figure 4.1.2: RamachandranΦ/Ψplots for the DrBphP dimer: a) CHARMM 27 cMD, b) CHARMM 27 aMD, c) AMBER ffSB14 cMD, d) AMBER ffSB14 GaMD.
Monomer’s dynamic
No changes of the secondary structure or significant protein unfolding is observed in any of the four simulations. The most significant large scale distortion is the bend of the helix spine connect-ing GAF and PHY domains, which changes the angle of the spine centered at Cαof Leu331, from about 140° to 100°(Figure 4.1.4b). This causes the approach of the PHY domain toward the GAF domain as depicted in Fig.4.1.3. Theα-helix bend happens after about 15 ns and remains, with some fluctuations, in this bent position until the end. As result of this displacement the distance between the centers of mass of the GAF and PHY domains drops from about 36 Å to 34 Å (Figure 4.1.4a).
This bend is reproduced in all four simulations, independent of the MD method and the force field used. The driving force for the displacement of the PHY domain towards the GAF are hydrophobic interactions, mainly involving Met1, Pro180, Pro399, Leu410, Leu411 and Pro482 while Gln403, Gln409 and Tyr479 interact with Asp181, that leads to the closure of the cavity between these two domains thus minimizing the solvent exposure by leaving most of the charged side chains at the surface.
Figure 4.1.3: Structure of DrBphP monomer after 500 ns cMD simulation with CHARMM27 force field (in red) compared with the crystal structure [PDB: 4q0j] (in gray), front and back view. The movement of the PHY towards the GAF domain, pivoting around L331, is highlighted by the green dashed line connecting the center of mass of the two domains (plotted in Fig. 4.1.4(a)).
Figure 4.1.4: Evolution of the Center Of Mass (COM) distance (a) and bend angle centered at Leu331 (b) for the monomer, cMD, computed using CHARMM27 (purple) and AMBER ffSB14 (green) force fields.
DCC analysis of the monomer trajectory shows that in the x direction (Fig. 4.1.5b) the motion of the PHY domain exhibits a strong positive self correlation, represented by the intense red col-ored box, meaning that it behaves practically as a rigid body, except for the tongue region, depicted in Figure 4.1.5 inset 1. Due to its high flexibility and distance from the parent domain, the move-ment of the tongue seems to be uncorrelated with respect to the rest of the PHY domain, showed in the white area (box 1 in Fig.4.1.5b) and in all other coordinates. In the y direction (Figure 4.1.5c) the overall behaviour is very similar to the x direction. In the z direction, however, (figure 4.1.5d) the region of the PHY domain spanning residues 367 to 443 (Figure 4.1.5, box 2) show
anticorre-lated movements relative to the rest of the domain. These residues are mostly solvent exposed and connected by flexible loops, which allows for independent movement.
Figure 4.1.5: Dynamic Cross Correlation (DCC) plots for DrBphP monomer throughout cMD using CHARMM27 after frame alignment with respect to GAF domain. a) Spatial orientation of the protein, b) DCC along the x direction, c) DCC along the y direction, d) DCC along the z direction, e) PHY domain spanning residues 367 to 443, colored according to the DCC values shown in the zoom-in box 2.
Dimer’s dynamic
Similarly the DrBphP dimer did not suffer any relevant change in the secondary structure nor a do-main unfolding in the course of the four MD simulations. Different is the dynamic compared to the monomer. All four MD simulations performed using the dimer structure predict the approach of the two PHY domains toward each other, closing the gap separating the C-terminus of the two molecules. This major structural change causes the bend of theα-helix spine, again at Leu331,
lead-ing to a side twist rearrangement, as shown in Figures 4.1.6. The bend angle of molecule A decreases from about 140° to 125°, and in molecule B from 160° to 125°, showing a partial asymmetry (Figure 4.1.7b). According to the simulations, the distance between the centers of mass of the two PHY do-mains changes from 39 Å to about 33 Å (Figure 4.1.7a). The trigger for this movement are electro-static interactions, mainly between Arg495 and Asp496 of the two chains and between two Glu432 facing each other, that draw an average of 4 sodium ions during the simulation. The sodium ions create an attraction site for the, otherwise repelling, glutamate residues (zoomed view in Fig.4.1.6).
This suggests that the ionic strength and ions content of the buffer are important parameters for computational simulations and experimental setup.
Figure 4.1.6: Structure of the DrBphP dimer, molecule A in blue and molecule B in red, after 500 ns cMD simulation with CHARMM27 (in color) compared with the crystal structure (in gray). Zoom in shows the interaction between Arg495 (orange sticks) and Asp496 (bronze sticks) and the two Glu432 (silver sticks) attracting sodium ions (yellow spheres) that allow close the approach of the two residues.
Figure 4.1.7: Evolution of Center Of Mass (COM) distance for DrBphP dimer (a) according to cMD simulations using CHARMM27 (purple) and AMBER ffSB14 (green) force fields. Evolution of the bend angles (b) of chain A (purple) and chain B (green) of DrBphP dimer from cMD with CHARMM27.
The bend angles where computed considering Cαof Leu331 as vertex.
DCC plots in Figure 4.1.8 illustrate the collective motion of all six domains in DrBphP dimer.
The discernible four red colored blocks in the DCC plot in the x direction (Figure 4.1.8b) are indica-tive of a considerable correlated- (diagonal blocks) and cross-correlated movements (off-diagonal blocks) involving the PHY domains of molecules A and B (Figure 4.1.8, box 1), as they move like rigid bodies. The tongue region is uncoupled from the parent domain as well, as reflected by the white stripes in the DCC plot. In molecule B, the tongue exhibit a slightly anticorrelated movement with respect to the PAS domain (Figure 4.1.8b, box 2). Otherwise the tongues show no correlation, due to their flexibility and unique placement. A similar behaviour is predicted in the y direction, (Figure 4.1.8c) albeit with lower correlation intensities between PHY domains and significantly higher anticorrelated motion between PAS and PHY domains of molecule B (Figure 4.1.8c, box 3).
Such dynamic behavior results from the opening of the cavity formed between PHY-GAF domains within one molecule and the closure of the gap between PHY domains of sister molecules. In z di-rection (4.1.8d), the highlighted box 4 shows a mix of both positive and negative cross-correlations, similarly box 5 shows positive as well as negative dynamic correlation within the PHY domain. This region of the PHY domain spanning residues 885 to 932, involves mostly solvent exposed residues which are connected to rest of the protein by flexible loops (Figure 4.1.8e). The dynamic behaviour of this region differs from that of chain A, showing some anticorrelation, similarly to the DCC plot in the monomeric structure.
Figure 4.1.8: Dynamic Cross Correlation (DCC) plots for DrBphP monomer throughout cMD using CHARMM27 after frame alignment with respect to GAF domain. a) Orientation of the DrBphP dimer.
b) DCC along the x direction. c) DCC along the y direction. d) DCC along the z direction. e) Zoom of the box 5 and the correspondent residues (885 to 932) colored according to the DCC values.
Interestingly this kink resemble the one present in a recent DrBphP construct linked to a Syne-chocystisadenylate cyclase (PaaCΔC) [97], Figure 4.1.9, also observed at Leu331, which suggest that accelerated MD can sample conformations that are not otherwise achievable in the absence of the output domain.
Figure 4.1.9: Kink of theα-helix spine of the final DrBphP dimer structure according to GaMD (AM-BER ffSB14), in red, compared to theα-helix spine of the crystal structure of the PaaCΔC construct, in gray (PDB entry: 6FHT).
The bend of theα-helix spine in DrBphP dimer was noticed in the same region in previous stud-ies confronting the structures of Pr and Pfr states [14,98,99]. While, Takala et al. deduce a kink of the helical spine at Val318 and Lys319 on the basis of solution X-ray scattering (SAXS) experi-ments [99], Burgie et al. observed a bow of the spine at Ala326 in the crystallographic structure of DrBphP [98].
Conclusions
While the main dynamic feature in the monomer is the movement of the PHY towards the GAF domain, pivoting around Leu331, a contrasting behaviour is observed in the dimer in which the two sister PHY domains interact with each other, closing the gap that separates them.
Similar movements of the PHY domains were reported forSynechocystissp. PCC6803 Cph1 during Pr→Pfr transition, but the opening motion was regarded only as a transient state in the photoin-duced structural evolution[100].
The kink in the backbone spine was reported, though less pronounced, also for other phytochromes, such as inXanthomonas campestris[25],Pseudomonas aeruginosa[86] andRhodopseudomonas
palus-tris[85]. Bow of the spine and subsequent reorientation of the PHY domains is thought to be a relevant functional step for signal transduction to the output domain [25,98]. There is a certain variability concerning magnitude and timescale of the opening/closing motion of PHY domains, ranging between 1 and 3 nm and in the micro- and milli-second. Although this effect is experi-mentally detected in various phytochromes [14,99–101] there is no complete agreement on its physiological function.
The importance of the dimeric construct, for the photocycle and signaling, was already experimen-tally reported [6,83–85,99]. Although, an essential dimerization interface is located between the GAF domains, as already discussed in previous studies on DrBph [6,14,85] and forXanthomonas campestrisbacteriophytochrome [25], electrostatic interactions between sister PHY domains are also relevant for stabilizing a closed dimer conformation.None of these effects can be observed dur-ing monomer only simulations, which emphasize the importance to conduct MD simulations based on the dimer in order to address biologically relevant large scale dynamic.
Other influential parameters are ionic strength and composition of the solvent that may influence the delicate dimeric interface interactions, both in computational and experimental studies, as can be seen in the formation of the Glu432 salt bridge mediated by sodium ions.
In summary, this study highlight the importance of considering the dimer phytochrome structure when performing atomistic MD, since the resulting thermal equilibrated structures and overall large scale dynamic is significantly different, based on results independent on method or force field used.
4.1.2 Knotless Phytochrome variant
Since resolution of the very first crystal structure of DrBph the knot structure sparked interest due to its unusual nature. Many hypothesis for its role have been advanced, from PAS-GAF rigidity en-hancer to N-terminus restrain[11], but despite the initial speculations, very little has been done to asses its function, but few studies of the effects of its presence[102]. Partly due to the difficulty to efficiently express and crystalize, so called, knotless constructs, i.e. a PAS-GAF photosensory do-main containing BV in which the knot is absent and the two lobes are linked without the interlacing N-terminus. One reason for this is that one of the few constructs that was expressed precipitated is solution, making impossible any further analysis. A reason may be traced back to the increased flex-ibility and freedom of the N-terminus, that entangling with other monomers’ N-terminus, favours aggregation. This already could suggests a function of the knot, to limit N-terminus mobility and aggregation.
While experiments are still challenging, building a model and combining it with MD calculations provides valuable information on the structural role of the knot. Here we present consequences and
While experiments are still challenging, building a model and combining it with MD calculations provides valuable information on the structural role of the knot. Here we present consequences and