Phycocyanobilin in Solution - A Solvent Triggered Molecular Switch

We present a computational investigation of the conformational response of phycocyanobilin (PCB) to the solvents ability to form hydrogen bonds. PCB is the chromophore of several proteins in light harvesting complexes. We determine the conformational distributions in different solvents (methanol and hexamethylphosphoramide HMPT) by means of ab-initio molecular dynamics simulations and characterize them via ab-initio calculations of NMR chemical shift patterns. The computed trajectories and spectroscopic ﬁngerprints illustrate that the energy landscape is very complex and exhibits various conformations of similar energy. We elucidate the strong inﬂuence of the solvent characteristics on the structural and spectroscopic parameters. Speciﬁcally, we predict a cis-trans isomerization of phycocyanobilin upon switching from the aprotic to the protic solvent, which explains an experimentally observed change in the NMR patterns. In the context of technological molecular recognition, solvent induced conformational switching can be considered a precursor mechanism to the recognition of single molecules.


Introduction
Phycocyanobilin (PCB) is an open chain tetrapyrrole chromophore, that can be found in the biliproteins allophycocyanin and phycocyanin.Similar open tetrapyrroles can be found in a vast variety of phytochrome photoreceptors as well as in phycobiliproteins that are found in various plants as well as in different bacteria and fungi.The general function of those photoreceptors is the reaction to light with the goal of collecting energy or triggering further processes in the larger biological system.The process is initialized by an isomerization of the chromophore, which then starts the propagation of the absorbed energy toward the photocenter or triggers the systems structural response on the excitation.While the exact microscopic processes of the latter two steps are still under discussion, the structure of the phytochromes as well as the isomerization process is well documented by different experimental and theoretical investigations.On the experimental side, structural analysis has been performed based on X-ray spectroscopy and neutron scattering experiments.Spectroscopic properties as well as the dynamics during the photocycle have been investigated using various techniques, among them ultrafast mid infrared spectroscopy [1][2][3] as well as Raman and NMR [4][5][6] spectroscopy.Theoretical investigations have been performed on a wide range of aspects employing various theoretical approaches reaching from classical molecular dynamics simulations, as well as QM/MM hybrid approaches.Ad-ditionally, various theoretical spectroscopic approaches have been used, including Raman 7 , UV/VIS via TDDFT [8][9][10][11] as well as NMR simulations 12 , i.e. giving insights into the different chromophore conformations of the single states during the photocycle.In many cases, only the isolated chromophore has been considered, modeling the missing protein structure either by PCM models or ignoring it completely.
In addition to the investigations on the whole protein or photosystem, both experimental and theoretical investigations have also been performed on the isolated chromophore.Naturally, the conformational distribution of the chromophore alone is difficult to assess with diffraction techniques.However, magnetic resonance methods, in particular NMR, are highly sensitive to small changes in structural parameters such as dihedral angles.
The conformational space of a chromophore embedded in the binding pocket of a protein framework is highly restricted.In contrast to that, an isolated chromophore in solution can adopt any configuration, which often leads to an extended distribution of structures.Here, a computational analysis provides a thermodynamically weighted ensemble of geometries from the molecular phase space; this ensemble can subsequently be validated by calculations of spectroscopic observables which can be compared to experimental data [13][14][15][16][17][18][19] .
Previous theoretical investigations on the conformational behavior of isolated similar chromophores have been performed employing Monte Carlo conformational sampling methods 20 as well as direct dihedral scanning 21 .
In a recent experimental NMR study on the isolated chromophore, performed by M. Röben et al. 22 , two fully assigned spectra of the isolated PCB in two different solvents have been obtained.In both 2D 1 H - 15 N spectra, we can observe typ-

Physical Chemistry Chemical Physics Accepted Manuscript
Fig. 1 Structure of the fully protonated (charge +1) phycocyanobilin with the 4 pyrrole rings labeled A to D, connected by 6 dihedrals labeled 1 to 6.The conformation is described by the three double (cis/trans) and three single (syn/anti) bonds between the pyrrole rings ical patterns.The B and C rings chemical shifts are in close proximity in both chemical shift directions, while the A and D rings shifts keep their relative distances of around 1 ppm in the hydrogen and 30 ppm in the nitrogen direction.When changing the solvent from the aprotic HMPT to the protic methanol, a general small downfield shift for the hydrogen chemical shifts can be observed.Additionally, the B and C rings nitrogen chemical shifts are strongly shifted downfield.On a first look, the hydrogen chemical shifts seem to indicate stronger hydrogen bonding for the HMPT solvent, while actually stronger hydrogen bonding can be anticipated for the methanol.The nitrogen chemical shifts suggest specific conformational changes, which on the first look also cannot be predicted from a solvent change.
The possible conformational changes as a response to the presence of specific molecules in the direct environment are a specific form of molecular recognition.In our case, it is not a specific key-lock machanism, but rather a modification of the average hydrogen bonding patterns.Nevertheless, this mechanism can be considered a special variant of molecular recognition [23][24][25][26] .
We investigate said system by a combination of free energy, molecular dynamics simulations and ab-initio NMR chemical shifts calculations to gain insight into the microscopic origin of these on the first look contradicting findings.
Since the nature of the system allows for a huge variety of different conformations and possible hydrogen bonding patterns, in a first step the underlying energy surface in regards the characteristic geometric variables of the system has to be explored.In a next step, we will perform molecular dynamics simulations based on the previously identified most important conformations to finally compute ensemble averaged NMR chemical shifts of the chromophore.Using experimen-tal data as benchmark, we will identify the preferred conformations, molecular structure and hydrogen bonding networks depending on the local surroundings of the chromophore.

Computational Details
Well tempered 2d metadynamics simulations 27,28 have been employed using the 3 dihedral pairs (1,2), (3,4) and (5,6).The simulation was based on a DFTB 29 approach including dispersion correction.The timestep was set to 0.5 fs.A canonical ensemble has been used, utilizing the canonical velocity rescaling (CSVR) thermostat 30 .Gaussian hills with a width of 20 degrees and an initial height of 0.001 Ha have been spawned every 200 timesteps.The temperature of the tempering algorithm was set to 3000 K.
The metadynamics algorithm allows the scanning of collective variables, which is performed by modifying the underlying potential energy while running a molecular dynamics simulation.This modification is achieved by depositing Gaussian hills in regular time intervals.By this, frequently visited regions of the collective variables energy surface are less likely to be visited again in the simulation, energy basins are slowly filled, enabling the molecular dynamics simulation to escape minima and explore regions of the energy surface that could otherwise not be reached on the limited simulation time and temperature.The energy surface can then be reconstructed from the deposited Gaussians.
For the additional sampling of the various found minima, 75 ps DFTB molecular dynamics simulations with dispersion correction 31 have been performed with a timestep of 0.5 fs and employing a canonical ensemble using the Nosé-Hoover 32 thermostat.The same settings have been used for the calculations including the explicit solvent molecules.Here 212 methanol molecules has been added in a fully periodic system with a box size of 25 Å.
We employed a combined gaussian plane wave method (GPW) 33 using the BLYP functional 34,35 with additional van der Waals dispersion 36 for performing the ab-initio molecular dynamics simulations, in order to incorporate hydrogen bonding effects as accurately as possible.We used a time step of 0.4 fs for the integration of the equations of motion, in combination with a Nosé-Hoover thermostat.All MD and metadynamics simulations have been performed using the CP2K program package.
QM/MM molecular dynamics simulations have been performed on two selected conformations, based on a BLYP/ForceField mixed approach.The coupling of the two regions has been done within the gaussian expansion of the electrostatic potential (GEEP) 37 NMR chemical shift averages have been calulated based on regular snapshots (50 snapshots per trajectory) from the molecular dynamics simulations.All calculations have been

Physical Chemistry Chemical Physics Accepted Manuscript
performed using the ORCA program 38 within the IGLO approach 39 .The PBE0 40 functional as well as the optimized IGLO-III 41 basis set have been used.For the explicit solvent calculations, solvent molecules residing within 3 Å of the chromophore have been considered.

Free Energy Surface & Preferred Conformations
The overall structure of phycocyanobilin can be described by means of the dihedral pairs connecting the four pyrrole rings A to D (see fig. 1).The determination of the equilibrium conformation for a given surrounding solvent thus requires a careful sampling of this conformational space.Since the energy barriers of the partly double bonded and possibly conjugated electronic systems are considerable, the molecule will cross them only rarely.Hence, the application of straight molecular dynamics simulations would result in an insufficient sampling or even in stable conformations for the entire simulation.A wealth of methods has been developed to overcome such sampling problems, and have proven to effectively escape local minima.Among them are the parallel tempering algorithm 42 , transition path sampling 43,44 , various basin hopping methods [45][46][47] , minima hopping 48,49 and different tabu search based strategies 50,51 .Other algorithms have been inspired from selective mechanisms found in nature, such as artificial bee colony algorithms [52][53][54] , evolutionary algorithms 55 and genetic algorithms 56 .
Recent studies on a different hydrogenation state of the PCB used Monte Carlo Sampling 20 , while others generated conformations by changing all degrees of the molecules dihedral angles 21 to tackle this problem.In this work, we used the technique of metadynamics simulations, which is an alternative approach that allows the explicit representation of the free energy hypersurface with regards to selected collective variables.This in turn enables us to detect local minima and the corresponding energy barriers between physically meaningful conformations.
We first performed three independent 2D metadynamics simulations of the chromophore in vacuo, each using one of the three pairs of dihedral angles connecting the rings A/B B/C and C/D as collective variables.For all three calculations, a linear (ZEZsss) conformation has been used to initialize the simulation.During the simulations for each dihedral pair, all doubly bonded dihedral angles not driven by the metadynamics remained stable throughout the simulation.Thus, the resulting energy landscape is only valid for the given conformation of the unbiased dihedral angles.However, the perturbance of the energy lancdscape by changes in the electronic structure in distant parts of the chromophore will only have a minor influence.
The resulting free energy surfaces are illustrated in fig. 2 for the three metadynamics simulations sampling the dihedral angles between A/B, B/C and C/D.In the resulting energy surfaces, four minima for each dihedral pair are located at (0, 0), (0, π), (π, 0) and (π, π).All four minima appear twofold degenerate; this effect can be traced back to steric hindering of the amide hydrogen atoms and the pyrrole rings sidegroups.The actual depth of the energy minima however cannot be explained as easily and delicately depends on possible steric repulsion, hydrogen bonding and dihedral energies.
Examining the energy hypersurface of the dihedral pair (1,2) connecting ring A and B, the previously described twofold degenerate minima can be observed at all four conformations leading to planar arrangement of the rings.The minima are separated by two clearly distinguishable energy barriers, a small barrier of approximately 10 kcal/mol in direction of 2 located at ± π 2 and a larger barrier of 30 kcal/mol in direction of dihedral 1 at a similar position.This result clearly shows the different bond characteristics, resembling a double bond character for the central bond of dihedral 1 and single bond characteristics for dihedral 2.
For the two central dihedral angles 3 and 4 connecting the rings B and C, again we can see four distinctive minimum energy conformations that again correspond to the planar conformations.Again, the local maxima within the energy basins are caused by the proximity of different groups and sidechains of the B and C rings.Compared to the previous dihedral pair however, a clear distinction in the bond character and rotational energy barriers for the two dihedrals can not be made.The height of the energy barriers lies in-between those of the rigid dihedral 1 and the rotatable dihedral 2. With a height of approximately 12 kcal/mol, the barrier is small enough to be overcome in the time-limit of actual experiments, resulting in a delicate equilibrium for the two lowest minima of similar energy ((0,0) and (0,π)) with the actual population ratio critically depending on the exact energy difference and external influences.
The remaining two dihedral angles behave similar to the dihedrals connecting ring A and B. The dihedral 5 can rotate easily, while in direction of dihedral 6 a high rotational barrier can be seen.However, here the energy barrier is even higher than in the first case.Also we see a clear energetic preference for a value of 0 deg for dihedral 6, with the lowest possible energy being found at (0,0).
As a main result from our metadynamics simulations, we see similar energetic properties for the different configurations of the (3,4) dihedral pair.The (π,π) and (-π,-π) states however are less favorable due to strong steric hindering.For the dihedral pair (1,2), a preference for the (π,0) conformation can be seen, while for the (5,6) pair the (0,0) state is energetically favorable.This information already heavily reduces the number of needed conformations for a sufficient sampling of the conformational space using molecular dynamics simulations.Additionally, induced by the delicate equilibrium between the two lowest energy minima, the central dihedral pair (3,4) can be expected to react strongly to small changes in the chromophore and its surrounding.Based on the energy surfaces of our metadynamics simulations, unbiased molecular dynamics simulations of the PCB have been employed to allow for a sampling of the unperturbed system.Based on the same DFTB level of theory, we performed simulations for a broad selection of energetic favorable conformations.In 75 ps of simulation time, a good sampling of the relevant regions in the energy surface can be achieved, i.e. positions in both sides twofold degenerate minima are visited, while the general conformation does not change during the simulation.This behavior can be observed for all other conformations and dihedral pairs.As a result, a multitude of differently initialized molecular dynamics simulations is needed to calculate the properties of the system.When evaluating the dihedral distributions during the simulation, the bond characteristics found in the metadynamics simulations can be confirmed (cp fig.3).While the outer dihedral angles 1 and 6 show a narrow double bond like distribution, the neighboring bonds 2 and 5 exhibit a broad distribution indicating a single bond like behavior.The width of the central dihedrals 3 and 4 lies in between the previous ones, pointing towards a conjugated electronic system.While the semi-empirical level of theory allows for the needed sampling, intramolecular hydrogen bonding however is underestimated; the rarely emerging hydrogen bonds only remain stable for parts of picoseconds.

Unbiased Molecular Dynamics Simulations
In a next step, a proper description of the hydrogen bonds needs to be achieved.Here a BLYP with the VdW correction level of theory has been used, which strongly improves the description of the intramolecular hydrogen bonds.Maximally hydrogen bonded conformers can be found in the ZZZsas and ZEZsss conformations, enabling the formation of 6 stable hydrogen bonds between the B and C rings carboxyl groups and the C and D or A and B rings NH and oxygens.

Explicit Solvent Calculations
For the conformational distribution of a solute in an aprotic solvent with only weak specific intermolecular interactions, the gas phase can usually be considered a very good first approximation.For other solvents, in particular hydrogenbonding liquids such as water or alcohols, an explicit representation of the solvent molecules is normally the method of choice in order to reproduce the sensitive equilibrium of intraand intermolecular interactions.Hence, we have considered the solvent molecules explicitly in our molecular dynamics The higher complexity of explicit solvent models results in considerably higher computational costs for metadynamics simulations which require several nanoseconds of simulation time even for an approximate convergence.An additional side effect is caused by the inertia of the explicitly moving solvent, which slows down the effective response time of conformational changes in the solute.This in turn requires a more conservative choice of specific metadynamics parameters (in particular either the amplitude and width of the spawned potential energy hills or the deposition rate).Again, this issue will result in longer simulation times.
When looking at the overall dynamical behavior of the chromophore, again no conformational change can be observed on the computationally accessible timescales.Judging from the dihedral angle distribution, there are no strong changes in the structure of the energy surface.Also the bond lengths suggest no strong changes in the electronic structure.
When analyzing the overall potential energy however, we see a clear change toward lower energies for the closed compared to the open conformation.This behavior can be explained by the formation of intermolecular hydrogen bonds between the chromophore and the surrounding solvent molecules.While the closed conformations in vacuo showed no intramolecular hydrogen bonds during the 75ps simulation time, we now observe an average of two to four intermolecular bonds to the solvent for the different conformations.In contrast, the open conformation was able to form two to six intramolecular hydrogen bonds in the in vacuo simulation, while in methanol the possible hydrogen bonding sites are blocked by intermolecular hydrogen bonding partners.Thus, in a solvent that forms strong hydrogen bonds, the energetic advantage of the open conformation vanishes.
In a conclusion we can assume a different probability for the open and closed classes of conformations, depending on the polarity of the surrounding solvent molecules.The higher the polarity of the solvent, the more intermolecular hydrogen bonds will counter the energy advantage of intramolecular hydrogen bonds.

NMR chemical shifts
In comparison to the recent NMR chemical shift measurements mentioned in the introduction, we perform ensemble averaged NMR chemical shift calculations on the different low energy conformations obtained from the metadynamics simulations.In the experiment, a full assignment of PCB NMR chemical shifts in the HMPT solvent has been achieved 22 .Additionally, a 2D NMR spectrum for the PCB in methanol has been provided by the group of P. Schmieder.While the spectra show a similar shape and ordering, on the first look contradicting changes can be identified.When taking the actual values of the changes between the two different solvents (see fig.

Physical Chemistry Chemical Physics Accepted Manuscript
has to be explained.Based again on the most promising conformations taken from our metadynamics calculations, we perform ensemble averaged 15 N as well as 1 H NMR chemical shift calculations (see tab. 1).When looking at the 15 N chemical shifts of the outer rings A and D, we can see the shifts vary only within a range of 5 ppm for all energetically favorable conformations, indicating only a minor distortion in the local electronic structure upon conformational changes.Generally, the chemical shift of ring D is the lowest, A the highest, while those of rings B and C are always within a range of 5-10 ppm to each other.This overall behavior already is in good agreement with experiments not only in solvent but also with previous experiments on the whole protein.Since the experimental nitrogen shifts are this similar even for hugely different systems or more specifically different surrounding hydrogen bonding networks, one can already suspect the main influence on the nitrogen shifts to be the conformation of the system, while the hydrogens on the other hand are well known to be strongly influenced by hydrogen bonds.
When changing the conformation by moving between the two lowest energy basins of the dihedrals 3 and 4, we see a change in the rings B and C 15 N chemical shifts by approximately 10 ppm.When we look at changes on the chemical shifts induced by rotation around the outer double and singlebonds, we observe a number of small changes in the range of 5 ppm in different directions for the chemical shifts of the central rings B and C.
By a combination of all these possible permutations, a multitude of different, within the sampling accuracy not clearly distinguishable distributions of chemical shifts arises.This quite ambiguous result however can be refined by using the previously attained information about the underlying energy surface and consider only the spectra of the energetically most probable conformations, as well as already comparing the large influence of the central dihedrals to the experimental data.
When considering explicit methanol solvent molecules on the same level of theory, we again see the pattern and changes upon change of the conformation known from the in vacuo calculations.The B and C rings nitrogen chemical shifts are in close proximity and move by approximately 10 ppm upon change of the central dihedrals, while the A and D rings positions remain stable.As a result from the hydrogen bonding to the solvent, the hydrogens chemical shifts slightly move downfield.However, on the DFTB level of theory combined with the range cutoff in the NMR calculations, those interactions are weakened, resulting in a lower hydrogen chemcial shift than expected.As the main result we notice that the influence of the surrounding explicit solvent can mainly be seen in the hydrogen chemical shifts, while the nitrogen shifts stay close to the in vacuo results, however within their usual broad distribution and thus higher possible error from low sampling rates.In no case an influence of the same strength as that of the change of the central dihedral angles on the nitrogen chemical shifts can be observed.

Effect of hydrogen bonding
To properly include hydrogen bonding effects in the suspected main conformations for the HMPT case, BLYP calculations with VdW correction have been performed.In each of those simulations the system has been prepared in a maximally hydrogen bonded state.Here, the hydrogen bonds of the carboxilic groups of the rings B to the C and D ring and C to the A and B ring remained stable during the simulation.While the experimentally observed spectrum cannot be reproduced by a single of those conformations, the influences of the single conformations can attribute for different parts of the spectrum.Combining all of them finally gives an explanation of the experimentally observed chemical shifts.
A "sufficient" sampling of the dihedral energy surface would require long trajectories in the range of 50 to 100 ps which is not unfeasible with present-day computational resources, but appears unnecessary in view of the required computational cost.On the other hand, examining short trajectories of different hydrogen bonding patterns can also give insights into the resulting spectra and possible origins of the different characteristics for the two solvents.
One low energy conformation is the EEEsss configuration.Here however, for geometric reasons only one stable intramolecular hydrogen bond can be observed.The result-  For the maximally hydrogen bonded ZEZsss conformation, we see a strong change in almost all rings positions in the spectrum.However, with exception of the only weakly hydrogen bonded B ring, the alignment of the chemical shifts is similar to the non hydrogen bonded case with only the hydrogen chemical shifts moved downfield by approx. 2 ppm.Most importantly, the non-or weakly hydrogen bonded B rings nitrogen chemical shift is strongly shifted downfield.This shows clearly, that the position of the nitrogen chemical shift is not induced by the hydrogen bonding, but has instead to be caused by conformational effects and deviation from the energy minima of the dihedral distributions found in the metadynamics simulations.

Conclusions
We have quantified the dependence of the conformation of a tetrapyrrole chromophore (phycocyanobilin) on the character of its chemical environment, specifically the solvent type.By means of semiempirical metadynamics simulations, we have sampled the free energy hypersurface of the chromophore, resulting in the identification of the geometric equilibrium structures at ambient temperatures as a function of the solvent.In particular, our simulations show a delicate equilibrium between helical and linear conformations that reacts very sensitively on the solvent characteristics.Based on ab-initio molecular dynamics simulations, we have computed the twodimensional NMR pattern of 15 N and 1 H NMR chemical shifts at the ab-initio level of theory.
The experimentally observed solvent shifts exhibit a counterintuitive trend (towards higher 1 H NMR chemical shifts for HMPT which is expected to form fewer hydrogen bonds than methanol).Our ab-initio molecular dynamics simulations have resolved this puzzling situation, which could be traced back to stronger intramolecular hydrogen bonding of the chro-mophore.These hydrogen bonds are broken in methanol due to the better electrostatic screening and the stronger intermolecular interactions with the solvent.
A striking secondary effect of the change in hydrogen bonding pattern is the isomerization of the chromophore which is equally visible in experimental as well as computed chemical shift pattern.This solvent-induced structural transition represents a mechanical switching function on the molecular scale.With a suitable chemical functionalization, this transition might be exploited for applications in the field of molecular sensors, molecular recognition or even externally triggered drug delivery.

1- 8 | 3 PageFig. 2
Fig. 2 Energy surface for the three dihedral pairs (1,2) to (5,6) connecting the pyrrol rings A to D. Twofold degenerate minima are visible for planar conformations, while the overall depth of the minima and heigth of the energy barriers varies to a great extent.

Fig. 3
Fig. 3 Dihedral angle distributions during a molecular dynamics simulation.All dihedral angles remain stable on the simulation time of 75 ps.The different bond characteristics are clearly visible.

Fig. 4
Fig. 4 Distribution of the four possible hydrogen bonds between the NH groups and the carboxilic groups lone oxygens, taken from a 10 ps BLYP-D trajectory of the ZEZsss conformation.Top left: ring A, right: ring B, bottom left: ring C, right: ring D

Fig. 5
Fig.5Experimental 2D 1 H and 15 N NMR spectra different solvents: HMPT22 (blue dots) and methanol57 (red squares) of the amine groups of rings A to D

5 Page
5) into account, we see that the A and D rings δ N change by 2 ppm, the δ H by 0.6 ppm.For the B and C ring, we see changes of δ N by 10 ppm and δ H by 1.2 ppm.Assuming the change is mainly induced by hydrogen bonds and exhibits a linear behavior, i.e. upon change in hydrogen bond strength the chemical shift values move on a straight line in the 2D spectrum, an additional change of the δ N by around 6 ppm 1-8 |

Fig. 6
Fig. 6 Overview of the different influences on the 1 H and 15 N chemical shifts of PCB.A conformational change in the central dihedral pair of ring B and C causes a change in the positions of the signals of ring B and C as indicated by the green and blue areas.The remaining signals A and D remain stable.Changes in hydrogen bonding strength mainly causes the hydrogen chemical shifts to change.

Table 1 1
15and15N NMR chemical shift values for a series of characteristic conformations.1Hreferenced to benchmark calculations of TMS,15N referenced to match the PCBs D ring in methanol.