Molecular Mechanisms of Activation in the Orange Carotenoid Protein Revealed by Molecular Dynamics.

Light-harvesting in photosynthesis is accompanied by photoprotective processes. In cyanobacteria, the photoprotective role is played by a specialized complex, the orange carotenoid protein, which is activated by strong blue-green light. This photoactivation involves a unique series of structural changes which terminate with an opening of the complex into two separate domains, one of which acts as a quencher for the light-harvesting complexes. Many experimental studies have tried to reveal the molecular mechanisms through which the energy absorbed by the carotenoid finally leads to the large conformational change of the complex. Here, for the first time, these mechanisms are revealed by simulating at the atomistic level the whole dynamics of the complex through an effective combination of enhanced sampling techniques. On the basis of our findings, we can conclude that the carotenoid does not act as a spring that, releasing its internal strain, induces the dissociation, as was previously proposed, but as a "latch" locking together the two domains. The photochemically triggered displacement of the carotenoid breaks this balance, allowing the complex to dissociate.

Introduction

Bacteria are the oldest type of photosynthetic organisms on earth. This makes their investigation even more important as they can be seen as ancestors of modern organisms and they could help us to understand how the photosynthetic process has appeared and evolved over time. Most photosynthetic bacteria have developed anoxygenic forms of photosynthesis, but there is one type, the cyanobacteria, which performs photosynthesis in a manner similar to that of plants using water as an electron donor and releases oxygen. Cyanobacteria have developed a photoprotection system that prevents photodamage when they are illuminated with more light than they can use for photosynthesis, exactly as in plants and algae. Photoprotection necessarily competes with the efficient collection of light energy which is the basis of photosynthesis, and as such, the system has to be switchable and there must be a reversibility between the so-called light-harvesting and photoprotective states. In the latter, the excess energy has to be rapidly dissipated into heat through what is generally known as nonphotochemical quenching (NPQ).[1−4] In cyanobacteria, NPQ is rather different with respect to plants and algae, where the same antenna complexes which perform light harvesting can switch to the quenching function.[5,6] The antenna of cyanobacteria are in fact phycobilisomes (PBSs), peripheral membrane complexes that incorporate several types of bilins, which are open-chain tetrapyrrole pigments. They are organized into disks which create large rod-like architectures attached radially to a core protein.[7] Excitation energy is funneled along the rods to the core which finally can attach to either photosystem I or II. NPQ occurs at the core level and involves a water-soluble protein called orange carotenoid protein (OCP).[8−10]

OCP can be found in two interchangeable conformational forms: orange (OCPO) and red (OCPR).[11,12] In the former, the embedded keto-carotenoid cofactor is placed in a pocket crossing the two domains (N- and C-terminal domains, or NTD and CTD) which constitute the protein, and it is firmly connected to the CTD domain through two hydrogen bonds with two protein residues (Figure ). This closed orange form has to switch to the red one to become active in NPQ: this activation requires the absorption of strong blue-green light. The final product (OCPR) is obtained through a translocation of the carotenoid within the NTD domain and a dissociation of the two domains which remain connected by a long, flexible (largely unstructured) loop linker.[12−19] The NTD (with its embedded carotenoid) can then interact with the PBS core, allowing the quenching of the bilins.[20] The red form slowly back-converts to the orange one in the dark, but the recovery can be accelerated through the interaction with another protein, the fluorescence recovery protein (FRP).[21,22]

Figure 1

Graphical representation of the OCP:CAN complex from its crystal structure (4XB5).[12] Colors are used to show different domains (magenta for NTE, yellow for NTD, gray for the linker, blue for CTD, and green for CTT). CAN and some relevant residues (Y44, W110, Y201, W277, and W288) are shown as licorice.

Experimental techniques have probed the evolution of OCP after the photoactivation with high resolution in either time[18,19,23,24] or space.[25] Nonetheless, a clear and complete characterization of the molecular mechanisms through which the carotenoid translocates between the two domains and the protein dissociates is not available.

The goal of this article is to reveal such mechanisms through atomistic simulations. This goal poses a major challenge as very different time scales are involved, from the nanosecond to microsecond for the carotenoid to translocate within the protein, to the millisecond for the domain dissociation.[18] Moreover, the light-driven nature of the process should require the use of a quantum-mechanical dynamical description that at the present time can deal only with very short time scales.

To achieve such a challenging goal, we exploit a combination of different enhanced sampling techniques[26−29] to overcome the initial energy barrier and to efficiently sample the equilibrium configurations visited by the CAN during the translocation, after the initial photochemical event is concluded.

Using this paradigm, we could reproduce the whole evolution of the complex, from its resting form (OCPO) to its final open form. The simulations demonstrated that the carotenoid acts as a “latch” locking together the two domains in OCPO, by establishing interactions of comparable strength with each of them. This latch, however, is lifted as soon as the carotenoid is allowed to move toward the NTD thanks to the electronic transition, and the dissociation can proceed. Indeed, with our simulations, we showed that only after carotenoid translocation can the dissociation of the two domains take place. Moreover, we could identify two relative free-energy minima in the translocation process that are related to the intermediates revealed by recent spectroscopic studies.[18] Finally, we revealed the fundamental role played by water in both the translocation (by stabilizing the intermediates) and the dissociation (through solvation effects), which we showed to occur, leaving the secondary and tertiary structures of the two domains almost unchanged.

Results

To achieve a clearer description of the results of the simulations, we divided their presentation into two parts, translocation and dissociation. Finally, we recollected all of the main findings of both parts in a general discussion. In the latter, a direct comparison with available experimental evidence is also reported.

Translocation

In our simulation of the translocation process, to enhance the sampling of the reaction path we performed an umbrella sampling (US)[30] simulation. At each point of the collective variable, to sample perpendicular degrees of freedom, instead of a plain MD, we performed a replica exchange molecular dynamics (REMD) simulation.[31] The combination of umbrella sampling and REMD (US-REMD) should allow for an efficient exploration of the transition space along both the main direction and the slow degrees of freedom orthogonal to the collective variable, such as structural rearrangement of the protein or conformational changes of the carotenoid. The collective variable (CV) chosen to represent the translocation process is the distance between the center of mass of the carotenoid and that of the NTD. A detailed discussion of the methods used for the simulation and their analysis is provided in the SI.

To asses the convergence of the simulation, we performed blocking of the data set, and for each block, we computed the free-energy surface (FES) along the CV (details in the Supporting Information). The analysis shows that the second half of the simulation has completely converged; therefore, we analyzed the structures (about 0.8 million) extracted from the second half of the 10 lowest-temperature replicas of our simulations and calculated the FES of the system along the CV. The resulting profile is reported in Figure C. To provide a first interpretation of the predicted free-energy minima, for each of them we calculated the average structure and the average RMSD distance of the NTD and CTD using OCP and RCP crystal structures[12] as references (Figure A,B,D).

Figure 2

(A) Average structure at CV ≈ 11.5 Å (in blue) superimposed on the crystal structure of OCPO (in orange) and (B) at CV ≈ 2.5 Å superimposed on the RCP crystal structure (in red). (C) Free-energy profile along the CV (in kcal/mol). (D) Average RMSD values (in Å) for the NTD computed using as reference the NTD of the crystal structures of OCPO (orange line) and RCP (red line) and for the CTD computed using as reference the CTD of the crystal structure of OCPO (blue line). Comparison of representative structures of OCPO and (E) T1 and (F) T2 states. A licorice representation is used for CAN and for the side chains of some relevant residues spotted either in the literature (particularly the two H-bonding Y201 and W288 and the π-stacking Y44, W110, and Y111)[32] or in this present work (the interface H-bond network E244, R155, N104, and W277). In order to give a better view of the complexes, we used a different orientation in panels A and B (roughly rotated by 180° around the horizontal axis).

In a previous work by our group, 1-μs-long unbiased MD simulations of OCPO and RCP were performed.[33] Considering that the selected CV, on this unbiased OCPO trajectory, has an average value of 11.4 Å, we can safely state that the deepest minimum (highlighted in orange in Figure C) corresponds to the OCPO state. Indeed, by comparing the RMSD values of the two main domains and the three-dimensional structure of the absolute minimum overlapped with the reference OCPO structure, it is clear that this first state is almost identical to the one known from OCPO crystals (Figure A).

Moving to the left of the FES, following the translocation of the carotenoid into the NTD, two other relative minima appear at CV ≈ 7.5 Å (highlighted in blue, T1) and at CV ≈ 2.5 Å (highlighted in red, T2).

According to our simulations, T1 represents an intermediate in the reaction path with the carotenoid only slightly displaced from its original position but still enough to lose its H-bonds with the β1 ionone ring and form other stabilizing interactions with other residues of the pocket, such as W41 and P124–P126 (Figure E).

On the contrary, in T2, the position of the carotenoid is strikingly similar to the one found for CAN in RCP (Figure B). We stress that the CV here used does not contain any information on the exact binding site in which the carotenoid should reach into the NTD. This finding confirms that the “tunnel” already spotted from the X-ray diffraction structure[12] acts as a guide to drive the carotenoid to a new position, which is very similar to the binding mode observed in RCP.

During the process, the NTD undergoes small but significant modifications. These can be quantified by looking at the RMSD of the NTD along the evolution of the CV when calculated with respect to OCPO and RCP (Figure D): while the former increases during the process, the latter decreases. It is important to stress here that these results are completely independent of any a priori knowledge of the structure of RCP and that they rely only on the crystal structure of OCPO and on the carotenoid–NTD distance (i.e., very general and unspecific information).

Let us now consider the RMSD of the CTD along the CV from the reference CTD of OCPO. This domain is in general more stable than the NTD during translocation, as it roughly shows a single jump from 1.25 to 1.75 Å at CV ≈ 6.5 Å. Visual inspection reveals that this change in the RMSD can be connected to a small conformational change localized in the highly conserved[34] T275–N281 loop, which will be discussed in more detail later. On the basis of these observations, we suggest that T2 is a reasonable representation of the final product of the translocation.

So far, we have explored the evolution of the structure of the system during the translocation process. Now we analyze the same evolution for the interactions among the various components (carotenoid, protein, and solvent).

As a first analysis, we dissected the average energy of the carotenoid in three parts: the interaction with the protein, the interaction with the solvent (i.e., all water molecules), and its internal energy. The values of each part along the CV are reported in Figure . We note that the energies reported there do not include the internal energy change of the protein and any entropic contribution to the free energy. As can be seen from Figure A, the solvent plays a central role: immediately after the first displacement of the carotenoid from its original binding pocket, strong favorable interactions appear that in T2 largely compensate for the observed decrease in the stabilizing effects of the protein (Figure A, green line). The latter behavior can be explained in terms of the loss of the two hydrogen bonds with Y201 and W288 and the appearance of interactions between the water molecules and the polar groups on the β-ionone rings of CAN.

Figure 3

(A) Interaction energies between CAN and OCP (green line) and CAN and water (blue line), internal energy of CAN (orange line), and total interaction energy (black line): all values are calculated with respect to OCPO. (B) Interaction energies between CTD and CAN (blue lines) and NTD and CAN (yellow lines); interactions among the CTT, linker, and NTE are negligible; dashed lines are the electrostatic interactions, dotted lines the Lennard-Jones interactions, and solid lines are the sum of the two contributions. (C) Electrostatic (ES) and Lennard-Jones (LJ) interaction energies among individual residues of NTD, CTD, and CAN. Values are reported as the difference from the OCPO state. Numerical values for ES and LJ interactions of CAN with the most relevant residues of the protein are reported in Tables S2 and S3, respectively.

If we analyze more in detail the profile of the interactions of CAN with the protein, then we see complex behavior. After the initial destabilizing effect due to the loss of electrostatic effects from the broken H-bonds, we see a first minimum roughly corresponding to T1 and a second shallow one corresponding to T2. To understand more clearly, we can further dissect the interaction energy into nonelectrostatic, namely, repulsion and dispersion as described by the Lennard-Jones (LJ) potential, and electrostatic (ES) contributions from the two main domains (Figure B). Contrary to the chemical intuition, in OCPO the carotenoid is bound somewhat more tightly to the NTD than to the CTD, despite the two H-bonds. Therefore, it is not surprising that, when the carotenoid is placed outside the H-bonding range of W288 and Y201, it is much more strongly bound to the NTD than to the CTD (about 20 kcal mol–1). This large effect has a strong contribution from new LJ interactions with the NTD and partially also from breaking those with the CTD.

The evolution of the LJ and ES interaction energies of CAN with single residues along the CV is reported in Figure C. Many residues in the NTD enhance their LJ interactions with CAN, passing from OCPO to T1. The most relevant ones are W41 (−2.1 kcal mol–1), R155 (−2.0 kcal mol–1), I125 (−1.5 kcal mol–1), and P124 (−1.2 kcal mol–1) (Figure E). This seems to indicate that the translocation process is not driven by a few specific residues but instead the whole region of the protein around the carotenoid contributes. It is also interesting that after the loss of the two H-bonds, no other strong electrostatic interactions arise in state T1.

As the translocation proceeds further and the carotenoid has moved deeper into the NTD, its LJ interactions with this domain keep increasing. In particular, we observe that a group of residues has a large stabilization effect (E34, L37, W41, T80, N104, W110, Y111, I125, P126, and Y129). The set of these residues defines a mainly apolar/aromatic binding pocket for the carotenoid. As expected, two of the tightest interactions between CAN and NTD that were present in OCPO state (the ones with W110 and R155) are considerably looser in T2 because of the new geometry of the complex in which the carotenoid is rearranged in a new binding site (Figure F). Moreover, two noticeable electrostatic interactions appear in T2. They can be classified as (relatively loose) H-bonds with the highly conserved W277 and T275 residues belonging to the CTD. The formation of these two H-bonds is connected to the largest local conformational change observed in the protein, corresponding to a rearrangement of the T275–N281 loop (Figure S9), which is one of the constituents of the NTD/CTD interface. To understand the importance of this conformational change, we should note that in the RCP crystal structure, as in those obtained from previous MD simulations,[33] a stable H-bond between R155 (involved in a salt bridge between NTD and CTD) and N104 is observed. In OCPO and in T1, N104 is stably bound to W277; after the rearrangement of loop T275–N281, this last H-bond breaks up and N104 can establish a new stabilizing long-range interaction for R155.

We finally note that the β1 ionone ring can form H-bond interactions with various residues belonging to both NTD (G57 and S60 (OH in the side chain)) and CTD domains (W277 (NH in the side chain and NH in the backbone) and T275 (OH in the side chain)). As T2 is expected to evolve into a dissociated state, it is reasonable to expect that the carotenoid could interact with different sites and obtain a global stabilization effect (clearly visible in Figure A), without forming very tight connections that could inhibit the further evolution of the system.

As already pointed out, water plays a central role in stabilizing intermediates T1 and T2. It is well known that the β1 carbonyl of CAN in OCPO is stably bound with Y201 and W288 in a solvent-excluded region of the protein, and it does not form any relevant interaction with water. However, as the system evolves toward T1, the number of waters that can form H-bonds with CAN starts increasing (Figure S4A). This is likely due to the combination of different effects: the rupture of the two stable H-bonds of CAN with the protein residues, the motion of CAN toward a more solvent-exposed region, and also the entrance of water in the pocket of the CTD that initially was occupied by the CAN. We investigated the latter process by calculating the change in the solvent-accessible surface (SAS) for each residue along the CV (Figure S4B in the SI).

Many CTD residues increase their solvent accessibility during the process. With few exceptions, this change takes place starting from T1. This is an indication that the first H-bonds between β1 carbonyl and solvent are formed with water molecules that are already present within the protein pocket. In the second part of the translocation process (from T1 to T2) instead, a great increase in the SAS of CTD residues is observed together with a significant global decrease in the SAS of NTD residues. This fact suggests that the motion of the carotenoid in the protein tunnel displaces the water molecules present in the NTD and simultaneously allows them to enter in the CTD cavity, probably through the NTD/CTD interface.

As a final important analysis of the translocation energetics, we note that the gain in the internal energy of the carotenoid from OCPO to T2 is small, and at T1, the CAN has a larger internal energy than in OCPO (Figure A). On the contrary, the interactions with solvent and protein provide a substantial stabilization of the carotenoid in T2 with respect to OCPO. All of these findings contradict a previously proposed hypothesis[35] that the distortion of the carotenoid acts as a spring that pushes the system from OCPO to OCPR.

To reach a more complete analysis of the role of internal geometrical changes, we have investigated the conformational freedom of the terminal rings of the carotenoid during the translocation. To do so, we have plotted the FES on the CV vs β1 and β2 dihedral angles (Figure A,B) and extracted some relevant conformations of the carotenoid from each minimum. It should be noted that since the torsion constants of the force field have been fine-tuned against DFT data, the terminal rings are normally able to rotate on the nanosecond time scale (Figure S8 and SI of ref (33)). The conformational rigidity observed in some states of the system is therefore determined by intermolecular interactions between the carotenoid and protein/solvent and does not depend on the force field parametrization itself.

Figure 4

Plot of the free-energy surface projected on the (A) CV-β1 and (B) CV-β2 subspaces obtained from reweighting of the US-REMD simulations. (C) Overlap of 10 configurations of CAN in T2, T1, and OCPO together with a single representation of the protein environment as a spatial reference. The upper panels show the β1 side of the carotenoid, and the lower ones, the β2 side.

Looking at β1, it is clear that this angle remains locked into an s-trans-like conformation until the system reaches T2. This conformation is the same as the one seen in the OCPO crystal structure (Figure C). No other relevant conformation of this dihedral angle is seen until the completion of the translocation, consistent with the picture of Konold et al.[18] where the photoproducts before translocation have the β1 ring in an s-trans conformation. β2 behaves similarly, but since in T1 this ring of the carotenoid is already outside the pocket defined by Y44 and W110, it has a larger conformational freedom that allows it to switch between the two s-cis conformations. Once the system is in T2, both terminal rings are able to rotate with remarkable freedom, as the disorder found in the configurations extracted from the MD clearly demonstrates. Nonetheless, we can notice a slight preference of β1 for s-cis conformations and of β2 for s-trans conformations. This fact is coherent with the picture of a very different conformational population of the carotenoid in RCP and OCPO, as already shown by our previous simulations and confirmed by subsequent experimental findings.[33,36]

Dissociation

There is general agreement on the fact that the dissociation is by far the slowest process in the photoactivation of OCP. This step is expected to have large entropic contributions (due to the larger conformational freedom of the two dissociated domains in OCPR) that compensate for the loss of binding interactions between minor and major interfaces. In light of these observations, we used well-tempered metadynamics (wt-MetaD) to enable a more efficient exploration of the region relevant to the dissociation process. The choice of the CV for wt-MetaD is based on the experimental evidence that upon dissociation both the minor and major interfaces of OCPO break up. We therefore defined a contact map, restricted to the specific contacts relevant to the interfaces: by definition the CV is zero when both major and minor interfaces are formed and reaches one when the dissociation is complete. (A rigorous definition of this variable is provided in the SI.) The wt-MetaD simulations were run by biasing both the radius of gyration (Rg), which is expected to increase, and the interface contact map.

We compared two equivalent wt-MetaD simulations starting from a representative structure of OCPO and T2 states. This comparison will allow us to understand if the translocation of the carotenoid facilitates the dissociation and the changes in the interactions that are responsible for this opening.

The results of the two simulations are collected in Figure : panels A and B show the configurations of the system sampled in each wt-MetaD simulation projected in the subspace of the biased collective variables, and panels C and D show the estimated free-energy surfaces for the dissociation process.

Figure 5

(A, B) Configurations of the system sampled in the wt-MetaD simulation projected in the subspace of the biased collective variables: the radius of gyration (Rg) and the interface contact map (CM). (C, D) Free-energy surfaces for the dissociation process projected in the subspace of the biased collective variables. Plots A and C refer to the simulation started from OCPO, and plots B and D refer to the one started from T2.

At a qualitative level, we note a clearly different behavior in the two simulations. Starting from OCPO we do not see any dissociation event even if the bias accumulated during the wt-MetaD simulation is very large (about 80 kcal mol–1); we recall that this bias can be seen as an artificial lowering of the free-energy barriers encountered by the system along the CVs. Since under these simulation conditions the system remains trapped in a set of configurations very close to the original minimum, the well-known stability of the closed form of the system is confirmed. On the contrary, the T2 state dissociates rather quickly: it explores a transition region where the two interfaces are only partially formed, and then it reaches the dissociated region (CM ≈ 1) where the two interfaces are broken and the system is able to move almost freely in the Rg dimension.

This result is not unexpected from what was reported in the previous section. According to our results, in fact, the interaction of the carotenoid with each of the domains of the protein (roughly 40 kcal mol–1) was almost perfectly balanced in OCPO and was fundamental to the stabilization of the closed form of OCP. Once the carotenoid is translocated into the NTD, the interaction energy with the CTD is largely reduced (less than 10 kcal mol–1) and the dissociation process becomes much easier, as the only barrier for the process is the cleavage of the inter-residue interactions at major and minor interfaces. From a thermodynamic point of view, we expect that the energy loss due to the disruption of the two interfaces can be compensated for by the solvation of the interfaces themselves and from the entropy gain of the dissociated state.

According to this model, the interactions of the carotenoid with NTD and CTD are crucial to the stabilization of the “closed”, globular form of the protein. To verify this hypothesis, we performed an additional wt-MetaD simulation under the same conditions detailed above, but starting from OCPO without the CAN molecule inside the cavity (details in the SI). In this case, the two domains of OCP are able to dissociate as they do in the simulation starting from T2 (Figure S5). This fact further proves that the dissociation is mainly enabled by a change in carotenoid–protein interactions and not by a modification of the protein structure induced by the photochemistry of the carotenoid itself.

Unfortunately, simulating transitions with such a large entropy production remains an extremely challenging task even when using advanced enhanced sampling techniques. In fact, even when extending the simulation starting from T2 up to 1 μs, no back-conversion is observed. Since the ergodic regime was not attained, the simulation cannot be used to estimate the relative free energy of the closed and open states. In spite of these limitations, our wt-MetaD simulations clearly show that the dissociation of the two domains follows the translocation of the carotenoid, confirming the hypothesis of the most recent working model for the photoactivation process.[18,24,25] Moreover, it becomes clear that this order of events is a consequence of the fact that the modification of the carotenoid–protein interactions caused by the translocation plays a great role in enabling the subsequent dissociation.

To gain information on the dissociated state, we analyzed the data collected after the dissociation event which involve several transitions from small and large Rg.

In particular, we focused on the structural stability of the different domains of the protein. To do so, we calculated the RMSD of structured domains over the wt-MetaD trajectory. Furthermore, for each structure extracted from the trajectory, we computed a structure assignment using the DSSP algorithm.[37,38] Considering the secondary structure data (Figure A), we note that the two NTD and CTD domains are quite stable; this is also confirmed by the fact that the RMSD of these domains remains around 2 Å for the whole simulation (Figure B) and by the direct observation of the structures extracted from the simulation (Figure C,D). A minor instability in CTD is observed in the region corresponding to the T275–N281 loop (Figure A,D). This loop, in the globular form of the protein, is localized at the interface between NTD and CTD, and its structure is already influenced by the translocation process. Therefore, it is not surprising that, upon dissociation, it can experience increased conformational flexibility. On the other hand, the structure of the two minor domains, namely, CTT and NTE, seems to be deeply influenced by the dissociation: the RMSD of these two domains already indicates a notable structural change with respect to the globular form of the protein (Figure B). Moreover, looking directly at the structures extracted from the simulations, it seems that the relative orientation of NTD and NTE (represented respectively in yellow and pink in Figure C) is not fixed after the dissociation and the NTE behaves as a disordered coil. This finding suggests that the formation of the minor interface is essential to the stability of the structure of the NTE.

Figure 6

(A) Secondary structure of OCP calculated on the frame extracted from wt-MetaD simulation. (B) RMSD of NTE (purple), NTD (yellow), CTD (blue), and CTT (green) for the structures extracted from wt-MetaD. Each RMSD was calculated using OCPO as a reference and aligning the domain to the reference structure.Superposition of 80 frames of the wt-MetaD simulation, showing NTE and NTD with (C) CAN inside and (D) CTD. (E) Sample structure of the full dissociated OCPO with the α-helix already formed in the linker region.

As a final remark, we note that the relative position of the carotenoid inside the NTD remains almost unchanged with respect to the T2 state, which closely resembles the binding site of RCP, as shown in Figure . This fact strongly supports the validity of RCP as a model for the photoactivated form of OCP.

A last interesting fact that arises from the analysis of the dissociated state is that a new α-helix motif appears in the region of the unstructured linker of the two domains (Figure E, with the linker region represented in gray). This motif seems to be quite stable during our simulation, suggesting that, even if it corresponds to a low conservation zone of the protein, it could have biological relevance.

Discussion

In our simulations, we cannot see the ultrafast photochemical steps after carotenoid excitation. However, irrespective of the mechanism involved, these steps lead the complex to a first intermediate where the hydrogen bonds are broken. This is an out-of-equilibrium state that is thermally not accessible, which should resemble the geometry of T1. The nonequilibrium population of this state should be able, at least in part, to overcome the barrier and evolve toward T2. On the contrary, if the system was equilibrated in state T1, then, considering the relative heights of the barriers found in our simulations, we would observe a fall back to OCPO. This is what happens during OCPR → OCPO thermal reversion, which follows a different path with respect to the photoactivation.[17,25]

The most recent experiments also support the fact that the translocation process precedes the dissociation of the NTD and CTD domains by orders of magnitude.[18,24,25] These findings are fully confirmed by our results. In fact, our wt-MetaD simulations show that the T2 state has a remarkable preference to undergo dissociation with respect to OCPO. In addition, our results indicate that the carotenoid interactions with the NTD and CTD are finely balanced in the dark-adapted state (Figure B), stabilizing the closed OCPO form. Similarly, carotenoids were shown to be essential to the stabilization of cyanobacterial PSII.[39] Instead, as soon as the carotenoid is displaced along its conjugation axis, toward the NTD, an imbalance in these interactions appears, slackening the binding of the two domains (Figure B). This change in the interaction energy with both domains allows them to break apart in the ensuing entropically driven thermal process.

This picture is in contrast with the model of a dissociation process pushed forward by a released strain of the carotenoid.[12,35] Indeed, according to our results, the internal energy of the carotenoid only slightly changes during the translocation, showing that the geometrical relaxation of the chromophore cannot be the driver of the process. Instead, here we propose that after excitation the carotenoid relaxes back to its ground state inside the cavity of an out-of-equilibrium system resembling the T1 state. In such a state, the structure of the protein is almost unchanged, the two H-bonds are already broken, but the carotenoid is still less than 5 Å away from its original binding site. Part of this population will then evolve into an equilibrium state (T2), where the carotenoid is in its final location and the NTD has rearranged, together with the T275–N281 loop. This hypothesis allows us to link the calculated FES found for the translocation process with recent experimental studies.[11,18,24] In particular, by focusing on the model proposed by Konold et al.,[18] we note a remarkable resemblance between our T1 state and the P2 proposed by the authors. In addition, the global rearrangement of the NTD and the loop that we observe along the translocation, between T1 and T2, could be responsible for the changes observed in the time-resolved IR spectra of the complex 0.5–10 μs after the photoexcitation.[18] The structure of the NTD shows several changes (Figure D) which reflect more than one rearrangement taking place along the translocation. Finally, it is easy to link our T2 state to the P3 intermediate spectroscopically characterized by the same authors.

We note that the largest component of the carotenoid–protein binding energy is nonelectrostatic in character (Figure C). This is not unexpected, as it was noted that many highly conserved residues with hydrophobic/aromatic side chains are found in the carotenoid binding pockets of both OCPO and RCP.[34] The importance of these interactions in the photoactivity of the protein has been demonstrated by several mutagenesis experiments. In fact, mutating aromatic residues, even absolutely conserved ones, around the binding tunnel preserve the photoactivity of the system as long as these are substituted by other aromatic residues.[12,34] On the other hand, all NTD single mutations affect the relative stability of OCPO/OCPR to some extent, enhancing the thermal recovery rate and reducing the photoactivation yield. A combination of these mutations, which taken alone are not able to switch off the OCPO photoactivity, makes the protein nonphotoresponsive.[12] This evidence, together with our results, leads us to propose that the photoactivity of OCP and its ability to thermally reform the orange state are the results of the fine-tuning of electrostatic and nonelectrostatic interactions of the carotenoid with the two main domains.

Another important open question that cannot be easily answered, even by exploiting the most advanced experimental techniques, is which structural changes of the protein are induced by the translocation. Whereas it is possible to gain general knowledge about the identity of structural rearrangements from IR studies, the insight obtained is far away from atomic-level resolution. On the other hand, XF-MS experiments, which have an intrinsically high spatial resolution, cannot provide time-resolved data at the required accuracy level.[12,14,24,25] Our model, on the contrary, provides a detailed atomistic description of the structural changes of the protein along the whole process. Surprisingly, localized changes in the NTD, such as helix unfolding, were not observed in our simulation. Instead, we found a global rearrangement of this domain, mostly affecting the relative orientation of helices αC, αD, αE, and αG and of the N-terminal side of αI (using the nomenclature proposed in ref (34)). In the CTD, we observed a localized rearrangement of the T275–N281 loop that connects the β3 and β4 strands. These changes in the tertiary structure of the protein allow for a rearrangement of the H-bond network at the interfaces of the two domains (Figure ). Moreover, the change in side-chain orientation of some residues, such as W277, provides stabilizing interactions for the carotenoid along the translocation. Remarkably, a very recent paper indeed found that the mutation of W277 with a phenylalanine results in a nonphotoactive complex, while the W277H mutation does not disrupt the photoactivity.[40] This result indirectly supports the importance of residue W277 not only for π-stacking interactions with the carotenoid moiety but also for its H-bond donor properties that allow its interaction with the carbonyl group of the carotenoid itself.

In our simulations, without any a priori bias in this direction, we observe a change in the NTD that makes it more similar to the crystal structure of RCP (Figure B,D). However, our simulations do not see the “bending” of helices αC and αD which is evidenced by the comparison of crystal structures of OCPO:CAN and RCP:CAN. This discrepancy can be explained as either an artifact of the crystal structure or a change that takes place after the dissociation of the CTD. To support the first hypothesis, we also note that the two crystal structures of the RCP:CAN complex (PDB IDs 4XB4 and 5FCX) show slightly different structures in αC and αD helices (Figure S7).

Moving to the dissociation step, our data provide a clear indication of the relative stability of the different domains of the protein after their separation. As expected, the structures of NTD and CTD do not change significantly upon dissociation, as their average RMSD values remain at around 2 Å. The main instability is located in CTD, and corresponds to a relative flexibility of the T275–N281 loop (Figure D). This loop was already noted to play an important role during the dissociation and is an active constituent of the NTD/CTD interface. It is therefore not surprising that its structure is less stable after the unbinding of the NTD.

Upon dissociation of the NTE, this polypeptide chain is not able to maintain its helical structure and appears disordered both in its secondary structure and in its relative position with the NTD. This result is indeed expected from IR light-minus-dark data provided by Konold et al. for WT-OCP and ΔNTE-OCP.[18] Our results support the fact that the presence of the minor interface is essential to stabilizing the helical structure of NTE.

On the other hand, we noted that the linker has a certain tendency to form an α-helical structure after dissociation. We think that the presence of this structural element in a region that has always been thought to be unstructured in both OCPO and OCPR should be considered when interpreting IR spectra. Our data do not allow a quantitative assessment of the stability of this helix. This fact, combined with the low conservation of the primary structure of this region, suggests that a more in-depth study is needed to asses the importance of the aforementioned folding in the actual OCPR solution structure.

Conclusions

By simulating for the first time, with atomistic techniques, the evolution of OCP after the photoactivation event has taken place, this study has shed light on the many still unclear features of this unique light-driven process, clarifying its most important structural and energetic aspects.

We found that only after carotenoid translocation could the dissociation of the two domains take place. Moreover, the latter process leaves the secondary and tertiary structures almost unchanged. We identified two relative free-energy minima in the translocation process, namely, T1 and T2, that are related to the intermediate products revealed by recent spectroscopic studies.[18]

The atomistic detail of our simulations allowed us to understand the reason for the reduced or absent photoactivity of OCP mutants recently reported in the literature.[40] In fact, we highlighted the crucial role played by W277, which stabilizes the T2 state via H-bonding to the carotenoid, and by many aromatic residues in the NTD cavity. Mutations of these residues are expected to alter the stability of the intermediate, making it short-lived and more likely to return back to OCPO. Therefore, nanosecond transient absorption of NTD mutants, such as W277F, would reveal the spectral response of the early translocation intermediates only. Comparing these transient absorption spectra among the mutants and with the wild type will then help to identify the intermediates along the photoactivation path.

On the basis of our findings, we can conclude that the carotenoid locks the two domains together in the resting form of the complex (OCPO) by establishing comparably strong interactions with both of them. The photochemically triggered displacement of the carotenoid breaks this balance, allowing the complex to dissociate. In this model, the carotenoid does not act as a “spring” that, releasing its internal strain, induces the dissociation, as was previously proposed.[35] Instead, the carotenoid seems to act as a “latch” that, when moving into the NTD, releases the binding interactions between the two domains, allowing them to split up.

Methods

We model OCP starting from the final structure of the 1 μs plain dynamics described in our previous work[33] and based on the high-resolution crystal structure of CAN-binding OCP (PDB: 4XB5).[12] All simulations were run with GROMACS 2019.4 suite using the classical MM parameter from the ff14SB AMBER force field.[41] CAN was described using modified GAFF parameters reported in a previous article.[33] Water was described using the TIP3P model. Unless differently specified, all of the simulations were run using a 2 fs integration time step, together with the LINCS algorithm to constrain every bond between a heavy atom and a hydrogen. A cutoff was applied after 12 Å on the LJ interactions, while long-range electrostatics was treated with a particle mesh Ewald algorithm. The simulation were performed in the NVT ensemble using the modified velocity rescaling thermostat,[42] applied to the whole system, with a coupling constant of 0.1 ps. Periodic boundary conditions were applied in the three dimensions.

US-REMD Simulations

We performed 11 REMD simulations, each restrained with a different umbrella potential as described in Table S1. Each REMD simulation was performed with 50 replicas in the temperature range of 300–400 K (details in the Supporting Information). This allows for a significantly improved exploration of the orthogonal degrees of freedom. The harmonic potential on the reaction coordinate was applied using PLUMED 2.5.[43] Each replica of the REMD simulation was first minimized and then equilibrated in the NVT ensemble (100 ps) and NPT ensemble (100 ps). The production simulation was then started. An exchange attempt was performed every 500 steps (1 ps). Average exchange probabilities of around 0.1 between neighboring replicas were attained (details on replica exchange results in the SI). Every 5 ps we saved the structure of the system for further analysis. Simulation lengths are reported in the SI. To generate the initial structures for the US-REMD simulations, a short ratchet-and-pawl MD (rMD) run was performed.[44,45] In an rMD simulation, the system evolves freely whenever it spontaneously approaches the target state (defined according to a suitable CV). Conversely, a harmonic biasing force switches on any time the system attemps to backtrack to the initial state. This simulation was started from a structure of the OCPO complex (taken from a previous unbiased simulation[33]) equilibrated in an appropriate dodecahedron solvent box (14 763 water molecules, 37 Na+, and 30 Cl– needed to obtain a neutral system with a salt concentration of 0.1 M). The rMD simulation was run for 10 ns. The harmonic ratchet force was applied to the same CV as already described for the US-REMD simulation (the distance between NTD’s and CAN’s centers of mass). The target value was set to 1.0 Å, and the force constant was set to 10 000 J mol–1 nm–2. The closest frame to the center of the umbrella potential was used as the initial configuration for the subsequent simulation.

wt-MetaD Simulations

wt-MetaD simulations were performed using PLUMED to apply the bias.[43,46] We used σ = 0.05 on the CM dimension and σ = 0.2 on the Rg dimension. The hill height was set to 0.5 kcal/mol, and the γ factor was set to 15. A hill was deposited every 2000 steps (4 ps). To speed up the bias calculation, the potential was calculated on a grid spanning from −0.1 to 1.1 in the CM dimension and from 10 to 60 Å in the Rg dimension; 50 points were used in each dimension. Two simulations were started from representative structures of OCPO and T2. Both of these structures were extracted from the previous US-REMD simulations. In order to avoid self-interaction of the system once it underwent dissociation, a larger cubic box was used with respect to the US-REMD simulations. To do so, we re-equilibrated the system in a larger cubic box (20 Å from the box walls instead of about 10 Å used before) composed of 40 850 water molecules and 86 Na+ and 79 Cl– ions. Before wt-MetaD was started, the system was equilibrated with harmonic restraints on the initial position of the heavy atoms of the complex in the NVT ensemble (100 ps) and in the NPT ensemble (100 ps) and finally for 10 ns without any restraint in the NVT ensemble. To further guarantee that self-interactions were avoided, harmonic upper and lower wall potentials were used to limit the system in a reasonable region of the phase space on the Rg dimension. The upper wall was placed at 45 Å, and the lower one was placed at 15 Å. Both harmonic potentials have a force constant of 5 kcal/mol Å–2.

Analysis of the Trajectories

An analysis of the MD trajectories was performed using in-house python scripts. Data manipulation was performed by exploiting numpy[47,48] and scipy[49] libraries, and graphics were generated with matplotlib[50] libraries. Direct manipulation of the MD trajectory (e.g., extraction of relevant configurations, calculation of geometrical properties, and interaction energies) was performed directly with either GROMACS, PLUMED,[43] or the MDAnalysis library.[51] To efficiently solve the MBAR/UWHAM equations in the analysis of the US-REMD trajectory, we used the fastMBAR[52] library, which allows us to use a modern GPU-accelerated computer to solve these equations with limited computational cost. Utilities from the PyEMMA[53] library were also used for trajectory binning and CV discretization. All graphical representations of the system’s structures were rendered with the open source version of PyMol.