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

class="Chemical">Tetrapyrroles showcase aclass="Chemical">n overwhelmiclass="Chemical">ng abuclass="Chemical">ndaclass="Chemical">nce iclass="Chemical">n liviclass="Chemical">ng species, accompaclass="Chemical">nied by aclass="Chemical">n impressive fuclass="Chemical">nctioclass="Chemical">nal versatility [1]. They caclass="Chemical">n be classified iclass="Chemical">nto cyclic aclass="Chemical">nd opeclass="Chemical">nclass="Chemical">n class="Chemical">-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 class="Chemical">bilin cofactor, phycocyaclass="Chemical">noclass="Chemical">n class="Chemical">bilin (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 tuning in 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 lenclass="Chemical">gth [44]. Receclass="Chemical">ntly, Wiebeler aclass="Chemical">nd co-workers showed the hypsochromic shift iclass="Chemical">n the greeclass="Chemical">n-absorbiclass="Chemical">ng photoproduct of Slr1393g3 (a PCBclass="Chemical">n class="Chemical">-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 nonclass="Chemical">-covaleclass="Chemical">nt iclass="Chemical">nteractioclass="Chemical">ns without losiclass="Chemical">ng accuracy at describiclass="Chemical">ng the chromophore aclass="Chemical">nd its iclass="Chemical">nstaclass="Chemical">ntaclass="Chemical">neous thermal fluctuatioclass="Chemical">ns. We focus oclass="Chemical">n two cyaclass="Chemical">nobacterial photoseclass="Chemical">nsor proteiclass="Chemical">ns (Figure 1): the red/far-red switchiclass="Chemical">ng class="Chemical">n class="Chemical">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 inclusion 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 lenclass="Chemical">gth aclass="Chemical">nd the C chemical shift of class="Chemical">n class="Disease">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 uncertainty in the obtained averages of C chemical shifts suggests a high sensitivity of the class="Chemical">carbon atoms’ shieldiclass="Chemical">ng towards aclass="Chemical">ny direct or iclass="Chemical">ndirect chaclass="Chemical">nge iclass="Chemical">n the eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">nt, the further aclass="Chemical">nalysis of their depeclass="Chemical">ndeclass="Chemical">nce oclass="Chemical">n chromophore–apoproteiclass="Chemical">n coclass="Chemical">nclass="Chemical">nectioclass="Chemical">ns aclass="Chemical">nd temperature effects eclass="Chemical">nables a more reficlass="Chemical">ned class="Chemical">notioclass="Chemical">n of the matter. Computatioclass="Chemical">nally, it is straightforward to repeat chemical shift calculatioclass="Chemical">ns for the same setup, gradually “switchiclass="Chemical">ng-oclass="Chemical">n” those proteiclass="Chemical">n residues that iclass="Chemical">nteract with the chromophore, or class="Chemical">numerically “cooliclass="Chemical">ng-dowclass="Chemical">n” the system, which forces the class="Chemical">n class="Disease">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 n class="Disease">carbon atoms C–C, though beiclass="Chemical">ng bliclass="Chemical">nd towards supramolecular effects, are aclass="Chemical">n importaclass="Chemical">nt probe of their immediate covaleclass="Chemical">nt situatioclass="Chemical">n.

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

The interconnection of nuclear shieldings and bond lenclass="Chemical">gths is class="Chemical">not particularly surprisiclass="Chemical">ng, but the remarkable “isolatioclass="Chemical">n” of C chemical shifts from supramolecular iclass="Chemical">nteractioclass="Chemical">ns deserves a closer look. Figure 3 preseclass="Chemical">nts the liclass="Chemical">near correlatioclass="Chemical">n matrix of C shifts with Cclass="Chemical">n class="Chemical">-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 sensitivity follows the clear rationale that an atom’s chemical shift increases with the lenclass="Chemical">gth of the origiclass="Chemical">naticlass="Chemical">ng boclass="Chemical">nds—aclass="Chemical">nd decreases with the class="Chemical">next, aclass="Chemical">nd iclass="Chemical">ncreases agaiclass="Chemical">n with the followiclass="Chemical">ng boclass="Chemical">nd aclass="Chemical">nd so oclass="Chemical">n. Iclass="Chemical">n Figure 3 this alterclass="Chemical">natioclass="Chemical">n is marked by red aclass="Chemical">nd blue bullets, respectively, shapiclass="Chemical">ng aclass="Chemical">n imagiclass="Chemical">ned double diagoclass="Chemical">nal; like Figure 2b, it is aclass="Chemical">nother represeclass="Chemical">ntatioclass="Chemical">n of the uclass="Chemical">nderlyiclass="Chemical">ng coclass="Chemical">njugated system. The symmetry of the geometric seclass="Chemical">nsitivity, however, is clearly disturbed leaviclass="Chemical">ng certaiclass="Chemical">n class="Chemical">nuclei with a “preferred” boclass="Chemical">nd they are correlated with. For iclass="Chemical">nstaclass="Chemical">nce, there is a stroclass="Chemical">ng coclass="Chemical">nclass="Chemical">nectioclass="Chemical">n of class="Chemical">n class="Disease">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 lenn class="Chemical">gth iclass="Chemical">n P aclass="Chemical">nd P states of PCB alike, assumiclass="Chemical">ng that the correlatioclass="Chemical">n matrix iclass="Chemical">n Figure 3 remaiclass="Chemical">ns uclass="Chemical">naltered iclass="Chemical">n either state.

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

In 2015, Song et al. reported a peculiar difference pattern of C nuclear shieldings () in PCB upon photoconversion in AnPixJg2 and class="Chemical">Cph12: Nuclei beloclass="Chemical">ngiclass="Chemical">ng to riclass="Chemical">ngs C aclass="Chemical">nd D exhibit chaclass="Chemical">nges iclass="Chemical">n chemical shifts of up to 10 ppm wheclass="Chemical">n traclass="Chemical">nsformiclass="Chemical">ng from the P state iclass="Chemical">nto the respective photoproduct [68]. The authors attribute this observatioclass="Chemical">n to a “large-scale rearraclass="Chemical">ngemeclass="Chemical">nt of the critical iclass="Chemical">nteractioclass="Chemical">ns [that] caclass="Chemical">nclass="Chemical">not be explaiclass="Chemical">ned oclass="Chemical">nly by the coclass="Chemical">nfiguratioclass="Chemical">nal chaclass="Chemical">nges of the chromophore.” However, recalliclass="Chemical">ng that C chemical shifts are bliclass="Chemical">nd towards supramolecular perturbatioclass="Chemical">ns (cf. Figure A1) the chaclass="Chemical">nges caclass="Chemical">nclass="Chemical">not be explaiclass="Chemical">ned by the removal of class="Chemical">n class="Chemical">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 lenn class="Chemical">gths of PCB iclass="Chemical">n the P state, calculated from the MD trajectories. They support the iclass="Chemical">nferred chaclass="Chemical">nge iclass="Chemical">n boclass="Chemical">nd order, for iclass="Chemical">nstaclass="Chemical">nce to what exteclass="Chemical">nt a siclass="Chemical">ngle or double boclass="Chemical">nd character is retaiclass="Chemical">ned. The giveclass="Chemical">n uclass="Chemical">ncertaiclass="Chemical">nties of predictioclass="Chemical">n help to disticlass="Chemical">nguish importaclass="Chemical">nt from less sigclass="Chemical">nificaclass="Chemical">nt treclass="Chemical">nds upoclass="Chemical">n photococlass="Chemical">nversioclass="Chemical">n.

Comparing the two P states of AnPixJg2 and class="Chemical">Cph12, we do class="Chemical">not ficlass="Chemical">nd big differeclass="Chemical">nces betweeclass="Chemical">n the two proteiclass="Chemical">ns, which correspoclass="Chemical">nds to the earlier argumeclass="Chemical">nt that PCB iclass="Chemical">n both proteiclass="Chemical">ns is chemically equivaleclass="Chemical">nt. The fouclass="Chemical">nd class="Chemical">numbers are maiclass="Chemical">nly due to the comparably large differeclass="Chemical">nce iclass="Chemical">n chemical shift for atom C, lyiclass="Chemical">ng at the joiclass="Chemical">nt of riclass="Chemical">ngs C aclass="Chemical">nd D aclass="Chemical">nd which is coclass="Chemical">nclass="Chemical">nected to maclass="Chemical">ny boclass="Chemical">nds (cf. Figure 3). Goiclass="Chemical">ng over to predicted chaclass="Chemical">nges iclass="Chemical">n boclass="Chemical">nd leclass="Chemical">nclass="Chemical">n class="Chemical">gth 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 sensitivity of C chemical shifts are based on experimental evidence reported by Song et al. [68] We inferred expected changes in bond lenclass="Chemical">gth of the PCB chromophore wheclass="Chemical">n goiclass="Chemical">ng from the dark state P to the photostates P for Aclass="Chemical">nPixJg2 aclass="Chemical">nd P for class="Chemical">n class="Chemical">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 lenclass="Chemical">gth. Oclass="Chemical">ne has to keep iclass="Chemical">n miclass="Chemical">nd that the alterclass="Chemical">natioclass="Chemical">n patterclass="Chemical">n represeclass="Chemical">nts the predomiclass="Chemical">naclass="Chemical">nt character of the Cclass="Chemical">n class="Chemical">-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

Considering that the shaping of delocalisation of double bonds by means of mesomeric scn class="Chemical">hemes has so far class="Chemical">not received due atteclass="Chemical">ntioclass="Chemical">n, we argue that workiclass="Chemical">ng out what is actually possible aclass="Chemical">nd favourable iclass="Chemical">n terms of boclass="Chemical">nd alterclass="Chemical">natioclass="Chemical">n is a key to fully uclass="Chemical">nderstaclass="Chemical">nd the colour tuclass="Chemical">niclass="Chemical">ng pheclass="Chemical">nomeclass="Chemical">noclass="Chemical">n iclass="Chemical">n 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 lenclass="Chemical">gth. We coclass="Chemical">nsider the coclass="Chemical">njugatioclass="Chemical">n betweeclass="Chemical">n riclass="Chemical">ngs B aclass="Chemical">nd C, which is maiclass="Chemical">ntaiclass="Chemical">ned iclass="Chemical">n all states of the chromophore, as core coclass="Chemical">njugatioclass="Chemical">n. Assumiclass="Chemical">ng that iclass="Chemical">n P state riclass="Chemical">ng D is twisted out of placlass="Chemical">ne, core coclass="Chemical">njugatioclass="Chemical">n aloclass="Chemical">ne carries the absorptioclass="Chemical">n of greeclass="Chemical">n light, which is iclass="Chemical">n liclass="Chemical">ne with what has beeclass="Chemical">n reported by Wiebeler et al. [17] There exist several mesomeric formulae describiclass="Chemical">ng this coclass="Chemical">njugatioclass="Chemical">n, the most importaclass="Chemical">nt of which are showclass="Chemical">n iclass="Chemical">n Scheme 2. Thereiclass="Chemical">n, the positive charge is delocalised aloclass="Chemical">ng with the chaclass="Chemical">nge of boclass="Chemical">nd alterclass="Chemical">natioclass="Chemical">n, but the high electroclass="Chemical">negativity of the class="Chemical">n class="Chemical">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 class="Chemical">amide resoclass="Chemical">naclass="Chemical">nce pulliclass="Chemical">ng electroclass="Chemical">n declass="Chemical">nsity from the class="Chemical">n class="Chemical">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 conclude that a plausible prediction of the geometric bond lenclass="Chemical">gth coclass="Chemical">njugatioclass="Chemical">n patterclass="Chemical">n of the PCB iclass="Chemical">n its P state is possible through measured aclass="Chemical">nd computed NMR chemical shifts (aclass="Chemical">nd their geometric seclass="Chemical">nsitivities). The resulticlass="Chemical">ng coclass="Chemical">njugatioclass="Chemical">n patterclass="Chemical">n, iclass="Chemical">n turclass="Chemical">n, provides a coclass="Chemical">nsisteclass="Chemical">nt explaclass="Chemical">natioclass="Chemical">n of the chaclass="Chemical">nge iclass="Chemical">n absorptioclass="Chemical">n waveleclass="Chemical">nclass="Chemical">n class="Chemical">gth 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 wavelenn class="Chemical">gth is due to the exteclass="Chemical">nsioclass="Chemical">n of coclass="Chemical">njugatioclass="Chemical">n, which class="Chemical">now iclass="Chemical">ncludes riclass="Chemical">ngs B, C, aclass="Chemical">nd D siclass="Chemical">nce the latter rotates iclass="Chemical">n the placlass="Chemical">ne with the core [44]. We shall call this full coclass="Chemical">njugatioclass="Chemical">n because it resembles the loclass="Chemical">ngest possible resoclass="Chemical">naclass="Chemical">nce iclass="Chemical">n PCB. Iclass="Chemical">ncludiclass="Chemical">ng riclass="Chemical">ng D iclass="Chemical">n the exchaclass="Chemical">nge of -electroclass="Chemical">ns has class="Chemical">not oclass="Chemical">nly boclass="Chemical">nds C–C aclass="Chemical">nd C–C toggle their siclass="Chemical">ngle–double boclass="Chemical">nd patterclass="Chemical">n, but also withdraws charge declass="Chemical">nsity from the riclass="Chemical">ng as the positive charge is class="Chemical">now shared as well.

Scheme 4 presents 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 class="Chemical">carbon atoms do class="Chemical">not directly “see” supramolecular effects; however, they are seclass="Chemical">nsitive to iclass="Chemical">ndirect effects, which iclass="Chemical">n the preseclass="Chemical">nt case is the stabilisatioclass="Chemical">n of the positive charge by the apoproteiclass="Chemical">n: Rohmer et al. argued that iclass="Chemical">n P the class="Chemical">negatively charged class="Chemical">n class="Chemical">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 in its P state a plausible prediction of the bond lenn class="Chemical">gth coclass="Chemical">njugatioclass="Chemical">n patterclass="Chemical">n is obtaiclass="Chemical">ned, based solely oclass="Chemical">n measured aclass="Chemical">nd computed NMR chemical shifts (aclass="Chemical">nd their geometric seclass="Chemical">nsitivities). Agaiclass="Chemical">n, the resulticlass="Chemical">ng coclass="Chemical">njugatioclass="Chemical">n patterclass="Chemical">n serves as explaclass="Chemical">natioclass="Chemical">n of the chaclass="Chemical">nge iclass="Chemical">n absorptioclass="Chemical">n waveleclass="Chemical">nclass="Chemical">n class="Chemical">gth through photoconversion.

4. Conclusions

We have addressed the colour tuning mechanism in the photosensor proteins AnPixJg2 and class="Chemical">Cph12, which are based oclass="Chemical">n the phycocyaclass="Chemical">noclass="Chemical">n class="Chemical">bilin (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 complementary 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 lenclass="Chemical">gth patterclass="Chemical">n. This caclass="Chemical">n be realised iclass="Chemical">n the form of a matrix represeclass="Chemical">nticlass="Chemical">ng the traclass="Chemical">nsformatioclass="Chemical">n of the boclass="Chemical">nd leclass="Chemical">nclass="Chemical">n class="Chemical">gth 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.

Showcasing the examples of cyanobacteriochrome AnPixJg2 and phytochrome n class="Chemical">Cph12, our approach represeclass="Chemical">nts aclass="Chemical">n alterclass="Chemical">native method to iclass="Chemical">nterpret NMR data that caclass="Chemical">n readily be traclass="Chemical">nsferred to related class="Chemical">n class="Chemical">tetrapyrrole chromophores in various protein environments.

5. Materials and Methods

5.1. Preparations

Starting geometries for P states of AnPixJg2 and class="Chemical">Cph12 were obtaiclass="Chemical">ned as X-ray structures from the Proteiclass="Chemical">n Data Baclass="Chemical">nk (PDB; Aclass="Chemical">nPixJg2: 3W2Z, [13] class="Chemical">n class="Chemical">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 class="Chemical">water layer of 30 Å thickclass="Chemical">ness. Boclass="Chemical">nd leclass="Chemical">nclass="Chemical">n class="Chemical">gths 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 ensemble using 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 class="Chemical">carbon atoms [83,84]. The MM part was described by the CHARMM22 force field (see previous sectioclass="Chemical">n). The size of the QM box was set to 30.0, 30.0, 30.0 Å. The BLYP fuclass="Chemical">nctioclass="Chemical">nal [85,86] was employed together with Grimme’s dispersioclass="Chemical">n correctioclass="Chemical">n (D3), ref. [87] usiclass="Chemical">ng the GPW scheme [88], class="Chemical">n class="Chemical">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). class="Chemical">Water molecules withiclass="Chemical">n a raclass="Chemical">nge of 2.6 from the chromophore were iclass="Chemical">ncluded iclass="Chemical">n the QM part as well. Sclass="Chemical">napshots were extracted from the AIMD trajectory every 200 fs. All atoms iclass="Chemical">n the QM regioclass="Chemical">n were treated at aclass="Chemical">n all-electroclass="Chemical">n level usiclass="Chemical">ng the Gaussiaclass="Chemical">n-augmeclass="Chemical">nted placlass="Chemical">ne waves (GAPW) method of CP2K 5.1 [95] with declass="Chemical">nsity cutoff of 400 Ry aclass="Chemical">nd pcS-3 basis set (pcS-2 oclass="Chemical">n class="Chemical">n class="Chemical">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 n class="Disease">atoms C–C were used as variables to accouclass="Chemical">nt for the importaclass="Chemical">nce of fluctuatioclass="Chemical">ns iclass="Chemical">n these features. PCA was class="Chemical">not used for dimeclass="Chemical">nsioclass="Chemical">nality reductioclass="Chemical">n; though, coclass="Chemical">nsequeclass="Chemical">ntly, the class="Chemical">number of compoclass="Chemical">neclass="Chemical">nts equalled the class="Chemical">number of iclass="Chemical">nput 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 Cclass="Chemical">-C boclass="Chemical">nd leclass="Chemical">nclass="Chemical">n class="Chemical">gths in rings C and D, shown in Figure 3, based on the relationship

being the change in C chemical shift (features), the change in bond lenn class="Chemical">gth (outcome/predictioclass="Chemical">n), aclass="Chemical">nd declass="Chemical">noticlass="Chemical">ng elemeclass="Chemical">nts of the liclass="Chemical">near correlatioclass="Chemical">n matrix (Figure 3); usiclass="Chemical">ng the previously obtaiclass="Chemical">ned priclass="Chemical">ncipal compoclass="Chemical">neclass="Chemical">nts.

6. Other

Python-based MDAnalysis 1.0.0 [99,100] was used for handling protein topologies. Geometric analysis, QM/MM setups, as well as general 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].