Quantifying electron transfer reactions in biological systems: what interactions play the major role?

Various biological processes involve the conversion of energy into forms that are usable for chemical transformations and are quantum mechanical in nature. Such processes involve light absorption, excited electronic states formation, excitation energy transfer, electrons and protons tunnelling which for example occur in photosynthesis, cellular respiration, DNA repair, and possibly magnetic field sensing. Quantum biology uses computation to model biological interactions in light of quantum mechanical effects and has primarily developed over the past decade as a result of convergence between quantum physics and biology. In this paper we consider electron transfer in biological processes, from a theoretical view-point; namely in terms of quantum mechanical and semi-classical models. We systematically characterize the interactions between the moving electron and its biological environment to deduce the driving force for the electron transfer reaction and to establish those interactions that play the major role in propelling the electron. The suggested approach is seen as a general recipe to treat electron transfer events in biological systems computationally, and we utilize it to describe specifically the electron transfer reactions in Arabidopsis thaliana cryptochrome–a signaling photoreceptor protein that became attractive recently due to its possible function as a biological magnetoreceptor.


Simulation protocol
The hybrid quantum mechanics/molecular mechanic (QM/MM) calculations with the PE approach are rather involved, as they include many technical steps, necessary to perform the calculations accurately. Below we summarize the crucial steps that are needed for such calculations. Each calculation constitutes four steps: 1. Fragmenting the molecular system into the core and the environment regions.
2. Obtaining the embedding potential for the environment region.
3. Handling the molecular interface between the core and the environment regions. 4. Calculating the desired properties of the core region.
First the protein is subdivided into two regions, the core region, or the active site, and the environment (here we denote the molecular subsystem outside the active site as the environment), by using the dedicated VMD 1 extension pe-fragmenter (available in house).
During this procedure, covalent bonds connecting the active site and the environment are broken, and the resulting dangling bonds are then capped with hydrogen atoms. In the present study, the active site consists of the side chains of the three tryptophans and the isoalloxazine moiety of the FAD, as shown in Fig. 2A. The environment includes the remainder of the protein and a shell of water molecules and ions surrounding the protein, with a layer of 12Å, see Fig. 2B Second, the embedding potential is calculated using an in-house script, developed by one of the authors. Here all residues (amino acids) of the environment are treated as overlapping fragments, according the MFCC approach, and the potential is represented through electrostatic components of the individual residues. For AtCry the protein is divided into S2 Figure S1: Computational description of the interface between the core and environment regions. Graphical depiction of an overlap between the active site and the environment, due to breaking of covalent bonds in the course of the PE calculations and termination of the dangling bonds with hydrogens. A) between W A backbone and its sidechain, B) between the isoalloxazine moiety and FAD remains. In both cases, the charges of the atoms marked purple are redistributed to the atoms marked green. 969 fragments. Next, the interface between the core and environment subsystems should be taken care of. Capping the cut covalent bonds between the active site and the environment leads to an near-overlap of atoms from the two regions, as for example illustrated in Fig. S1. This overlap is handled by redistributing the charges from the atoms that clash with the active site, to the two nearest heavy atoms, equally, to the two nearest heavy atoms, i.e. non-hydrogen atoms. Figure S1 shows examples of such redistributions in the case of tryptophan and isoalloxazine cofactor. The remaining contributions to the embedding potential, i.e. higher-order multipole moments and polarizabilities for the conflicting atoms are deleted. Finally, the properties of the active site in the embedding potential are calculated using quantum chemistry software, e.g. Dalton 2,3 .

Environment stabilizes radical pair states.
The electronic excitation spectra of the AtCry active site are studied for the four optimized structural configurations of the active site, namely,  Here the interactions of the active site with the molecular environment have been decomposed into five contributions, representing the environment through: partial charges (q 0 ), dipole moments (d) and quadrupole moments (Q) ascribed to each atom of the protein; point charges (q ESP ) fitted to reproduce the electrostatic potential; induced dipole moments on each atom of the environment, calculated from the ground state of the active site (α 0 ); induced dipole moments of all atoms of the environment, that take into account also the charge redistribution in the system upon electronic excitation (α 1 ).

S5
is more stabilized in the case of the RP-B opt structure, while RP-C stabilization leads to a dramatic energy decrease in the case of RP-C opt structure. This stabilization is absent in the vacuum model, however, it is of significant importance since it indicates that environment leads to the energetic favoring of a certain electronic state.
Electrostatic interactions and polarization capture the physics of electron transfer.

Electronic properties of AtCry could be modeled with CAM-B3LYP.
The calculations for the environment model are carried out with TDDFT method using the CAM-B3LYP exchange-correlation functional, as described in Methods. To examine whether TDDFT can adequately describe the electron transfers in AtCry we have considered the electronic excitation spectra of its active site, optimized in four structural configurations in S6 vacuum (the CS opt state, the radical pair states RP-A opt , RP-B opt and RP-C opt , see Fig. 2A and Fig. 3). For the optimized configurations we then have computed additionally the electronic excitation spectra by using the XMCQDPT-2 9 , CAM-B3LYP 10-13 and CASSCF 14,15 methods. Figure S2 (left part) summarizes the comparison between the different calculations.
It can be seen that for all the structural configurations of the active site, the order of electronic excitations is the same for the three methods, indicating that overall CAM-B3LYP yields qualitatively correct results. For all structural configurations, except the RP-A opt , the energy of all radical pair states computed with CAM-B3LYP is lowered, as compared to the energies calculated with XMCQDPT-2 and CASSCF methods. It should be noted that CASSCF is known to overestimate the energy of excited states 16 , which accounts for some of the discrepancies between the CASSCF and CAM-B3LYP approaches. S7