Inhomogeneous ensembles of radical pairs in chemical compasses

The biophysical basis for the ability of animals to detect the geomagnetic field and to use it for finding directions remains a mystery of sensory biology. One much debated hypothesis suggests that an ensemble of specialized light-induced radical pair reactions can provide the primary signal for a magnetic compass sensor. The question arises what features of such a radical pair ensemble could be optimized by evolution so as to improve the detection of the direction of weak magnetic fields. Here, we focus on the overlooked aspect of the noise arising from inhomogeneity of copies of biomolecules in a realistic biological environment. Such inhomogeneity leads to variations of the radical pair parameters, thereby deteriorating the signal arising from an ensemble and providing a source of noise. We investigate the effect of variations in hyperfine interactions between different copies of simple radical pairs on the directional response of a compass system. We find that the choice of radical pair parameters greatly influences how strongly the directional response of an ensemble is affected by inhomogeneity.

landscapes 41 , proteins motions occur over many time scales, and one expects protein fluctuations generally being larger for slower time scales. Such motion on time scales larger than radical pair reaction times will result in varied conformational states even in an ensemble of initially completely identical protein complexes. The protein movement affects quantum states of cofactors, as observed as inhomogeneous broadening of energy lines, e.g. in photosynthetic proteins 42 . Fluctuations of hyperfine interactions in a flavin cofactor were estimated on short time scale 43 , and similarly sized variations are expected to occur on larger time scale as well. Here, we assume a relatively modest level of variation, i.e. 2% of the value of hyperfine interactions.
In addition to the unavoidable inhomogeneities due to molecular motion at physiological temperatures, additional variations of hyperfine couplings may occur if the magnetoreceptor ensemble involves inherently heterogeneous copies. For the candidate photo-magnetoreceptor cryptochrome specifically, such additional inhomogeneity could arise from different dimerization states [44][45][46] , different states due to rearrangements of the C-terminus of cryptochrome for binding signaling proteins [47][48][49] , or due to metabolites binding [50][51][52] , as illustrated in Fig. 1(d).
Variations in hyperfine interactions of the flavin cofactor due to changes in the protein environment, e.g. due to substrate binding, altered geometries of hydrogen bonds, surrounding amino acids have been shown experimentally [53][54][55][56][57][58][59][60][61] and theoretically 43,62,63 . If such inherent heterogeneities existed, they could potentially lead to significantly larger levels of variation than the 2% level assumed here. We therefore also investigate ensembles with larger levels of variations in hyperfine couplings.
By evaluating changes in the directional response of several radical pair magnetic compass systems, we investigate whether some nuclear spin environments lead to directional responses that are more robust to variations of hyperfine interactions than others. We find, interestingly, that some ensembles are dramatically less affected by such hyperfine variations than others. Figure 2 shows an ensemble with interesting, quadrumodal dependence on the geomagnetic field that is, however, a dramatic example of a non-robust system. If one considers variations of hyperfine couplings in the ensemble, the seeming directional sensitivity ends up being buried in the noise. In contrast, one can find ensembles that are much less affected by variations of hyperfine couplings. We find this observation to be true in radical pairs with spin-1/2 nuclei that have been used in many models in the literature as well as radical pairs with spin-1 nuclei, as they occur e.g. in nitrogens of flavin cofactors. Likewise, we see separation into robust and non-robust ensembles in two types of multinuclear spin systems with more realistic hyperfine parameters taken from flavin cofactor and tryptophan radicals. This qualitative conclusion does not depend on the exact size of variations assumed. Certainly, the effects on the magnetoreceptor responses become more pronounced for larger size of variations, but the observation of a robust class of ensembles remains that displays four or more orders of magnitude less signal uncertainty than non-robust cases (cf. Table 1). This significant difference is noteworthy as it might provide an evolutionary target by favoring ensembles that are robust to unavoidable hyperfine coupling variations. Moreover, we find that an intriguing interplay exists between optimal and robust radical pair ensembles and suggests that evolution can select radical pairs that are both optimal and robust.  obtained by averaging over all N anisotropic yields Φ n . Noise due to inhomogeneity leads to a signal uncertainty σ (or standard deviation), represented as a grey grid. (e) Section of the average singlet yield Φ for φ = 0, with average signal strength ∆Φ and uncertainty represented by σ. The average signal uncertainty (not shown here) 0 012). The number of radical pairs is N = 3000. Blue and yellow colors represent singlet yield values, respectively, larger and smaller than the spherical average. The radical pair model is the N-N model, where each radical is coupled to a nitrogen nucleus. The two hyperfine values are the same, correspond to N5 in flavin 30 , and are arranged perpendicularly. The hyperfine tensors components used for N5, in its principal axis, are A x = − 104.9, A y = − 99.6, A z = 1382.6 in μT (for details see ref. 30). The magnitude of variations is 5%, the lifetime of the radical pair is τ = 0.1 μs and the strength of the external magnetic field is B = 0.05 mT. The color bar reported in panel (a) applies to panels (a,c,d). Some anisotropic patterns have been enlarged by quoted scale factors relative to the first pattern (from the left) in panel (c).

Theory.
A single radical pair comprises a pair of radicals created by photoinduced electron transfer in a spin-correlated singlet or triplet state. In its basic model 5 each radical has an unpaired electron spin coupled to an external magnetic field via Zeeman interaction and a few local nuclei via hyperfine interactions, described by the following spin Hamiltonian: where θ φ θ φ θ = B (sin cos , sin sin , cos ) is the external magnetic field, with θ and φ being the polar and azimuth angles, which define the direction of B with respect to the radical pair frame. γ e is the electron gyromagnetic ratio, Ŝ A , Ŝ B the electron spin operators, Î Ai and Î Bj the nuclear spin operators, and A i and A j the hyperfine coupling tensors of the radical A and B. Dipolar and exchange interactions are neglected as is typical for most radical pair compass studies (but see ref. 27).
A simplified radical pair reaction model combines coherent quantum spin dynamics of singlet-triplet interconversion, driven by the spin Hamiltonian Ĥ , with the kinetics of spin-selective reactions, in which singlet and triplet pairs decay into distinct singlet and triplet reaction products 64 .
Following the Haberkorn approach, we describe spin-selective reactions by using phenomenological rate equations with first-order rate constants k S and k T for the singlet and triplet decay 65 . We point out that other approaches based on the theory of quantum measurements are currently under discussion [66][67][68][69][70] . The time evolution of a radical pair reaction is thus governed by the following phenomenological master equation: where P S and P T are, respectively, the singlet and triplet projection operators. We consider a spin-independent decay rate k S = k T = k 25,71,72 , and assume that a radical pair is created in a singlet state, ρ(0) = P S /Tr[P S ]. The singlet yield Φ , i.e. the amount of products decaying via the singlet channel, is calculated according to where ν mn = ν m − ν n are energy gap between pairs of eigenstates |n〉 and |m〉 of the spin Hamiltonian Ĥ . M is the number of nuclear spin states and = P mP n mn S S , with P S being the singlet projection operator. k is the spin-independent rate constant capturing all processes leading to a loss of spin correlation. We call τ the lifetime of the radical pair, which is given by τ = k −1 .
We do not explicitly model spin relaxation mechanisms.  Table 1. One-HFI radical pair models with different percentage-based (2%, 5% and 10%) and fixed-width (0.02, 0.06 and 0.1 mT) standard deviations Δ of normally distibuted hyperfine variations. Signal strength ∆Φ, robustness in logarithm scale σ log( ) and signal-to-noise ratio S/N of selected radical pairs from the optimal, plateau and least robust regions, for (a) the H-model and (b) the N-model. For both models, the selected radical pair from the optimal region has hyperfine strength a = 5.0 mT and axiality α = 1.0, from the plateau region has strength a = 5.0 mT, and axiality α = 0.5 for the H-model and α = − 0.5 for the N-model, from the least robust region has strength a = 0.05 mT, and α = 0.4 for the H-model and α = − 0.22 for the N-model (see Fig. 3). In all calculations the lifetime is τ = 10 μs, and the geomagnetic field is set to B = 0.05 mT. The signal-no-noise ratio is obtained considering ∆Φ as the signal S, and σ as the noise N.

Model.
We consider an ensemble of N radical pairs arranged to form perfectly ordered crystals with all radicals aligned and immobilized within photoreceptor cells. In a homogeneous ensemble, as generally considered previously in the literature, all N radical pairs have identical hyperfine tensors {A n } = A o , n = 1, N, and thus the same singlet yield Φ n (θ, φ) = Φ o (θ, φ), and the same anisotropic pattern Λ o (θ, φ), which in turn it represents the directional response of the ensemble. In an inhomogeneous ensemble, hyperfine tensors vary between different copies of radical pairs. Lacking detailed information about the type of variations in hyperfine couplings in protein-cofactor systems, we model the variations of N hyperfine tensors as being drawn from a normal distribution around A o . In particular, for each hyperfine interaction included in the model, we draw a Gaussian distribution  for each hyperfine tensor component We call Δ the standard deviation of a Gaussian distribution. We compute N singlet yields Φ n (θ, φ) and anisotropic patterns Λ n (θ, φ), corresponding to the N normally distributed hyperfine tensors. The directional response of the inhomogeneous ensemble is the average singlet yield θ φ Φ( , ), obtained by averaging over the N individual singlet yields, with a corresponding average pattern Λ θ φ ( , ). The signal strength of the inhomogeneous ensemble is obtained as max m in . We call the standard deviation of the average singlet yield, σ(θ, φ), signal uncertainty, whose average over all directions σ is used to evaluate the overall uncertainty of the directional response of the inhomogeneous ensemble. We define the signal-to-noise ratio as S/N, and consider the signal S being the signal strength θ φ ∆Φ( , ) and the noise N being the signal uncertainty σ(θ, φ).

Estimates of hyperfine coupling variations.
It is not yet possible to perform quantum mechanics-molecular mechanics (QM-MM) simulations on a 100 microsecond -1 millisecond timescale so as to evaluate the effect of conformational changes on these timescales on hyperfine couplings of flavin cofactors in cryptochromes. Such calculations are furthermore complicated by the absence of detailed structural information of the active form of cryptochrome. Finally, experimental values of hyperfine interaction variations for the desired protein systems are not readily available. We therefore base our estimate on two different approaches.
First, we can compare hyperfine coupling values between copies of a flavin radical of the same protein when substrates are added through experimental measurements. While such measurements do not exist for cryptochromes from birds, we can use values available for other flavoproteins.
The spin density of the ring system (i.e. isoalloxazine moiety) of a flavin cofactor radical of flavoproteins, including cryptochromes 61 , is highly sensitive to change in the protein environment. This has allowed hyperfine spectroscopy methods to use flavin cofactor radicals as a spin probe [59][60][61][62][63] . Particularly, electron nuclear double resonance (ENDOR) studies of flavoproteins have shown changes of proton hyperfine couplings of the flavin's isoalloxazine moiety upon substrate binding 53,54,56 . The measured proton hyperfine couplings changed on average by about 2% 56 , 5% 54 and 10% 53 in different flavoproteins. The 2% change was for enzyme-substrate binding in Escherichia coli DNA photolyase, which is highly homologous to cryptochrome 75 .
Secondly, through QM-MM simulations hyperfine interactions of flavin radicals were estimated to fluctuate on short time scales (picoseconds) with an average magnitude of 0.06 mT 43 . The smallest magnitude was 0.02 mT while the largest 0.1 mT. These short fluctuations would be averaged out in a radical-pair-based compass magnetoreception 38 , but we use them as a lower-bound estimate for the size of the relevant long time fluctuations that are beyond the reach of QM-MM calculations.
Based on such estimates, we model variations of hyperfine interactions with Gaussian distributions and consider a percentage-based standard deviation of 2% and a fixed-width standard deviation of 0.06 mT.
In the following, we show the results for a Gaussian with a width of 2% of the initial hyperfine components A k o , k = x, y, z. We have performed additional calculations with percentage values of 5% and 10%, as well as with fixed widths of 0.02 and 0.1 mT for all components and find the results to be qualitatively unchanged (see Tangle 1,  Table 2). The only exception occurs for radical pairs with the limiting axiality value of α = 1, which will be discussed in more detail below (see Table 1). Figure 2 depicts the anisotropic pattern, signal strength, and signal uncertainty of a particular inhomogeneous ensemble of radical pair reactions and represents the rationale for this study. For comparison, Fig. 2(a) depicts an intriguing anisotropic pattern Λ o with quadrumodal directional sensitivity for a homogeneous ensemble of N radical pairs with identical hyperfine tensors. Figure 2(b) shows a section (azimuth angle φ = 0) of the corresponding singlet yield as a function of the polar angle θ, indicating a sizeable signal strength. When considering inhomogeneity, as expected in realistic biological environments, variations of hyperfine tensors lead to different and complex anisotropic patterns for each radical pair as illustrated in Fig. 2(c). The directional information from the inhomogeneous ensemble stems from the average of the singlet yields, shown (for φ = 0) in (e) and as a three-dimensional anisotropic pattern in (d). It is noteworthy that the average pattern has qualitatively different features from the homogeneous pattern, as the quadrumodal character is now averaged out. More importantly, there remains a significant random component, uncertainty in the signal, as illustrated in both Fig. 2(d,e) by showing the average pattern together with one standard deviation σ of the signal variation. In the particular radical pair ensemble chosen here, the signal uncertainty σ is of comparable size to the maximum angular sensitivity Scientific RepoRts | 6:35443 | DOI: 10.1038/srep35443 or signal strength θ φ ∆Φ( , ). Clearly, in order to achieve a large signal-to-noise ratio, a reliable compass sensor not only needs to maximize the signal strength, but also needs to minimize the signal uncertainty.

Results
In the following, we ask whether there are choices of hyperfine values that reduce the signal uncertainty while also resulting in large signal strengths, thus leading to a robust and optimal compass sensitivity. To evaluate the robustness of the whole signal we use the average signal uncertainty σ, obtained by averaging σ(θ, φ) over all directions. We study simple radical pair models with up to two hyperfine interactions and only axial anisotropy, thereby limiting the parameter space and computational cost, and vary hyperfine parameters in a range of values typically encountered in organic radicals. We also investigate the effects of inhomogeneity between radical pair copies in more realistic chemical compasses by studying two types of multinuclear radical pair systems with realistic hyperfine interactions taken from flavin and tryptophan radicals 30 . Moreover, in all calculations we use radical pair lifetimes of the order of microseconds as suggested by theoretical interpretations of experimental results 38 . In particular, we report results of the simple radical pair models with a lifetime of τ = 10 μs in Figs 3, 4 and 5, and of the more realistic models with a lifetime of τ = 1 μs in Fig. 6.
One-HFI radical pair models. First, we consider the simplest radical pair model that can account for directional magnetic field effects, in which only one of the two electron spins is coupled to a nucleus through an axial hyperfine interaction, whereas the other radical is devoid of hyperfine interactions. We consider a spin-1/2 nucleus as in hydrogen. This H(ydrogen)-model has been extensively studied as a proof of principle in the literature 25,36,39,73 . We also consider a N(itrogen)-model with a spin-1 nucleus, motivated by the fact that the largest hyperfine coupling in cryptochromes is of this type. The Hamiltonian in Eq. 1 for these models has i = 1, j = 0, and the hyperfine tensor is axial. We parametrize the axial symmetry as follows: A x = A y = a − aα and A z = a + 2aα with a being the strength and α the axiality of the hyperfine interaction. The axial symmetry allows us to set the azimuth angle φ = 0, without loss of generality. We consider an ensemble of N = 3000 inhomogeneous copies of radical pairs, where each radical pair has hyperfine tensor components drawn from Gaussian distributions with averages A x = A y = a(1 − α) and A z = a(1 + 2α) and standard deviations being the 2% of respectively A x and A z . We compute the signal strength ∆Φ and signal uncertainty σ of the inhomogeneous for average hyperfine coupling strengths a between 0.05 mT and 5 mT, and all possible values for the axiality α. The strength of the external magnetic field is set to B = 0.05 mT.
The signal strength ∆Φ as a function of a and α, and for a lifetime τ = 10 μs, is depicted in Fig. 3(a) for the H-model and (c) for the N-model. The two models show similar features, and we can identify three regions. The optimal region of maximum signal strength occurs for axiality α = 1, with a signal strength of ∆Φ = .
0 166, for a ≫ 0.2 mT, for the N-model, both in agreement with previous   studies on a single radical pair reaction 36,37 . For the H-model, a large plateau region, depicted in green, with a constant signal strength of ∆Φ = .
0 125, can be observed for sufficiently large hyperfine strengths (a ≫ 0.5 mT) and away from the special axiality values α ≠ − 2, 0, 1. A similar plateau region is harder to identify for the N-model, but may occur with a less pronounced transition in signal strength and a smaller constant signal strength of ∆Φ = .
0 04, or a ≫ 0.15 mT and − 1.0 < α < 0 (within the blue area). Outside of the maximum and plateau regions, the signal strength depends strongly on the axiality. In general, larger axiality enhances the signal strength, while smaller axiality reduces signal strength. Figure 3 shows the signal uncertainty σ, in logarithmic scale, as a function of the hyperfine parameters for (b) the H-model and (d) the N-model, respectively.
For the H-model, an increase in signal strength appears to correlate with a decrease in signal uncertainty at first sight, as indicated by the similarity in qualitative features in Fig. 3(a,b). However, when plotting a histogram of signal uncertainty σ versus signal strength values ∆Φ (see Fig. 4(a)), we find that there is no simple correlation between signal strength and signal uncertainty. Rather, a strongly diminished signal uncertainty can be observed for two specific values of signal strength, namely ∆Φ = .
0 25 and ∆Φ = . 0 125, i.e. the optimal signal strength and the signal strength in the plateau region, respectively. Similarly, in the N-model, we find robust regions of greatly suppressed noise -or signal uncertainty -to occur at the two signal strength values of ∆Φ = .
Keeping a lifetime of τ = 10 μs, we repeated the calculations with 5% and 10% variations, and observed qualitatively the same results. That is, we could identify in every case a plateau and optimal region in signal strength of models, in the optimal and plateau regions the signal uncertainty is greatly suppressed, leading to a more robust compass system. The Gaussian distribution of the hyperfine tensors has a percentagebased standard deviation of 2%. The lifetime of the radical pair is τ = 10 μs and the strength of the external magnetic field is B = 0.05 mT. In the case of a constant-based standard deviation of e.g. 0.06 mT, the signal strength of the optimal region, for both models, decreases and the signal uncertainty significantly increases, thus leading to a loss of robustness, as described in the next section, and reported in Table 1. the same size as with the 2%, i.e. 0.125 and 0.25 for the H-model and 0.04 and 0.166 for the N-model, and found that both of these regions coincided with a dramatic decrease in signal uncertainty. Table 1 reports the signal strength ∆Φ, signal uncertainty in logarithm scale σ log( ), and signal-to-noise ratio S/N of selected radical pairs for different percentage-based standard deviations Δ , for (a) the H-model and (b) the N-model. The selected radical pairs are representative for the optimal, plateau and least robust regions. Interestingly, in all calculations (see Fig. 4(a,d) and Table 1) we find that, while the signal strength varies only within a comparatively small range, the signal uncertainty changes by four to five orders of magnitude for different selections of hyperfine parameters. This difference is noteworthy because leads to radical pairs in the optimal and plateau regions with a signal-to-noise ratio of at least four order of magnitude larger than radical pairs in the least robust region (see Table 1).
Moreover, we repeated the calculations with 2% variations for lifetimes of 1 μs and 100 μs, and find the same qualitative results as with a lifetime of 10 μs. If the lifetime is decreased to 0.1 μs, we begin to see a reduction of signal strength. For the H-model, such a reduction in signal strength as a function of the rate constant can be analytically predicted. Indeed, the signal strength values in the plateau and optimal regions can be determined analytically in the limit of sufficiently large hyperfine interaction a ≫ ω and sufficiently slow rate constant a ≫ k, where ω = γ e B. The signal strength for the plateau region (α ≠ − 2, 0, 1) is , and for the optimal region (α = 1) is θ ω ∆Φ − +  k k ( ) 1/4 ( / ) /4 2 2 2 , which reduces, respectively, to 0.125 and 0.25 for a ≫ ω ≫ k, i.e. when the rate constant is slow enough to allow the geomagnetic field to affect the spin dynamics 76 .
The optimal case of α = 1. As mentioned above and in ref. 36, a one-HFI radical pair shows the maximal signal strength for the unique case of an anisotropy α = 1, i.e. when two components of the hyperfine tensor in its diagonal form satisfy A x = A y = 0. Regardless of whether this case can be actually realized in a biological cofactor, one can ask whether this optimal case also has a desirable low signal uncertainty. We find that the answer to this question depends on the details of the distribution of hyperfine couplings for inhomogeneous ensembles. If we assume that variations will vary by 2%, 5% and 10% of the initial hyperfine value of each component, then such variations will leave the lateral components unchanged at zero, A x = A y = 0. In this case, the optimal case is also robust, because is unaffected by inhomogeneous noise, leading to minimal signal uncertainty and maximum signal-to-noise ratio as seen in Figs 3(b,d) and 4(a,d) and in Table 1.
However, there is no a priori reason why hyperfine couplings should not also vary in x and y directions and using a percentage-based variation may therefore lead to an underestimation of inhomogeneous noise in the case where the hyperfine coupling (or a component of it) is very small to begin with. To model an alternative distribution, we considered one-HFI models with hyperfine couplings drawn from a Gaussian distribution with a constant (i.e. 0.02, 0.06 and 0.1 mT), rather than percentage-based standard deviation for all hyperfine tensor components. Table 1 summarizes the results of the one-HFI models for a percentage-based and a fixed-width standard deviation. We find that inhomogeneity has two effects on a radical pair ensemble centered around hyperfine interactions with α = 1 (i.e. optimal region). It reduces the signal strength, because the variations cause the hyperfine couplings to deviate from the optimal anisotropy of α = 1. In addition, it increases the signal uncertainty, and consequently reduces the signal-to-noise ratio. For example, comparing the 5% with the 0.06 variation, the signal strength in the H-model reduces to 0.2 and the signal uncertainty increases from σ − .  log( ) 6 0 to σ −  log( ) 2, with a concomitant reduction of about four order of magnitude signal-to-noise ratio. For larger variations such a loss of robustness can increase signal uncertainty until it becomes comparable in size to the signal strength. It is noteworthy that the plateau regions discussed above are insensitive to the exact form of hyperfine variations. Changing from a percentage-based to a constant standard deviation in the distribution has no effect on either signal strength and signal uncertainty, as long as the magnitude of the standard deviations is comparable (see Table 1). Both numerical and analytical evaluations bear out that the plateau regions with signal strength of ∆Φ = .
0 125 and ∆Φ = . 0 04 in the H-model and N-model, respectively, correspond to the minimal signal uncertainty, regardless of the details of hyperfine coupling variations. This "robustness" appears to be a special feature of the plateau regions in the one-HFI models. For non robust radical pairs, changing from a percentage-based to a constant standard deviation further increases the signal uncertainty, which in some cases becomes larger than the signal strength. These results suggest that selection of noise-suppressing radical pairs can be an evolutionary target for compass sensors. Two-HFIs radical pair models. We extend our study to the slightly more complex situation in which we consider two axial hyperfine interactions. We study two arrangements. In one arrangement, one radical has two hyperfine interactions while the other is devoid of hyperfine interactions, and the Hamiltonian is given by Eq. 1 with i = 2 and j = 0. We denote the model with nuclear spin-1/2, 2H-0 model, and with nuclear spin-1, 2N-0 model. In the other arrangement, each radical has one hyperfine interaction and the Hamiltonian is given by Eq. 1 with i = 1 and j = 1, and we denote the model with nuclear spin-1/2, H-H model, and with nuclear spin-1, N-N model. In all cases, we parametrize the axial symmetry of the hyperfine interactions as with the one-HFI models, i.e. A x = A y = a − aα and A z = a + 2aα with a being the strength and α the axiality. Moreover, we define ψ as the polar angle (which is not shown in Eq. 1) between the two hyperfine tensors A 1 and A 2 , and without loss of generality we consider the azimuth angle φ = 0.  Φ o (θ, Φ). The values of σ and ξ are both in logarithmic scale. There is a correlation between σ and ξ, with both decreasing or increasing. This indicates that for a non-robust compass system, in which the noise leads to a large signal uncertainty, the average pattern of an inhomogeneous ensemble is very different from the pattern of the corresponding homogeneous ensemble (see Fig. 2). Conversely, for a robust system average and homogeneous patterns are very similar. The model used is the H-model (see Fig. 3(b)). We find the same trend for all models studied here.
Scientific RepoRts | 6:35443 | DOI: 10.1038/srep35443 We randomly generate 10 5 samples of nuclear spin environments with average hyperfine strengths, a 1 , a 2 , ranging from 0.05 to 5 mT, axiality α 1 , α 2 from − 2 to 2, and the angle ψ between the two hyperfine tensors from 0° to 90°. As with the one-HFI models, for each nuclear spin environment we consider an ensemble of 3000 radical pairs, where each radical pair has two hyperfine tensors drawn from Gaussian distributions with average A xm = A ym = a m (1 − α m ) and A zm = a m (1 + 2α m ), m = 1, 2. Using a lifetime of τ = 10 μs and a magnetic field strength of B = 0.05 mT, we calculate the average signal strength ∆Φ and signal uncertainty σ for each nuclear spin environment and for all models. We depict the results of the effects of variations of hyperfine interactions with magnitude 2%.  Figure 4 depicts how this relationship is affected when another hyperfine interaction is added to the one-HFI models. We note that a crucial factor for both signal strength and signal uncertainty is whether or not one radical remains devoid of hyperfine interactions. For the 2H-0 model (see Fig. 4(b)), the qualitative features resemble those of the H-model (see Fig. 4(a)) with the most robust compass sensors being for signal strengths of 0.125 and 0.25. The signal uncertainty for these values in the 2H-0 model is slightly higher than in the H-model, albeit still much lower than at other signal strength values. Furthermore, the signal uncertainty is slightly reduced for most other values of signal strengths. The trade-off is that it appears harder to realize specialized cases with extremely low signal uncertainty and optimal signal strength in such radical pairs. The features of the 2N-0 model (see Fig. 4(e)) confirm the trends seen in the 2H-0 model, but also show some new features. Regions of very low signal uncertainty can barely be discerned at signal strength of 0.04 (plateau region) and of 0.22 (optimal region), but similarly low uncertainty regions occur for other signal strengths as well. Overall, signal uncertainty is reduced again for most values of signal strengths, trading off with an increase in signal uncertainty for the extremely low uncertainty regions. We also note that two nitrogen nuclei on the same radicals give rise to a compass sensor with a larger directional sensitivity than the one with only one nitrogen, in agreement with a previous study 37 . Analyzing the parameters of radical pairs with large signal strengths for models with two hyperfine interactions on the same radical, we find that optimal and robust sensitivity can be realized when a nuclear spin environment comprises two hyperfine interactions of average strengths a 1 , a 2 ≫ 0.05 mT, both with a strong anisotropy α α   1 1 2 , and nearly parallel axes ψ  0. In both models the robust sensitivity of 0.125 and 0.04 is realized when one radical has one hyperfine interaction in the plateau region of the one-HFI models, and the other has a much smaller hyperfine strength.
We repeated the calculations with a lifetime of τ = 10 μs and with larger variation of 5% and 10%, and find similar features as with the 2% variation. Table 2 reports the results of selected radical pairs from the optimal, plateau and least robust regions for the different percentage-based standard deviations Δ . The optimal and plateau regions remain the most robust with unchanged signal strength and signal uncertainty that slightly decreases for larger variations. Similarly to the one-HFI model, the noise-suppressing 2HFI-0 model radical pairs have the advantage of a signal-to-noise ratio of at least three order of magnitude larger than non robust radical pairs.
Differing from the one-HFI models (see Tables 1 and 2), when changing from a proportional to a constant standard deviation in the Gaussian distribution we observe similar features as in the proportional case. Although the signal uncertainty slightly decreases, we can find regions that combine large signal strengths (0.125, and 0.25) with extremely low signal uncertainty (see Table 2). As in the proportional case, such optimal and robust compass sensors can be realized with hyperfine axiality close or equal to one, i.e. α α   1 1 2 . The signal uncertainty of the plateau region slightly decreases as well when changing from a proportional to a constant standard deviation, and the least robust region has at least two order of magnitude less signal-to-noise ratio than optimal and plateau regions.
When both radicals have one hyperfine interaction each as for the H-H and N-N models (see Fig. 4(c,f), respectively), the situation changes. Directional sensitivity is generally reduced. Regions of low signal uncertainty appear to correlate with regions of low signal strength, but clearly such radical pairs do not result in a useful compass sensors.
The fact that there are much fewer hyperfine parameters for the H-H and N-N models that result in combined large signal strength and very low signal uncertainty, indicates that these models are much less suitable for a compass sensor. One exception appears to occur in the N-N model at a signal strength of about 0.11 that shows a signal uncertainty of σ = − . log( ) 4 3. Such exception remains for larger variations of 5% and 10%, where the signal strength remains unchanged to 0.11, and the signal uncertainty resulted, respectively, σ = − . log( ) 3 8, and σ = − . log( ) 3 5. Such radical pair is obtained with two hyperfine strengths much larger that 0.05 mT, axiality close to one, and parallel axis (e.g. a 1 = 0.9 mT, a 2 = 4.3 mT, α 1 = 0.9, α 2 = 1.0 and ψ = 10.0). Noteworthy, such a radical pair shows robust feature when changing to a fixed-width standard deviation. For a standard deviation of 0.02 mT, the signal uncertainty resulted σ = − . log( ) 3 6, for 0.06 mT σ = − . log( ) 3 2, and for 0.1 mT σ = − . log( ) 3 0. We also repeated the set of calculations with 2% variation for the following lifetimes τ = 0.1, 1, 100 μs, and obtained qualitatively analogous results.
Inhomogeneous vs. homogeneous ensemble. The uncertainty of the signal due to inhomogeneity provides a fundamental problem and the main purpose of our study is to investigate how it could be overcome in building a reliable magnetic compass sensor. However, it is worth pointing out that considering inhomogeneity in a realistic biological ensemble also raises a serious caveat for interpreting patterns of directional sensitivity based on a single set of hyperfine parameters alone, as was routinely done in the literature. This caveat is illustrated already in Fig. 2. The constant set of hyperfine parameters, neglecting inhomogeneity, indicates a quadrumodal directional signal of the magnetic compass with similar singlet yields in four directions, as shown in Fig. 2(a). However, when considering inhomogeneity, the average singlet yield of the inhomogeneous ensemble, shown in Fig. 2(d) displays no more quadrumodal character, but shows a bimodal pattern. Thus, any conclusions about e.g. behavior of an animal based on the homogeneous singlet yield would be invalid, at least for the currently used estimate for the variation strength of hyperfine interactions. Figure 2 shows an example with a particularly high signal uncertainty of σ −  log( ) 2. We now ask how such deviation of the average singlet yield θ φ Φ( , ) from the yield of the homogeneous ensemble Φ o (θ, φ) correlates with the signal uncertainty. We quantify the deviation between average and homogeneous singlet yield as ξ Plotting the signal uncertainty σ versus the deviation between average and homogeneous pattern ξ in Fig. 5, we find a near linear correlation in a double logarithmic scale up to log(ξ) = − 2. Figure 5 uses 10 5 ensembles of the H-model with different hyperfine parameters and 3000 inhomogeneous copies for each ensemble, as above. Similar results are found for the other models. Thus we find that for a large class of radical pair models, neglecting variation of hyperfine interactions may be rather misleading, because, when including such variations, the average pattern will differ significantly from the pattern of the homogeneous ensemble. The exception to this caveat occurs when a robust radical pair ensemble with low signal uncertainty is selected. In this case, the pattern does not change noticeably from copy to copy and one can well approximate it with the pattern of a homogeneous ensemble. This robust pattern has a bimodal form with a yield approximating Φ ≈ (1 + cos 2θ) in the limit of a ≫ ω ≫ k.

Multinuclear Radical Pair Systems
The one-HFI and two-HFIs radical pair models, despite being simplified models of realistic chemical compasses, have allowed to explore many nuclear spin environments for finding optimal and robust radical pair parameters. It is computationally prohibitively expensive to perform the same exhaustive study on more realistic radical pairs, with a larger number of hyperfine interactions. To limit the parameter space, we choose therefore to study the effects of inhomogeneity on two types of multinuclear radical pairs with eight hyperfine interactions and realistic hyperfine values of flavin and tryptophan radicals obtained from DFT calculations 30 .
One type of eight-nuclei radical pair system comprises the hyperfine values of two nitrogen nuclei (N5 and N10) of the flavin adenine dinucleotide radical (FADH • ), and of six hydrogen nuclei (H1, H2, H4, H5, H6, H7) of the tryptophan radical (Trp). For simplicity we call this type FADH • -Trp. The Hamiltonian for FADH • -Trp is given by Eq. 1 with i = 2, and j = 6. The other type of radical pair system comprises the hyperfine values of two nitrogen nuclei (N5 and N10) and six hydrogen nuclei (H5, H6, H8α, H8α, H1′ , H1′ ) of the flavin radical, while the other radical is supposed to be devoid of hyperfine interactions. We call this type FADH • -0. The Hamiltonian for FADH • -0 is given by Eq. 1 with i = 8, and j = 0. These two types of radical pair systems can be seen as an extension of the 1HFI-1HFI and 2HFIs-0 models, previously studied, when more nuclei are added to, respectively, both radicals or to only one radical and thereby allow comparisons with previous results. The choice of the nitrogen nuclei arises from studies that shown N5 and N10 being the predominant nuclei in flavin 37 , while hydrogen nuclei are chosen based on their strength.
To model variations of hyperfine interactions we use the method described in the Methods section. We consider the hyperfine tensor components taken from flavin and tryptophan radicals 30 as average values of Gaussian distributions with fixed-width of 0.06 mT. We choose a fixed-width of 0.06 G because this is the average size of variations of many hyperfine tensors in a flavin radical available from the literature 43 . For each hyperfine tensor component we generate 3000 samples of normally distributed values, and then compute 3000 singlet yields of each eight-nuclei radical pair system. We set the lifetime of the two types of radical pairs to τ = 1 μs, and the external magnetic field to B = 0.05 mT. Figure 6 reports (a) and (d) the homogeneous pattern Λ o , (b) and (e) the average pattern Λ, and (c) and (f) several individual patterns Λ n for (top) FADH • -Trp, and (bottom) FADH • -0. The deterioration of the signal due to inhomogeneous noise, or signal uncertainty σ, is represented by the grey grid in Fig. 6(b,e). The average uncertainty resulted σ = .
0 002 for FADH • -Trp, and almost one order of magnitude less, i.e. σ = . 0 01, for FADH • -0. Both values are within the range of uncertainties of the simplified models (Fig. 4). These results show that the size of the signal uncertainty is affected by the choice of hyperfine parameters even in more extended hyperfine coupling systems, and indicate that multinuclear spin systems are not always more robust to the deteriorating effects of the inhomogeneous noise. Furthermore, the fact that FADH • -Trp is more robust than FADH • -0 is in agreement with our previous results (see Fig. 4), that show that radical pairs with hyperfine interactions in both radicals are generally more robust to inhomogeneous noise.
Interestingly, the noise affects the shape of the anisotropy differently in the two systems. While in FADH • -Trp inhomogeneity changes the shape of (b) the average pattern and (c) the individual patterns with respect to (a) the homogeneous pattern, in FADH • -0 the shape of (e) the average and (f) individual patterns remains unchanged with respect to (d) the homogeneous pattern. The shape of the anisotropy is affected by inhomogeneous noise because of a loss of axial symmetry in one radical (given by N5 and N10) when hyperfine interactions are added on the other radical. This is further confirmed by the complex patterns depicted in Fig. 2 where the radical pair system lacks of axial symmetry. Conversely, in all axial systems studied in the one-HFI and two-HFIs models, the shape remains unaffected by inhomogeneous noise and all patterns, i.e. homogeneous, average and individuals, have similar shape to those of Fig. 6 (bottom).
Moreover, the average signal strength is larger in FADH • -0 (∆Φ = . 0 0835) than in FADH • -Trp (∆Φ = . 0 0019), as expected from previous studies 37 . While FADH • -0 results more directional sensitive, it is also less robust to inhomogeneous noise than FADH • -Trp. The larger signal strength of FADH • -0 compensates its smaller robustness, thereby resulting with a signal-to-noise ratio larger by a factor of ten than FADH • -Trp.
We repeated the calculations only for the maximum and minimum singlet yield with different fixed-width and percentage based variations, and calculated the signal strength ∆Φ, the signal uncertainty σ s and signal-to-noise ratio, for both models. Table 3 reports the results, and we note that signal uncertainty slightly decreases for larger variations, and reduces when changing to a percentage-based variations. Regardless of the form and magnitude of variations the signal uncertainty of FADH-0 is about ten fold larger than FADH-Tr. However, because FADH-0 has a larger signal strength, the signal-to-noise ratio of FADH-0 results ten fold larger than FADH-Tr.

Discussion
We have investigated the effects of variations of hyperfine interactions on the directional response of an inhomogeneous ensemble of radical pairs. Evaluating these effects on simple radical pair models, up to two axial hyperfine interactions, each with a wide parameter range, reveals that the nuclear spin environment has a significant influence on the effects of such variations. In fact, some nuclear spin environments can make a radical pair-based compass virtually immune to the influence of variations of hyperfine interactions. Such robustness is a design feature distinct from optimality. An optimal compass sensor can be found when the compass responds to the external magnetic field with a maximum directional sensitivity, whereas a robust sensor is one where the directional response remains unaffected by variations regardless of the magnitude of the directional sensitivity. Both features ultimately affect the signal-to-noise ratio of a magnetic sensory system. The strength of the directional sensitivity can be thought of as the signal and the change of the signal due to variations between copies in the ensemble can be thought of as a source of noise.
Intriguingly, the noise stemming from inhomogeneity is inextricably linked with the choice of a particular radical pair. In other words, evolution cannot change the signal of a radical pair independently from changing the robustness to this source of noise because selection of radical pair parameters will affect both signal and noise robustness at the same time.
The question arises whether it is possible to find radical pairs that are both optimal (i.e. show a maximum signal strength) and robust (i.e. minimize the signal uncertainty) at the same time. For the one-HFI models studied here, the optimal radical pairs are realized for axiality α = 1, i.e. all but one component of the hyperfine tensor, when expressed in its principal axis, equals zero. We find that this radical pair is robust to variations of the hyperfine interaction along its primary axis, but not robust to variations in directions perpendicular to it. Assuming that it is possible to realize the theoretical limiting case of α = 1, it would thus also be necessary to shield this radical pair from variations perpendicular to its primary axis so as to maximize the signal-to-noise ratio. By sacrificing a relatively modest amount of signal strength (e.g. a factor of two compared to the α = 1 for the one-Hydrogen model), one can identify a large robust region where the directional sensitivity is independent of the hyperfine parameters. Such robustness is a general feature for one-HFI radical pairs and can be understood through analytical calculations.
The relationship between optimality and robustness becomes, unsurprisingly, more complex for radical pairs with two hyperfine interactions. In keeping with earlier observations 37 , we find that the response is stronger if one radical is devoid of hyperfine interactions (2HFI-0 model). Moreover, it is more common to find robust radical pairs with large directional sensitivities in 2HFI-0 radical pairs than in 1HFI-1HFI radical pairs. However, even in the less ideal classes, e.g. in the N-N radical pairs, one can still find a sizable group of radical pairs with a directional response larger than 0.1 that are robust, providing no fundamental obstacles for evolutionary selection of radical pairs with favorable signal-to-noise ratios.
Surprisingly, more realistic multinuclear spin systems are not naturally more robust to inhomogeneous noise, and again the choice of hyperfine parameters affects both signal and robustness to inhomogeneous noise at the same time. Indeed, we find that one type of eight-nuclei radical pair system (FADH • -0) has a larger directional sensitivity but smaller robustness (of ten fold) than the other type (FADH • -Trp) of eight-nuclei system studied here. Overall the signal-to-noise ratio was larger in FADH • -0 than in FADH • -Trp.
As suggested in Fig. 1, different dimerization states, different states due to rearrangements of the C-terminus of cryptochrome for binding signaling proteins or due to metabolites binding, may lead to significantly larger variations of hyperfine couplings, greatly exacerbating the noise in an inhomogeneous ensemble. Of course, it may be possible that nature avoids such additional heterogeneities by favoring one form (e.g. the monomer form over the dimer form) for the design of functional magnetoreceptor ensembles. It might seem intuitive that evolution would favor such ensembles over ensembles involving multiple receptor states. However, a counterargument is that such differential states open up the possibility of adding potentially desirable control elements, such as feedback loops. This can be seen in another example of a receptor system operating near physical detection limits, bacterial chemotaxis. Instead of avoiding variation, nature here purposely varies the signal strength of its receptors by methylation to establish a negative feedback loop critical for the detection of weak chemical gradients 77 . This example is not unique, but phosphorylation and methylation, i.e. two processes that lead to receptor states with differing sensitivity provide in many signaling systems the key feedback step to make a receptor system work 24 . Given that little is known about how a potential magnetosensitive step in the cryptochrome photocycle is embedded in biochemical regulatory framework, it appears premature to dismiss the possibility that regulation  Table 3. Multi-nuclear radical pair systems. Signal strength ∆Φ, signal uncertainty σ log( ) s and signal-tonoise ratio S/N for different percentage-based (2%, 5%, 10%) and fixed-width (0.02, 0.06, 0.1 mT) standard deviation Δ of hyperfine variations, of (a) the FADH-Tr model, and (b) the FADH-0 model. The signal strength The signal uncertainty σ log( ) s is calculated by averaging the signal uncertainty of Φ max and Φ min . The signal-to-noise ratio S/N is calculated considering ∆Φ being the signal S, and σ log( ) s being the noise N. The lifetime is τ = 1 μs, and the geomagnetic field is set to B = 0.05 mT.
Scientific RepoRts | 6:35443 | DOI: 10.1038/srep35443 favors an ensemble with multiple receptor states. Beyond phosphorylation, a process occurring in cryptochromes, such different regulatory states could be realized also through different levels of metabolite binding, such as ATP. If such multiple receptor states exist in an ensemble, it is possible that the size of variations of hyperfine couplings is increased further from the level considered in our calculations. The observation of robust ensembles, put forth here, would provide a possible way for nature to avoid the conundrum of either having to favor a reduction of variations or of favoring the existence of additional regulatory states; it may be possible to realize both.
One might argue that the signal uncertainty does not pose a serious obstacle because using a larger ensemble would allow to reduce this noise by averaging. Indeed, one expects the noise to be reduced by a factor of N , where N is the number of copies of sensory molecules involved in producing a magnetic signal. However, one needs to recognize that at any given time only a small subset of sensory molecules is actively involved in signaling. For cryptochromes, with a molar extinction coefficient on the order of ε = 10 4 M −1 cm −1 and irradiation levels as used in bird magnetic orientation experiments, one can estimate that out of one million copies of cryptochromes in an outer segment of a photoreceptor cell, fewer than 100 would be activated by light every second. Reducing noise by averaging would thus either come at a high metabolic cost and ultimately be limited by the number of cryptochromes that could be expressed in a cell or would require averaging over longer times, thus slowing down the speed of the magnetic sensory system. Both of these undesirable effects could be avoided by selecting a sensory molecule with a small signal uncertainty. In a different sensory modality, vision, the remarkable sensitivity to a single photon is only possible because the sensory biomolecule rhodopsin has an extremely low probability of spontaneous isomerization 78 , i.e. an extremely low background noise level.
If the variations arise through thermal fluctuations, a certain degree of averaging will occur due to fluctuations on a time scale faster compared to the radical pair reaction times. Our treatment does not calculate the effect of such averaging. Such calculations are computationally expensive and are very uncertain without knowledge of the fluctuation sizes over different time scales. Our current study is only a first step into the investigation of noise through radical pair parameter variations. Instead of using Gaussian distributions to model variations, one would ideally evaluate the effect of particular conformational changes or aminoacid substitutions of a protein-radical pair complex on the hyperfine interactions. Combining such calculations with more detailed information would provide a better description of the variations in the radical pair ensembles and may allow one to arrive at experimentally testable predictions.
There are other sources of noises that need to be considered for a full analysis, such as molecular shot noise in later signal transduction steps 33 , noise arising from imperfect geometrical ordering of radical pairs 22,28,29 , and spin relaxation mechanisms 38 . Unlike other sources of noise, spin relaxation mechanisms and hyperfine inhomogeneity stem from protein motions. While relaxation processes stem from protein motions that are on time scales shorter than spin-correlation times, hyperfine inhomogeneity, discussed here, is realized from protein motions that occur at time scales larger than spin-correlation times. It would certainly be worth investigating whether the librational motion found to increase compass sensitivity 38,79 could occur in a protein that is also robust against hyperfine inhomogeneity or whether nature would have to "choose" between the two and, if so, which one would provide the greater benefit for compass magnetoreception.
Other interesting scenario can emerge by investigating possible realizations of hyperfine inhomogeneity in a biological environment. For example, it was shown that substrate binding affects the geometry of the flavin ring 80 , thereby raising the question whether different conformational states of the protein-cofactor complex would lead to different motional constraints of the flavin ring, realizing in this way inhomogeneity on short time scales.
If evolution over millions of years has optimized radical-pair-based magnetic compasses, what did it act on? Were ways found to better fix molecules in ordered geometrical arrangements? Was the signal transduction network improved to amplify the initial small signals from the primary sensors? Or was the biomolecule for the primary response improved to yield a maximal signal, possibly through evolving motions causing spin relaxation? Here, we add a new idea to the discussion: selection may find biomolecules that are particularly good in preventing the primary signal from deteriorating because of the inhomogeneity present in a biological environment.