How orange carotenoid protein controls the excited state dynamics of canthaxanthin

Orange Carotenoid Protein (OCP) is a ketocarotenoid-binding protein essential for photoprotection in cyanobacteria. The main steps of the photoactivated conversion which converts OCP from its resting state to the active one have been extensively investigated. However, the initial photochemical event in the ketocarotenoid which triggers the large structural changes finally leading to the active state is still not understood. Here we employ QM/MM surface hopping nonadiabatic dynamics to investigate the excited-state decay of canthaxanthin in OCP, both in the ultrafast S2 to S1 internal conversion and the slower decay leading back to the ground state. For the former step we show the involvement of an additional excited state, which in the literature has been often named the SX state, and we characterize its nature. For the latter step, we reveal an excited state decay characterized by multiple timescales, which are related to the ground-state conformational heterogeneity of the ketocarotenoid. We assigned the slowly decaying population to the so-called S* state. Finally, we identify a minor decay pathway involving double-bond photoisomerization, which could be the initial trigger to photoactivation of OCP.

The calculations for the model were performed to have a comparison with reference data obtained with high level ab initio calculations, feasible for highly symmetric systems, in a previous work of the group [1].In this benchmark the excitation energies were computed at the GS minimum geometry in a Frank-Condon (FC) framework.The reference is reported next to the semiempirical calculations with different active spaces in Table S1A.For these calculations, the geometry of the model was optimized at the MP2 level since it can provide a sufficiently correct treatment of dynamic correlation, which is necessary for long conjugated chains to recover the correct degree of delocalization.The purpose of these calculations was to select a well suited active space to be used as a reference for the following calculations on CAN (for which it is not possible to have a reliable ab initio treatment).We therefore chose the orbital active space (14,12) as reference since it manages to treat at the same time in a sufficiently accurate way both the 2A − g and the 1B + u vertical energies, also including the n orbitals and allowing to compute the n → π * states (A u and B g ).
As done for the model system, we computed the excitation energies of CAN, both at S 0 and S 1 minimum geometries, this time optimized at the semiempirical level (Table S1B).We could verify

S1
Electronic Supplementary Material (ESI) for Chemical Science.This journal is © The Royal Society of Chemistry 2023 that the active space (14,12) previously chosen for the model and a smaller active space of (8,10) gave very similar results, which were also consistent with some experimental values from the literature for the vertical excitation energy of 1B + u [2], and the adiabatic energy of 2A − g [3].We therefore decided to adopt an active space of (8,10) in the simulations to further reduce the computational cost and because the n → π * states are not likely to be relevant in the excited state process.
Table S1: (A) Comparison between the Benchmark FC excitation energies (in eV) and the R-AM1/FOMO-CISD FC excitation energies (in eV) computed with different combinations of the AS and the w parameter.The GS minimum geometry is in every case obtained with MP2.Some of the ASs selected contain n type orbitals, therefore the states Au and Bg can be computed.(B) Excitation energies (in eV) for CAN's low-lying excited states computed at the R-AM1/FOMO-CISD level using the active spaces (14,12) and (8,10), and a Gaussian width (w) of 0.2 and 0.1 Hartree, respectively.On the left are listed the Frank-Condon excitation energies computed at the S0 minimum geometry, while on the right are displayed the excitation energies at S1 minimum geometry, and in both cases the state energies are reported relative to the S0 energy at S0 minimum geometry.All the geometries were optimized with the same R-AM1/FOMO-CISD method used for the calculation of the state energies.(C) The experimental data for the FC energy of 1B + u and the adiabatic energy of 2A − g are reported.
(A) model system data

S2 Thermal equilibrations in the ground state
An ensemble of nuclear coordinates and velocities of CAN was generated by running thermal equilibrations in the ground state (S 0 ).For the simulations of CAN in OCP, we performed ten QM/MM thermal equilibrations, propagated for 10 ps using the Bussi-Parrinello thermostat [4], with a relaxation time constant of 10 fs.The starting coordinates for the QM/MM thermal equilibrations were extrated from a classical molecular dynamics trajectory of OCP, previously computed in our group [5].For CAN in the gas phase, one thermal equilibration was propagated for 25 ps using the van Gunsteren-Berendsen thermostat [6], with a friction coefficient of 5.0 × 10 13 s −1 for all the atoms.In all the thermal equilibrations, we used a temperature of 300 K and a time step of 0.5 fs.In the sampling of initial conditions for the surface hopping simulations, we used the last 5 ps of the ten QM/MM thermal equilibrations of CAN in OCP, and the last 20 ps of the gas phase thermalization.
The data on the thermal equilibrations is reported below.

S3 Diabatization
A diabatic wavefunction [7,8,9,10] can be defined as a linear combination of adiabatic states and determine a unitary transformation from the adiabatic basis to the diabatic one: with |η⟩ is the vector of the diabatic basis functions and |φ⟩ is the one of the adiabatic ones.
T is a unitary rotation matrix built using some arbitrarily chosen reference states |R⟩ that must have a well defined character (ionic or neutral, localized or charge transfer excitation, etc.).By definition, a diabatic basis annihilates the derivative couplings, therefore it allows the preservation of the physical character and symmetry of each state near strong-coupling regions.
Instead of constructing a diabatic basis where the non-adiabatic couplings vanish completely (the so-called crude diabatic basis), it is possible to find a unitary matrix T transforming a truncated adiabatic basis into a strictly diabatic one only for one internal coordinate at a time or along a given path, and to obtain a basis where the couplings are not completely cancelled but they can become negligibly small.We use a diabatization technique where the reference states R chosen to build the rotation matrix T are the adiabatic states at S 0 minimum geometry: T is defined so that the sum of the overlaps between the diabatic wave functions and the reference is maximized, so that the the obtained diabatic functions are as similar as possible to the references.This is done by projecting the reference states onto the space spanned by the adiabatic states and then applying a Löwdin orthogonalization [11]: S is the overlap matrix between the projections of the reference states onto the adiabatic basis: The diabatic energies are the diagonal elements of the Hamiltonian matrix in the diabatic repre-sentation (which does not diagonalize the Hamiltonian, unlike the adiabatic basis): S4 Diabatic analysis of the ultrafast process To understand the nature of the excited states we adopt a quasi-diabatic representation (from now on, diabatic), defined on the basis of the overlap to reference states [7,10].Here the reference states of CAN were calculated at the S 0 minimum geometry.Figure S11 shows the diabatic states along the bond-length-alternation (BLA) coordinate with their adiabatic counterpart, painting a picture analogous to what found for Lutein [12].BLA is defined by the following equation: where N s and N d are respectively the number of single and double bonds and d s and d d are their respective lengths.At the Franck-Condon point, states S 1 -S 3 correspond to 2A − g , 1B + u , and 1B − u , respectively.Vertical excitation populates 1B + u (S 2 ), and then CAN rearranges its geometry towards smaller BLA values, following the excited-state gradient.The nuclear wavepacket oscillates back and forth along the BLA coordinate, encountering the crossing region between 1B + u and 1B − u (See Figure S12).Passing through the crossing, the nuclear wavepacket can either remain in the 1B + u state or switch to 1B − u .From an adiabatic point of view, in the former case the population of S 3 increases when the BLA reaches smaller values, while in the latter the population remains on S 2 .
Notably, the one-dimensional picture of Figure S11 shows a higher energy for the 1B − u minimum than 1B + u , but along the dynamics other degrees of freedom are activated, lowering the energy of the 1B − u state until its population becomes prevalent.After ∼100 fs, the S 2 state has a prevalent 1B − u character, and the population oscillations stop as 1B + u is not reachable anymore.This is confirmed    Figure S16: Distribution of the ∆E value (eV) between S1 and S0 at S1 → S0 hopping times of the surface hopping trajectories.On the left we show the simulations in OCP, the distribution in red corresponds to a sample of trajectories that were ran with no constraint on the energy gap for allowing the hops (i.e. the hops are allowed only if the energy difference between the involved states is ≤ than the imposed threshold), while in orange such threshold is of 1.0 eV.On the right we show the simulations in vacuum, where no constrained were applied.Table S2: Average values of several relevant properties for each group obtained from the puckering-based separation.From top to bottom, lifetime of the S1 excited state (fitting according to Equations S8 and S9 for p + and p − , respectively), similarity index computed for the dihedral DS in the initial conditions, similarity index computed for the dihedral DD at S1 → S0 hopping times, energy gap between S1 and S0 at hops and percentage of trajectories that undergo the hula-twist isomerization during the S1 → S0 decay.a The two different lifetimes of S 1 (∼3 ps and ∼6-7 ps) were ascribed to structural heterogeneity of the carotenoid in OCP O (at least two different carotenoid configurations are present) [16].

S7 Supporting figures
b Two different excitation wavelengths were used, namely 540 nm and 470 nm, giving rise to slightly different decay times (shorter for the 540-nm excitation) [16].c As in Ref. [22], the two extracted times were ascribed to the S 1 /ICT lifetimes of two different hECN populations in OCP, differing in hydrogen bonding via the carbonyl group [19].d These two different decay times were both attributed to a state in which S 1 is coupled with an intramolecular chargetransfer (ICT) state, indicated as S 1 /ICT [22].
e The ICT and S 1 names were assigned on the basis of the dominant contribution to the spectrum, although the two spectroscopic signatures contain both ICT and S 1 character [17,21].

Figure S1 :Figure S2 :Figure S3 :Figure S4 :Figure S5 :Figure S6 :Figure S7 :Figure S8 :
FigureS1: Data on the thermal equilibration trajectory of CAN in vacuum.From top to bottom: potential energy (Hartree), kinetic energy (Hartree), bond length alternation (BLA) coordinate value (angstrom), excitation energies of the first 5 low-lying excited states (Hartree).The dimmed section is the discarded non-equilibrium portion of the trajectory, the rest was used for the sampling of the SH initial conditions.

Figure S9 :Figure S10 :
Figure S9: Values of the terminal C-C dihedral angles (degrees) of CAN during the 10 thermal equilibrations of CAN in OCP.The numbers refer to the frames sampled from a classical MD trajectory [5].

Figure S11 :
Figure S11: Diabatic (whole lines) and adiabatic (dashed lines) potential energy curves for the first four low-lying states of CAN in vacuo.The plotted lines are second-degree polynomial function fittings of the data obtained from the relaxed scan along the BLA computed at the R-AM1/FOMO-CISD level.The arrows highlight the path followed by the electronic population.On the right there is a focused image of the avoided crossing between S2 and S3 (dashed lines); the solid lines are the diabatic states 1B + u and 1B − u .

Figure S12 :
Figure S12: Ultrafast portion of CAN excited-state dynamics in vacuo.Results are obtained in the same way discussed in the main text.The panels on the left show the adiabatic state populations of CAN during the first 200 fs of the SH simulations.Below are reported, respectively, the excitation energies from the ground state and the BLA values.On the right, the populations and the energies are shown in their diabatic representation, obtained through the diabatizaton discussed in Section S3.

Figure S13 :
Figure S13: Correlation between the Cremer and Pople's θ and ϕ coordinate in the case of CAN's β1 ring, a six-term ring with three sp2 carbons.

Figure S14 :
Figure S14: Two representative structures and distribution of the puckering coordinate (ϕ) of the β1 ring in the initial conditions of the SH trajectories.Structures with ϕ < 0 are associated with the p − conformation, while the ones with ϕ ≥ 0 define the p + conformation.

Figure S15 :
Figure S15: Absorption spectrum of CAN in vacuum (left) and in OCP (right).Every plotted value except the ones relative to S2 (the blue line) is multiplied by a factor of 5.

Figure S17 :
Figure S17: Distribution of the ∆E absolute value (eV) between S2 and S3.The purple distribution is referred to the ground state initial conditions of the surface hopping trajectories, the blue distribution is the energy gap between S2 and S3 at S2 → S3 and S3 → S2 hopping times.On the left we show the data from the simulations in OCP, on the right from the simulations in vacuum.

Figure S18 :
Figure S18: Correlation between the puckering conformation and the distributions of the distances of 5 CTD residues.

Figure S19 :
Figure S19: Left panel: scatter plot of the PCA observations highlighting the value of ϕ (puckering angle) of each geometry with the red/light-blue color gradient (p − in red and p + in light-blue).Right panel: the same plot of the PCA observations in which is reported with the purple/yellow color gradient the value of sDS (purple corresponds to more distorted DS and yellow to more flat DS).