Computational Investigations on Phycocyanobilin ^†

Abstract

Phycocyanobilin was computationally investigated by means of DFT calculations in combination with implicit solvation starting from X-ray data. Different conformations and degrees of protonation were considered, and the acidity constants were estimated. The computed data suggest a syn-syn-syn conformation for the molecule, with the two carboxylic groups deprotonated under physiological conditions and weak acidic behavior of one of the pyrrolone heterocycles. The absorption transitions in the visible range were studied by means of TD-DFT calculations, focusing on the molecular orbitals involved. The frontier orbitals have a dominant role in the lowest energy absorption.

1. Introduction

Phycocyanobilin (PCB) is an antenna pigment found in cyanobacteria such as Arthrospira Platensis or Synechococcus Elongatus and in some algae. From a structural point of view, PCB is a linear tetrapyrrole with two carboxyl groups and some hydrocarbon substituents, the terminal heterocycles actually being a pyrrolidinone and a pyrrolone, given the presence of a carbonyl function in the alpha position with respect to the nitrogen atoms (Figure 1) [1,2,3,4,5,6]. PCB is a member of the biline family, which are pigments such as phycoerythrobilin or biliverdin [7,8]. In nature, PCB is covalently bound via thioether bonding with cysteine, thus resulting in the prosthetic group of some proteins such as phycocyanin-C (C-PC), allophycocyanin, or phycoerythrin. The interaction with sulfur involves the exocyclic double bond present in one of the terminal pyrrolones of the free PCB [9]. In general, phycobilins are bound to proteins with masses between 30 and 35 kDa. These proteins are then organized into higher structures called phycobilisomes located in the lamellar membranes of cyanobacteria and algae [10].

The main biological role of pigments such as PCB is to assist the photosynthetic process by acting as antennae, i.e., collecting radiant energy from the sun and transferring it to the photosynthetic reaction centers. The commonly accepted mechanism is the Förster resonance energy transfer [11]. Such a role is evidenced by the lowest energy absorption band in the visible range of PCB-containing species, which is located in the red region of the spectrum and exhibits a pH-dependent nature. The luminescence of PCB in C-PC corresponds to an emission band centered at about 650 nm, with scarce Stokes shift with respect to the main absorption band [12,13]. Unlike porphyrins and related cyclic systems, the open tetrapyrrole structure of PCB has not been found to be associated with metal centers in biological derivatives subjected to structural studies. Investigations carried out on PCB obtained by a solvolitic approach showed a progressive quenching of luminescence by the addition of Hg(II) [14].

Given the intriguing features of PCB as a bio-based chromophore of commercial interest [15,16,17], the electronic structure of the molecule and selected related spectroscopic properties were investigated using computational methods based on the DFT and TDDFT theories and summarized here.

2. Computational Methods

The optimization of the hydrogen atoms added to the {C₃₃N₄O₆} skeleton derived from experimental data was carried out with the Merck Molecular Force Field (MMFF) [18], keeping the other atoms frozen. The software used was Spartan ‘16 version 2.0.3 [19]. The geometry optimizations were initially performed using the PBEh-3c method, a reparametrized version of the PBE0 hybrid functional (42% HF exchange) [20] that uses a split-valence double-zeta basis set (def2-mSVP) [21,22] and adds three corrections considering dispersion, basis set superposition, and other basis set incompleteness effects [23,24,25]. Further calculations were carried out with the r²SCAN-3c method [26], based on the meta-GGA r²SCAN functional [27] combined with a tailor-made triple-ζ Gaussian atomic orbital basis set. The method also includes refitted D4 and geometrical counter-poise corrections for London dispersion and basis set superposition error [23,28]. A third DFT method used was based on the TPSSh hybrid meta-GGA DFT functional [29] in combination with Ahlrichs’ def2-TZVP basis set [21]. The C-PCM implicit solvation model was added, considering water as a continuous medium [30]. IR simulations were carried out using the harmonic approximation, from which zero-point vibrational energies and thermal corrections (T = 298.15 K) were obtained. The electronic transitions were simulated by means of TDDFT calculations, considering 12 singlet and 12 triplet roots and including the spin-orbit correction [31]. Calculations were carried out using ORCA version 5.0.3 [32,33] and the output files were analyzed with Multiwfn version 3.8 [34].

3. Results and Discussion

In order to construct the molecular models of PCB, the deposited structure of C-PC in Synechococcus Elongatus (1.45 Å resolution) was considered as starting point [35]. The PCB chromophore was isolated from the surrounding protein chains and the hydrogen atoms were added coherently with the accepted Lewis structure. The positions of the hydrogen atoms were initially optimized by means of the MMFF method, keeping the other atoms frozen. The relative orientations of the four heterocycles in the starting structure determine an anti-syn-anti conformation (Figure 1).

All the DFT calculations were carried out in combination with the C-PCM implicit solvation model, considering water as continuous medium. The molecule was initially optimized in its neutral form, indicated as [PCB_ASA]⁰, by means of the PBEh-3c, r²SCAN-3c, and TPSSh/def2-TZVP methods. The geometry optimizations afforded comparable stationary points and maintained the initial anti-syn-anti conformation. Given the presence of four acidic protons, all the possible [PCB_ASA]¹⁻ monoanions derived from the formal deprotonation of [PCB_ASA]⁰ were optimized with the PBEh-3c method. The structures are sketched in Figure 2 with their relative Gibbs energy values. The most stable [PCB_ASA]¹⁻ isomer was optimized again using the TPSSh/def2-TZVP approach and its Gibbs energy value, together with that of [PCB_ASA]⁰, was used for the estimation of the pK_a1 value, calculated on the basis of a formal acid-base reaction between PCB and a chemically related species with known pK_a. Given the presence of the carboxylic groups in PCB, the acetic acid (HOAc)/acetate (OAc⁻) couple was considered as a proper reference (pK_a^HOAc = 4.756). The general reaction is thus [PCB_ASA]ⁿ⁻+ OAc⁻ → [PCB_ASA]⁽ⁿ⁺¹⁾⁻ + HOAc. The acidity constants of PCB were obtained according to equation 1, where ΔG represents the calculated Gibbs energy variation for the reaction, R is the gas constant, and T is the temperature (298 K). Starting from the most stable [PCB_ASA]¹⁻ isomer, the [PCB_ASA]²⁻ dianions were then optimized at the PBEh-3c level (Figure 2) and the isomer with the lowest Gibbs energy was used for the estimation of the pK_a2 value after further optimization with the TPSSh/def2-TZVP method. The same procedure was repeated for the third deprotonation, with the determination of the most stable [PCB_ASA]³⁻ isomer (Figure 2).

Table 1. Computed pK_a values for PCB_ASA and PCB_SSS. C-PCM/TPSSh/def2-TZVP calculations.

	PCBASA	PCBSSS
pKa1	3.1	3.3
pKa2	4.1	4.3
pKa3	10.5	6.9

collects the calculated pK_a values. According to the pK_a values, the two carboxyl groups of the anti-syn-anti conformer of PCB should be deprotonated at physiological pH.

pK_a^PCB = 4.756 + ΔG/(2.303RT)

(1)

The value of root mean square deviation (RMSD) between the TPSSh/def2-TZVP optimized {C₃₃N₄O₆} skeletons of [PCB_ASA]²⁻ and the corresponding X-ray data is 1.00 Å. r²SCAN-3c calculations afforded a strictly comparable stationary point, with a RMSD of 0.83 Å with respect to the experimental structure. The differences occurring between the computed structures and the starting geometry used are attributable to the removal of the apoprotein before the calculations and to the flexibility of some of the substituents on the heterocycles.

Relaxed surface scan calculations indicated that the energy barriers for the rotation of the external methenylpyrrolone units with respect to the central methenyldipyrrolic fragment are around 5 kcal mol⁻¹, marking other conformations of the free molecule accessible at room temperature. The plots of the relative energy values vs. the dihedral angles are shown in Figure 3. Moreover, the syn-syn-syn conformers of PCB, generally indicated with [PCB_SSS]ⁿ⁻ (n = 0, 1, 2, 3), resulted in being more thermodynamically stable than the corresponding [PCB_ASA]ⁿ⁻ conformers, as already stated in previous studies [36]. The Gibbs energy differences (ΔG) between [PCB_SSS]ⁿ⁻ and [PCB_ASA]ⁿ⁻ range between 6.4 and 11.6 kcal mol⁻¹, depending upon the charge of the species. The ΔG values are collected in

Table 2. Gibbs energy differences (kcal mol⁻¹) between the [PCB_SSS]ⁿ⁻ and [PCB_ASA]ⁿ⁻ conformers (n, 0, 1,2, 3). C-PCM/TPSSh/def2-TZVP calculations.

n	ΔG
0	−8.9
1	−6.4
2	−6.7
3	−11.6

.

Since the computational outcomes suggest the syn-syn-syn conformation as the most stable for the free molecule in solution, the estimation of the acidity constants was repeated. The pK_a values are reported in Table 1. The change in conformation strongly lowers the pK_a3 value, calculated at around 6.9 for the syn-syn-syn conformer. Free PCB should therefore be present under physiological conditions in equilibrium between its doubly and triply deprotonated forms.

The C-PCM/TPSSh/def2-TZVP optimized structures of [PCB_SSS]²⁻ and [PCB_SSS]³⁻ are shown in Figure 4. As for all the DFT-optimized structures considered here, the stationary points were characterized as local minima by means of IR simulations. On considering selected computed spectral features, the simulated bands with the highest wavenumbers, excluding the C-H and N-H stretchings, are associated with the pyrrolone carbonyl stretchings both for [PCB_SSS]²⁻ and [PCB_SSS]³⁻. The unscaled computed wavenumbers depend upon the charge of the molecule, being 1725 and 1662 cm⁻¹ for [PCB_SSS]²⁻ and 1650 and 1649 cm⁻¹ for [PCB_SSS]³⁻.

TDDFT calculations at the C-PCM/TPSSh/def2-TZVP level predicted the lowest energy S₁←S₀ absorptions at 610 nm (oscillator strength = 0.71) for [PCB_SSS]²⁻ and at 647 nm (oscillator strength = 0.66) for [PCB_SSS]³⁻, in line with the experimental interval [12,13]. In both cases, the transition has 100% LUMO←HOMO character. The red shift of the absorption moving from the dianion to the trianion agrees with the calculated 0.1 eV lowering of the HOMO-LUMO gap. The frontier molecular orbitals are shown in Figure 4. The transitions involve the π-delocalized electronic structure on the conjugated tetrapyrrolic structure, as clearly observable also from the plots of the hole and electron distributions in Figure 5 [37]. The PCB molecule can be considered aromatic, as indicated by the Shannon aromaticity index (SA) computed from the electron density values at the C-C, C-N and C-O bond critical points [38]. The SA values are equal to 5 × 10⁻³ both for [PCB_SSS]²⁻ and [PCB_SSS]³⁻. The alterations of the electronic structures caused by the S₁←S₀ transitions are better highlighted by plotting the differences between the hole and the electron distributions, an approach known as Charge Density Difference [37] (Figure 5). For instance, the lowest energy absorption makes the methine group bonded to the terminal pyrrolidinone electron poorer, while the opposite situation occurs for the other two bridging {CH} fragments.

It is worth noting that the lowest energy absorptions are predicted with a blue shift around 45 nm for the related [PCB_ASA]ⁿ⁻ conformers, highlighting that the conformation assumed by the tetrapyrrolic unit influences the spectral features of PCB-containing species [39]. On considering other degrees of protonation, the lowest energy S₁←S₀ absorption is calculated at 609 nm for both [PCB_SSS]⁰ and [PCB_SSS]¹⁻. The TDDFT calculations reported here suggest that the absorption features of PCB become pH-dependent when the acid-base equilibrium involves one of the N-heterocycles, while the protonation of the carboxylic groups negligibly influences the absorption maximum.

The addition of SOC to the calculations caused negligible effects to the predicted transitions, a result expected given the lack of heavy atoms in the structure. The lowest energy forbidden T←S₀ absorption was predicted for both [PCB_SSS]²⁻ and [PCB_SSS]³⁻ in the NIR range, at 1131 nm for the dianion and 1200 nm for the trianion. The luminescence quenching observed by adding a heavy metal ion such as Hg²⁺ [14] can be explained by admitting an increased rate of intersystem crossing toward triplet excited states, followed by non-radiative decay given the quite low energy gap. To better investigate this point, the triplet state geometries of [PCB_SSS]²⁻ and [PCB_SSS]³⁻ were optimized. The RMSD values with respect to the singlet ground state geometries are low, with values equal to 0.193 Å for [PCB_SSS]²⁻ and 0.075 Å for [PCB_SSS]³⁻. The hypothetic T→S₀ phosphorescent decays are predicted to have low energy, being the computed wavelengths equal to 1520 nm for the dianion and 1495 nm for the trianion. As for the S₁←S₀ absorptions, these transitions involve only the two frontier orbitals.

4. Conclusions

The computational outcomes provided here support the comprehension of the acid-base behavior of PCB, which influences the absorption features of the molecule. Another factor that affects the absorption maximum in the red range is the conformation assumed by the tetrapyrrolic fragment, which resulted in noticeable differences when comparing the free molecule and PCB-containing proteins. All these parameters should be considered when PCB is embedded in a matrix and applied as a bio-based pigment. Moreover, the use of PCB as a luminescent sensor for heavy metal ions can be rationalized by considering the low energy values computed for the triplet excited states, which favor non-radiative decay routes.