Carbon Atoms Speaking Out: How the Geometric Sensitivity of 13C Chemical Shifts Leads to Understanding the Colour Tuning of Phycocyanobilin in Cph1 and AnPixJ.

We present a combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics-statistical approach for the interpretation of nuclear magnetic resonance (NMR) chemical shift patterns in phycocyanobilin (PCB). These were originally associated with colour tuning upon photoproduct formation in red/green-absorbing cyanobacteriochrome AnPixJg2 and red/far-red-absorbing phytochrome Cph1Δ2. We pursue an indirect approach without computation of the absorption frequencies since the molecular geometry of cofactor and protein are not accurately known. Instead, we resort to a heuristic determination of the conjugation length in PCB through the experimental NMR chemical shift patterns, supported by quantum chemical calculations. We have found a characteristic correlation pattern of 13C chemical shifts to specific bond orders within the π-conjugated system, which rests on the relative position of carbon atoms with respect to electron-withdrawing groups and the polarisation of covalent bonds. We propose the inversion of this regioselective relationship using multivariate statistics and to apply it to the known experimental NMR chemical shifts in order to predict changes in the bond alternation pattern. Therefrom the extent of electronic conjugation, and eventually the change in absorption frequency, can be derived. In the process, the consultation of explicit mesomeric formulae plays an important role to qualitatively account for possible conjugation scenarios of the chromophore. While we are able to consistently associate the NMR chemical shifts with hypsochromic and bathochromic shifts in the Pg and Pfr, our approach represents an alternative method to increase the explanatory power of NMR spectroscopic data in proteins.

1. Introduction

Tetrapyrroles showcase an overwhelming abundance in living species, accompanied by an impressive functional versatility [1]. They can be classified into cyclic and open-chain forms, the former comprising the well-known porphyrins that build the precursors of hemes and chlorophylls, renowned for their important role in oxygen transport, oxidoreductase enzymes, and photosynthesis [2]. The second group, open-chain tetrapyrroles, also referred to as bilins, are ubiquitous in living species as well. Organised in biliproteins, they take the role of a chromophore and serve phototrophic organisms as accessory pigments [3]. An important class of biliproteins is formed by phytochromes—photosensors that typically entail histidine kinase activity to regulate photomorphogenesis and photoperiodism [4,5,6,7]. The occurrence of phytochromes was long assumed to be restricted to higher plants, but over the years members of this class have also been found in cyanobacteria, bacteria, and fungi [6,7,8,9,10,11]. Recently, variants of (bacterio)phytochromes were found, genuine to cyanobacteria, that are therefore called cyanobacteriochromes (CBCRs) [12,13]. They fulfil a widespread range of tasks in the cell connected to the response towards light, such as positive phototaxis [12,14,15]. Although being phylogenetically related to phytochromes, CBCRs bear unique features: while the sensory module of phytochromes needs to encompass at least a period circadian/aryl hydrocarbon receptor nuclear translocator/single-minded (PAS)–cGMP phosphodiesterase/adenylyl cyclase/FhlA (GAF)–phytochrome specific (PHY) tridomain (Figure 1), CBCRs require only a chromophore-binding GAF domain to maintain their biochemical properties [16]. Most intriguing, however, is the ability of CBCRs to tune their absorption over the entire range of the visible spectrum, including near-IR and near-UV [12,17,18,19,20]. The treasure of photosensors provided by nature with which organisms adapt to their environment bears high utility for biochemical applications in bioimaging and optogenetics, but also in therapy and catalysis [14,21,22,23,24,25,26,27,28,29,30,31,32,33,34].

Figure 1

(Left) Sensory module of Cph1, a phytochrome of cyanobacterial origin. Phytochromes comprise a PAS (yellow)–GAF (green)–PHY (blue) tridomain, [16] wherein the latter acts as signal transmitter [10,35]. (Right) Photoactive GAF domain (green) of AnPixJ, a cyanobacteriochrome (CBCR) [13]. In both cases, the cofactor phycocyanobilin (PCB) is covalently bound to the GAF domain (red). Note: PAS—period circadian/aryl hydrocarbon receptor nuclear translocator/single-minded; GAF—cGMP phosphodiesterase/adenylyl cyclase/FhlA; PHY—phytochrome specific.

Many cyanobacterial phytochromes and CBCRs carry as bilin cofactor, phycocyanobilin (PCB), which operates as a photoreversible switch, toggling between the ZZZssa and ZZEssa configuration (Scheme 1) [4]. Exciting the chromophore at its absorption maximum, photochemical Z/E isomerisation of the C–C double bond initiates a flip of the D ring, which triggers a biochemical response [6,36,37,38,39,40,41]. The two states, Z and E, are denoted as dark state and photoproduct, respectively. The photoproduct is back-transformed into the dark state by photochemical E/Z isomerisation after excitation at its, typically shifted, absorption maximum. Depending on the colour (x) of absorbed light, the states are labelled as “P”. Phytochromes switch between red-absorbing (P) dark states and far-red-absorbing (P) photoproducts, [4,8] while the colour diversity in CBCRs is broader. Typical examples of CBCR absorption colour pairs are red/green, green/red, blue/green, red/orange, green/orange, violet/orange, blue/yellow, or blue/teal [12,14,18,19,20,34,42].

Scheme 1

Structural formula of phycocyanobilin (PCB) in ZZZssa (left) and ZZEssa (right) configuration. Reversible photoisomerisation of the C–C double bond triggers the photoswitch.

The colour tuninpan>g inpan> CBCRs is exceptional, regarding the very limited number of chromophore species they bind. We recently showed that the degree of solvation due to a sealed chromophore binding pocket cannot serve as an explanation model [43]. However, Peng and co-workers experimentally related the dihedral angle of ring D of PCB in allophycocyanin B (a phycobilisome) with the absorption maxima ranging from 590 to 670 nm and thus proved the bathochromic shift’s originating in an increased conjugation length [44]. Recently, Wiebeler and co-workers showed the hypsochromic shift in the green-absorbing photoproduct of Slr1393g3 (a PCB-carrying CBCR) to emerge from a decrease of conjugation, originating from the out-of-plane rotation of ring D (“trapped-twist” model) [17,45]. They underlined the important role of the apoprotein to keep the chromophore in the intended conformational state.

In this work, we present a quantum mechanics/molecular mechanics (QM/MM) study of PCB in a protein environment by means of extensive molecular dynamics (MD) and ab initio C-nuclear magnetic resonance (NMR) chemical shift calculations, based on density functional theory (DFT). Comprising the entire protein together with the solvent, our theoretical model is able to account for physical meaningful non-covalent interactions without losing accuracy at describing the chromophore and its instantaneous thermal fluctuations. We focus on two cyanobacterial photosensor proteins (Figure 1): the red/far-red switching Cph12 (see “Abbreviations”), a well-studied member of the phytochrome superfamily [10,18,35,38,39,40,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66] and the red/green switching AnPixJg2 (see “Abbreviations”), a cyanobacteriochrome that has been discovered later [12,13,43,45,67,68,69,70,71,72,73]. Both proteins have PCB as a cofactor and exhibit a red-absorbing dark state (P). Photoconversion yields a far-red absorbing photoproduct (P) in the case of Cph12, whereas AnPixJg2 renders a green-absorbing state (P). Although the spectral shift in opposite directions amounts to a difference of more than 100 nm, photoproduct formation neither implies chemical changes to the chromophore nor a changed protonation state [45,70,72]. However, according to the studies by Peng and Wiebeler it can be concluded that the D-ring twist and thus the effective conjugation length determine the absorption properties [17,44]. To date, no crystal structures of the photoproduct states of Cph12 and AnPixJg2 exist that could build the basis of theoretical structure models, but we show that relying on the two P states, structural fluctuations can profitably be correlated with calculated spectroscopic patterns and increase the explanatory power of NMR spectroscopy towards conjugation effects in the P and P state, respectively.

2. Results

2.1. How Bond Alternation and Thermal Fluctuations Determine 13C-NMR Chemical Shifts

Extensive MD sampling under the QM/MM regime allows for the inclusionpan> of thermal fluctuations of the molecular structure on the femtosecond and picosecond time scale. Figure 2a shows the computed time evolution of the C–C bond length and the C chemical shift of atoms C and C—exemplary for all similar features in PCB—as they are found in the simulations. The thermal motion of the atoms at 300 K leads to bond length oscillations in the range 1.4–1.5 Å, and to corresponding oscillations of the instantaneous carbon NMR chemical shifts within a range of about 20 ppm. These oscillations are too fast for the NMR experiment; only their averages (dashed lines in Figure 2a, cf. also Table 1) can be compared to corresponding measurements. Figure 2b matches computed C chemical shifts of PCB’s carbon atom scaffold (C through C) with reported MAS-NMR data from Matysik and co-workers for the P states of AnPixJg2 and Cph12 at 233 K [68]. The simulations clearly reproduce the characteristic alternating pattern of C chemical shifts mirroring the -conjugated system, which only proves that PCB in the model represents the physical state of the chromophore in the proteins. At the same time, though, the uncertainty in the computed averages impedes distinguishing between the two protein environments of the cofactor. In fact, with chemical shifts differing by less than 5 ppm, AnPixJg2 and Cph12 are hardly distinguishable experimentally let alone by the QM/MM model given the large fluctuations shown in Figure 2a.

Figure 2

(a) Thermal fluctuations exemplified by bond length and C chemical shifts of carbon atoms C and C in PCB along a quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulation of AnPixJg2 in its P form (solid line) and running average (dashed line). For atom numbering see Scheme 1. (b) C chemical shifts of PCB in AnPixJg2 (black) and Cph12 (blue). Crosses mark the values computed from QM/MM MD trajectories while full circles indicate the experimental reference (dashed lines to guide the eye) [68].

Table 1

Results of the inversion of the geometric sensitivity of C chemical shifts shown in Figure 3 based on extensive QM/MM MD simulations and reported NMR measurements by Song et al. [68] for AnPixJg2 (switches from P to P) and Cph12 (switches from P to P). (P–P) is the experimental difference pattern of C chemical shifts going from P to P; d(P) denotes averaged bond lengths calculated from the MD trajectories of both proteins in their respective P state; (P–P) is the change in bond length (together with uncertainty) going from P to P, predicted from principal component analysis and multivariate regression.

		C_13–C_14		C_14–C_15		C_15–C_16		C_16–C_17		C_17–C_18		C_18–C_19
	C_13		C_14		C_15		C_16		C_17		C_18
Δδ^c(P_r–P_r) ^†,*	−0.7		+1.2		+3.8		−0.5		−0.9		−2.0
Δδ^c(P_g–P_r) *	−6.1		−7.4		+2.3		−8.4		+4.7		+9.3
Δδ^c(P_fr–P_r) *	+4.3		+6.1		−1.6		+5.7		−6.6		+6.5
Δδ^c(P_g–P_fr) *	−11.1		−12.3		+7.7		−14.6		+10.4		+0.8
d(P_r) **		143		144		137		148		137		149
Δd(P_r–P_r) ^†,**		−1.3±0.5		+4.1±0.6		+2.1±0.5		−3.8±0.6		−3.1±0.4		−1.6±0.6
Δd(P_g–P_r) **		−16.6±2.0		+4.8±2.3		−14.4±1.7		−4.5±2.3		+19.1±1.3		−6.2±2.1
Δd(P_fr–P_r) **		+12.8±1.6		−3.2±1.9		+7.8±1.3		−4.6±1.8		−4.1±1.0		+16.0±1.7
Δd(P_g–P_fr) **		−30.7±3.0		+12.0±3.6		−20.1±2.6		−3.7±3.5		+20.1±2.0		−11.4±3.3

P−P, * in ppm; from reference [68], ** in pm.

Although it seems that the uncertainpan>ty inpan> the obtainpan>ed averages of C chemical shifts suggests a high sensitivity of the carbon atoms’ shielding towards any direct or indirect change in the environment, the further analysis of their dependence on chromophore–apoprotein connections and temperature effects enables a more refined notion of the matter. Computationally, it is straightforward to repeat chemical shift calculations for the same setup, gradually “switching-on” those protein residues that interact with the chromophore, or numerically “cooling-down” the system, which forces the atoms back to their equilibrium positions. Corresponding results can be found in Appendix A (Figure A1). These technical modifications shed light into the susceptibility of C chemical shifts in PCB—and into those of other nuclei. It turns out that they are largely insensitive to temperature and non-covalent interactions the chromophore is subjected to due to their mainly quartary nature. This contrasts with the high sensitivity of H and N nuclei, the chemical shifts of which change strongly with both temperature and supramolecular interactions.

Figure A1

Normalised RMSE of calculated shieldings with increasing size of the QM part (left) and after subsequent annealing and equilibration of FFMD snapshots (right) at the example of AnPixJg2.

However, Figure 2b undoubtedly indicates that the conjugation pattern and charge distribution in the chromophore strongly determine the extent of C chemical shifts. This striking correspondence points to the fact that instantaneous fluctuations like those shown in Figure 2a actually incorporate an intimate connection of bond parameters and nuclear shieldings. How this manifests in PCB will be shown and discussed in the following section; we so far conclude that the chromophore’s carbon atoms C–C, though being blind towards supramolecular effects, are an important probe of their immediate covalent situation.

2.2. Correlations of 13C Chemical Shifts and Their Geometric Sensitivity

The interconpan>nection of nuclear shieldings and bond lengths is not particularly surprising, but the remarkable “isolation” of C chemical shifts from supramolecular interactions deserves a closer look. Figure 3 presents the linear correlation matrix of C shifts with C-C bond lengths in rings C and D of PCB in P state. It is based on plotting the instantaneous values obtained from the QM/MM trajectories against each other and extracting statistically significant connections through linear regression. It therefore reveals whether chemical shifts decrease or increase upon increasing the individual bond lengths (see Figure A2 for the explicit correlation diagram). It is important to note that we are entitled to carry out a combined analysis of both proteins alike, given the P state of both proteins, AnPixJg2 and Cph12, having PCB in the same or a very similar configuration (the matrix shown in Figure 3 looks essentially the same for either protein). It is the central result of this report, because it identifies not only pivotal correlations, but furthermore reveals a particular pattern, which in the following we refer to as the geometric sensitivity of C chemical shifts. There are numerous representations that describe a sensitivity, yet in this study we are interested in those chromophore parts that are involved in colour tuning, that is, the nuclei describing the conjugation between rings C and D (C–C).

Figure 3

Correlation matrix showcasing the geometric sensitivity of C chemical shifts towards changes in selected C-C bond lengths in rings C and D of PCB in AnPixJg2 and Cph12 (P state). Colour and size of the circles correspond to the Pearson correlation coefficient (see colour bar on the right), whereas white numbers denote the slope of linear regression (in ppm Å). Mint cells are significant at ; snapshots.

Figure A2

Explicit correlation plot of C-C bond lengths and C chemical shifts in PCB.

Qualitatively, the geometric senpan>sitivity follows the clear rationale that an atom’s chemical shift increases with the length of the originating bonds—and decreases with the next, and increases again with the following bond and so on. In Figure 3 this alternation is marked by red and blue bullets, respectively, shaping an imagined double diagonal; like Figure 2b, it is another representation of the underlying conjugated system. The symmetry of the geometric sensitivity, however, is clearly disturbed leaving certain nuclei with a “preferred” bond they are correlated with. For instance, there is a strong connection of atoms C and C to bond C–C, but a rather weak one to bonds C–C and C–C, respectively. The same observation can be made for C and C, while the latter shows a strong connection to bond C–C, whereof it is not even a member. Expectedly, the two pyrrole systems mainly exhibit within-ring correlations, but atoms C and C are connected to both parts. Amide C, in turn, appears completely independent from the conjugated system.

The characteristic sensitivity in Figure 3 evidently is determined not only by the bond alternation of the conjugated system, but also by bond polarisation due to the relative position of the nuclei to electron-withdrawing groups and substituents. Nevertheless, the direct interpretation of the geometric pattern is difficult and we suggest employing tools from multivariate statistics. As will be shown in the following section, we can use the geometric sensitivity together with experimental evidence to predict the conjugation length in P and P states of PCB alike, assuming that the correlation matrix in Figure 3 remains unaltered in either state.

2.3. Inversion of the Geometric Sensitivity Points to Pg and Pfr States

In 2015, Song et al. reported a peculiar differenpan>ce pattern of C nuclear shieldings () in PCB upon photoconversion in AnPixJg2 and Cph12: Nuclei belonging to rings C and D exhibit changes in chemical shifts of up to 10 ppm when transforming from the P state into the respective photoproduct [68]. The authors attribute this observation to a “large-scale rearrangement of the critical interactions [that] cannot be explained only by the configurational changes of the chromophore.” However, recalling that C chemical shifts are blind towards supramolecular perturbations (cf. Figure A1) the changes cannot be explained by the removal of hydrogen bonds or the emersion of a new charged side chain—unless this involves changing the covalent situation of the carbon atoms. Most intriguing, however, is the observation by Song of an opposite trend of the pattern for most atoms upon switching to either P or P state with P lying in the middle. Following the correspondence principle, this trend must be found in the causing process (i.e., photoconversion) as well, which is supported by what is by now commonly understood as the colour tuning mechanism: [17,44] the extent of conjugation in rings B, C and D of the chromophore, where the P and P states are characterised by a decreasing and increasing conjugation length relative to the P state, respectively.

The NMR experimental observations are made amid scarce information on the atomic coordinates of the photoproduct states since no structural data from X-ray diffraction are available. In these terms, the geometric sensitivity may serve to retrieve structural data from the C-NMR experiment and to predict the single–double bond conjugation pattern for the P and the P state. To endow our theoretical observations with quantitative information, we inverted the correlation matrix in Figure 3 and applied the result on the patterns reported by Song et al. [68] We carried out a principal component analysis (PCA) to account for the importance of fluctuations in the feature variables before using multivariate regression to obtain predicted patterns that correspond to the experimental patterns, as shown in Table 1 (see Appendix A for computational details). We report as a reference the average bond lengths of PCB in the P state, calculated from the MD trajectories. They support the inferred change in bond order, for instance to what extent a single or double bond character is retained. The given uncertainties of prediction help to distinguish important from less significant trends upon photoconversion.

Comparinpan>g the two P states of Anpan>PixJg2 and Cph12, we do not find big differences between the two proteins, which corresponds to the earlier argument that PCB in both proteins is chemically equivalent. The found numbers are mainly due to the comparably large difference in chemical shift for atom C, lying at the joint of rings C and D and which is connected to many bonds (cf. Figure 3). Going over to predicted changes in bond length upon photoconversion, that is, differences going from P to P and P, respectively, we are able to distinguish key sites among other, less significant changes. Bond C–C, for which we hardly see a difference in the P state, undergoes important changes and exhibits a clearly opposite behaviour in the process: while in P state it is shortened by about 17 pm, its length increases by about 13 pm in P. Bond C–C is already longer in P of AnPixJg2 and becomes even longer in P, but shorter in P of Cph12, amounting to a net difference of 12 pm between the photoproducts. Bond C–C again shows opposite behaviour being shorter in the P state (–14 pm) and longer in the P state (+8 pm). For the bonds forming ring D the prediction tends to assign changes to either one or the other photoproduct, but changes in bond C–C are small and not significant. Bond C–C is clearly lengthened in P (+19 pm) as is bond C–C in P (+16 pm) with respect to P. Consequently, this analysis reveals that it is possible turning the asymmetry of the geometric sensitivity of C chemical shifts to precisely account for subtle changes in the covalent structure of PCB, thereby allowing regioselective predictions. In the following discussion we make the attempt to interpret the pattern on a molecular scale.

3. Discussion

The presented results of the inversion of the geometric senpan>sitivity of C chemical shifts are based on experimental evidence reported by Song et al. [68] We inferred expected changes in bond length of the PCB chromophore when going from the dark state P to the photostates P for AnPixJg2 and P for Cph12. We pointed out—in accordance with the common understanding of the colour tuning mechanism in these forms—the antagonistic nature of the changes involved. The results in Table 1 can be used to predict the predominant character of the investigated C-C bonds of PCB in P and P states. Associating bond lengths with bond orders, we propose the following alternation pattern:

Table 2 suggests that the bond alternation pattern in the two photostates is opposite and we present in the following an explanatory model to translate this into the effective conjugation length. One has to keep in mind that the alternation pattern represents the predominant character of the C-C bonds, but does not compare to pure single or double bonds since pyrrole rings form a -conjugated system not only between, but also within themselves.

Table 2

Predicted C-C bond character trends in the conjugated system of rings C and D of PCB in P and P state based on Table 1 and a threshold of 144 pm. There is no reliable prediction possible for bond C–C.

	C_13–C_14	C_14–C_15	C_15–C_16	C_16–C_17	C_17–C_18	C_18–C_19
P_g	double	single	double	(single)	single	double
P_fr	single	double	single	(single)	double	single

Conpan>siderinpan>g that the shapinpan>g of delocalisation of double bonds by means of mesomeric schemes has so far not received due attention, we argue that working out what is actually possible and favourable in terms of bond alternation is a key to fully understand the colour tuning phenomenon in PCB.

3.1. Pg: Short Conjugation

Scheme 2 and Scheme 3 illustrate the -resonance of PCB in P state where conjugation has been found to be of short length. We consider the conjugation between rings B and C, which is maintained in all states of the chromophore, as core conjugation. Assuming that in P state ring D is twisted out of plane, core conjugation alone carries the absorption of green light, which is in line with what has been reported by Wiebeler et al. [17] There exist several mesomeric formulae describing this conjugation, the most important of which are shown in Scheme 2. Therein, the positive charge is delocalised along with the change of bond alternation, but the high electronegativity of the nitrogen atoms shifts the distribution towards the formation of stabilised tertiary carbocations centered at C and C, thereby rendering bond C–C with an increased bond order. A stabilisation of the positive charge to some extent at the pyrrole nitrogen atoms seems reasonable since both proteins present residues at this location that may serve as hydrogen bond acceptors, but have to compete with water molecules—also in AnPixJg2 [43]. In any case this would not swivel the predicted bond character at C–C.

Scheme 2

Mesomeric formulae showcasing the delocalisation of positive charge in PCB between rings B and C. It represents the core conjugation, which is prevalent in the P state of the chromophore. The investigated bond C–C, marked in red, inherits a double bond character due to the importance of mesomeric forms with a (stabilised) tertiary carbocation.

Scheme 3

Negative mesomeric effect () exerted by the carbonyl group in ring D of PCB. It corresponds with amide resonance, which is prevalent in the P state of the chromophore. The formation of a (stabilised) tertiary carbocation renders bonds C–C and C–C with a single–double alternation.

Simultaneously, ring D, being fully detached from the conjugated core, is subjected to amide resonance pulling electron density from the pyrrole ring and generating a negative partial charge at the carbonyl oxygen atom. However, in pyrrole, resonance does not only include the NH group of the amide bond, but the entire aromatic system. As a consequence another mesomeric form can be found, which is again stabilised through formation of a tertiary carbocation far from the electronegative centres. This has the important consequence that bond C–C no longer retains a double bond character, which is shifted to C–C, but that of a single bond with atom C carrying the positive (partial) charge.

We want to note that the methylene bridges are not favourable sites to put a positive partial charge, due to their secondary type, but they still take part in conjugation. Furthermore, it has been argued that a certain single bond character of C–C is the underlying cause of “dark conversion” [6].

We conpan>clude that a plausible prediction of the geometric bond length conjugation pattern of the PCB in its P state is possible through measured and computed NMR chemical shifts (and their geometric sensitivities). The resulting conjugation pattern, in turn, provides a consistent explanation of the change in absorption wavelength after photoconversion, which is an explanation for the unusual colour change of the bilin chromophore in this protein.

3.2. Pfr: Long Conjugation

-Resonance of PCB in P state follows different principles than in the case for the P state. As was pointed out by Peng et al., the increase in absorption wavelength is due to the extension of conjugation, which now includes rings B, C, and D since the latter rotates in the plane with the core [44]. We shall call this full conjugation because it resembles the longest possible resonance in PCB. Including ring D in the exchange of -electrons has not only bonds C–C and C–C toggle their single–double bond pattern, but also withdraws charge density from the ring as the positive charge is now shared as well.

Scheme 4 presenpan>ts the most important mesomeric formulae in the P state and an immediate question might be why there is no favourable carbocation formation in ring D as has been adduced for rings B and C. The reason behind this is that in the case of P, an explicit interaction with the apoprotein makes a difference. We have argued before that the carbon atoms do not directly “see” supramolecular effects; however, they are sensitive to indirect effects, which in the present case is the stabilisation of the positive charge by the apoprotein: Rohmer et al. argued that in P the negatively charged Asp207 is in position to form a salt bridge with the pyrrole nitrogen atom carrying the positive charge [38]. Based on our analysis, we strongly suggest that this interaction is the dominating factor in the -resonance of the P state. The new stable mesomeric form that can be found in the P state corresponds with the predicted bond alternation pattern in Table 2.

Scheme 4

Mesomeric formulae showcasing the delocalisation of positive charge in PCB between rings B, C, and D. It represents the full conjugation, which is prevalent in the P state of the chromophore. Ring D carrying the positive charge is stabilised by a salt bridge to an aspartate residue (Asp207) of the apoprotein, thereby endowing this mesomeric form with high importance [37,38]; bonds C–C, C–C, and C–C are rendered with single–double–single alternation. The presence of positive charge in ring D also suppresses amide conjugation (cf. Scheme 3), which has bonds C–C and C–C retain their double–single alternation.

In summary, also for PCB inpan> its P state a plausible prediction of the bond length conjugation pattern is obtained, based solely on measured and computed NMR chemical shifts (and their geometric sensitivities). Again, the resulting conjugation pattern serves as explanation of the change in absorption wavelength through photoconversion.

4. Conclusions

We have addressed the colour tuninpan>g mechanism in the photosensor proteins AnPixJg2 and Cph12, which are based on the phycocyanobilin (PCB) chromophore, via a quantum chemical method using an indirect approach. While the direct approach (i.e., the calculation of absorption energies in the different configurations) is possible, but difficult due to the absence of reliable structural data, we have chosen to use reported experimental C-NMR chemical shifts of the chromophore backbone. These data are available for both P and P state, together with the chemical shift changes with respect to the dark state conformation, P, of known structure. We infer structural information and derive the respective bond length alternation pattern, which points towards the effective conjugation in the unknown P and P states. This can eventually be linked to the changes in absorption wavelength via simple confinement considerations.

Our complemenpan>tary approach represents an uncommon computational route that is interesting for situations in which clear and reliable structural information (e.g., from X-ray diffraction) is difficult to obtain due to experimental constraints (here, the lifetime of the photoproducts P and P). The approach requires the quantitative determination and subsequent inversion of a particular structure–property relationship; namely, the dependence of the NMR chemical shifts on the bond length pattern. This can be realised in the form of a matrix representing the transformation of the bond length vector into the NMR chemical shift vector at linear order. For the special case of AnPixJg2 and Cph12, we used extensive QM/MM MD simulations and ab initio computations of C chemical shifts to determine the geometric sensitivity of the latter. We could show that with the help of multivariate regression this approach leads to a semi-quantitative understanding of the colour tuning mechanism, without the need to have the explicit molecular conformations of the chromophore.

A series of molecular conjugation scenarios using mesomeric formulae that explicitly describe the possibilities of bond alternation in PCB helps to connect the statistical predictions with the actual conjugation patterns found in P and P state. An important role is played by the stabilisation of the positive charge in tertiary carbocations or by ionic interaction with the protein.

Showcasinpan>g the examples of cyanobacteriochrome AnPixJg2 and phytochrome Cph12, our approach represents an alternative method to interpret NMR data that can readily be transferred to related tetrapyrrole chromophores in various protein environments.

5. Materials and Methods

5.1. Preparations

Starting geometries for P states of AnPixJg2 and Cph12 were obtainpan>ed as X-ray structures from the Protein Data Bank (PDB; AnPixJg2: 3W2Z, [13] Cph12: 2VEA [35]). The simulation cell was set up using VMD, [74] solvating the protein in a box of TIP3P water [75]. Possible protonation states of protein residues were set to correspond to pH 7.0. The protonation state of histidine was chosen as follows: AnPixJg2: His293-E, His322-P, rest: D; Cph12: His260-E, His290-D, rest: D. The PCB chromophore was chosen with fully protonated pyrrole rings according to the experimental evidence [68].

5.2. Force Field Molecular Dynamics Simulations

FFMD simulations were carried out with the NAMD package [76] under periodic boundary conditions (PBC) using the CHARMM22 force field [77,78] in the isothermal-isobaric (NPT) ensemble and the combined Nosé–Hoover Langevin piston method with a period of 200 fs [79,80]. For PCB, force field parameters obtained by Mroginski and co-workers were used [58]. For van der Waals interactions a cutoff of 10 Å was employed. The initial simulation cell dimensions were chosen to include a water layer of 30 Å thickness. Bond lengths between heavy atoms and hydrogen atoms were held fixed using the SHAKE algorithm [81] with a time step of 2 fs. The system was equilibrated, first, by optimisation of the water cell, followed by heat up and equilibration keeping the protein positions fixed. Then, the protein structure was optimised with fixed water positions. Final equilibration for 200 ps was carried out after heating up the entire simulation cell to the desired temperature within 200 ps, followed by the production run (see Table 3 for more details).

Table 3

Computational details of force field and ab initio molecular dynamics calculations for the protein study.

Parameter	FFMD	AIMD	NMR
Ensemble	NPT	NVT	NVT
No. of atoms	37,651/104,960 *	≈100_QM	≈280_QM
Temperature	300 K	300 K	–
Pressure	1.01325 bar	–	–
Barostat	Langevin	–	–
Thermostat	Nosé–Hoover	Nosé–Hoover	–
ABC_QM (Å)	– **	30.0, 30.0, 30.0	35.0, 35.0, 35.0
Force Field/Functional	CHARMM22	BLYP-D3	BLYP-D3
Basis Set	–	TZVP-GTH	pcS2/3
Density Cutoff	–	320 Ry	400 Ry
No. of Snapshots	–	–	405
t_NVP/t_NVT	≈1 μs	>20 ps	–
Other	SHAKE [81]

* AnPixJg2/Cph1Δ2 + solvent box (water); ** solvation layer thickness: 30 Å.

5.3. QM/MM Molecular Dynamics Simulations

Snapshots from FFMD simulations were transferred into a QM/MM setup. AIMD simulations were performed as the Born-Oppenheimer type in the NVT enpan>semble usinpan>g CP2K 5.1 [82]. The partitioning of the protein into QM and MM part, as well as employed capping atoms, is listed in Table 4. QM/MM bond interfaces were handled using an optimised capping potential introduced between for the bond between and carbon atoms [83,84]. The MM part was described by the CHARMM22 force field (see previous section). The size of the QM box was set to 30.0, 30.0, 30.0 Å. The BLYP functional [85,86] was employed together with Grimme’s dispersion correction (D3), ref. [87] using the GPW scheme [88], GTH pseudo-potentials, [89,90,91] a density cutoff of 320 Ry, and the TZVP-GTH basis set. The QM/MM interaction term was evaluated within the GEEP scheme presented by Laino and co-workers [92,93]. Equilibration runs of length 5–10 ps were carried out under massive thermostatting using Nosé–Hoover chains of length 5 at 330 K and a coupling constant of 10 fs [80,94]. For subsequent production runs, QM and MM subsystems were coupled to separate thermostat units with a coupling constant of 500 fs.

Table 4

Information on protein residues included into the QM part, and corresponding capping atoms for response calculations.

AnPixJg2				Cph1Δ2
Residue	AIMD	NMR	Capping Atom	Residue	AIMD	NMR	Capping Atom
PCB	x	x	–	PCB	x	x	–
Trp289	x	x	CA	Asp207	x	x	CB/CA **
Asp291	x	x	CA	Ile208	x	x	C/C,CB **
His293	–	x	CA	Arg222	–	x	CB
Arg301	–	x	CA	Arg254	–	x	CB
Tyr302	–	x	CA	Ala256	–	x	CA
Phe308	–	x	CA	Tyr257	–	x	–
His318	–	x	CA	His258	x	x	CA/CB **
Phe319	–	x	CA	Cys259	x	x	–
Ser320	x	x	CA	His260	x	x	C
Cys321	x	x	–	Tyr263	–	x	CA
His322	x	x	C	His290	–	x	CA
Tyr334	–	x	CA	Water *	x	x	–
Water *	–	x	–

* Water molecules within a distance of 2.6 Å to PCB. ** AIMD /NMR.

5.4. Calculation of Nuclear Shieldings

NMR calculations were performed in the same setup as AIMD simulations, but the size of the QM/MM part was chosen adaptively such that RMSE of chemical shifts were below 1% (cf. Figure A1). Water molecules within a range of 2.6 from the chromophore were included in the QM part as well. Snapshots were extracted from the AIMD trajectory every 200 fs. All atoms in the QM region were treated at an all-electron level using the Gaussian-augmented plane waves (GAPW) method of CP2K 5.1 [95] with density cutoff of 400 Ry and pcS-3 basis set (pcS-2 on oxygen and sulphur atoms) [96]. The size of the QM box was set to 35.0, 35.0, 35.0 Å. Gauge origin was treated using the IGAIM method of CP2K, [95] which is an adaption of the CSGT method based on Atoms in Molecules [97].

5.5. Statistical Analysis

Principal component analysis (PCA) and multivariate regression was carried out with R using the packages tidyverse, ade4, and GGally [98]. For PCA the centred, unscaled C chemical shifts of PCB atoms C–C were used as variables to account for the importance of fluctuations in these features. PCA was not used for dimensionality reduction; though, consequently, the number of components equalled the number of input variables (6; cf. Figure 4).

Figure 4

Principal component analysis (PCA) of C chemical shift data obtained from the QM/MM MD simulations. From left to right: eigenvalues, components with PC 1 as abscissa and PC 2, 3, 4 as ordinate, respectively. PC 1–4 contain a cumulative projected inertia of 79%.

Multivariate regression and prediction was carried out, inverting the correlation matrix of C shifts with C-C bond lengths in rings C and D, shown in Figure 3, based on the relationship

being the change inpan> C chemical shift (features), the change in bond length (outcome/prediction), and denoting elements of the linear correlation matrix (Figure 3); using the previously obtained principal components.

6. Other

Python-based MDAnalysis 1.0.0 [99,100] was used for handling protein topologies. Geometric analysis, QM/MM setups, as well as genpan>eral post-processing of data was carried out using our own implementations in Python 3.8 [101,102,103,104] and numerical libraries NumPy 1.19.1 and SciPy 1.5.2 [105,106]. Plots were generated with Matplotlib 3.3.0 [107]. Molecular visualisations were created with VMD [74] using the Tachyon Ray Tracer [108] All structural formulae were created with using chemfig 1.56 [109].