Distribution of reaction yields

To derive the distribution of singlet yields for a finite number of single-radical-pair reactions, we first consider the case where N = 1. The outcome of this reaction depends on the singlet probability, $P_{\rm S}^\alpha \lpar {t_{\rm R}} \rpar$, of the radical pair at the time when it reacts, t _R, given the initial nuclear spin configuration, labelled α. The probability that this radical pair forms the singlet product is just

(1)$$P\lpar {\Phi_{\rm S}\lpar 1 \rpar = 1\vert {t_{\rm R},\alpha} } \rpar = P_{\rm S}^\alpha \lpar {t_{\rm R}} \rpar \equiv \left\langle {\Psi^\alpha \lpar {t_{\rm R}} \rpar } \right \vert \hat{P}_{\rm S}\vert {\Psi^\alpha \lpar {t_{\rm R}} \rpar } \rangle $$

where $\hat{P}_{\rm S}$ is the singlet projection operator and |Ψ^α(t _R)〉 is the state of the radical pair at time t _R given an initial nuclear spin configuration α. In order to compute the overall probability of forming the singlet product, we need to average over t _R and α by combining Eq. (1) with the probabilities that the radical pair reacts at time t _R, p(t _R), and has initial state α, P(α). Assuming there is no nuclear polarisation or coherence at the moment of formation of the radical pair, P(α) = 1/Z where Z is the number of nuclear spin configurations. For symmetric recombination (i.e. k _S = k _T = k in Fig. 1), $p\lpar {t_{\rm R}} \rpar = k\,{\rm e}^{-k t_{\rm R}}$ (Timmel et al., Reference Timmel, Till, Brocklehurst, Mclauchlan and Hore1998).

Integrating this combined probability over all possible combinations of t _R and α (i.e. using the law of total probability) gives the overall probability that a single radical pair forms the singlet product:

(2)$$\eqalign{P\lpar {\Phi_{\rm S}\lpar 1 \rpar = 1} \rpar = & \sum\limits_\alpha {\int_0^\infty {P\lpar {\Phi_{\rm S}\lpar 1 \rpar = 1\vert {t_{\rm R},\alpha} } \rpar P\lpar \alpha \rpar p\lpar {t_{\rm R}} \rpar {\rm d}t_{\rm R}}} \cr = & \sum\limits_\alpha {\int_0^\infty {P_{\rm S}^\alpha \lpar {t_{\rm R}} \rpar P\lpar \alpha \rpar p\lpar {t_{\rm R}} \rpar {\rm d}t_{\rm R}}} \cr = \;& E_{\alpha, \,t_{\rm R}}\lsqb {P_{\rm S}^\alpha \lpar {t_{\rm R}} \rpar } \rsqb = \Phi _{\rm S}.} $$

The probability that a single reaction results in the singlet product is simply the expectation value of the singlet probability of the radical pair, which is exactly the average singlet yield, Φ_S.

Using this result, we can now derive the distribution that describes the singlet yield for a finite number of reactions, Φ_S(N). For N reactions, because each reaction is independent and identically distributed, we have N identical trials, each with probability Φ_S. The probability that $N\phi _{ {\rm S}}$ of these reactions result in the singlet product is described by the binomial distribution ${{\cal B}}\lpar {N,\Phi_{\rm S}} \rpar$ (normalised by 1/N):

(3)$$\eqalign{ {\Phi _{\rm{S}}}\left( N \right)\sim{1 \over N}{{\cal B}}\left( {N,{\Phi _{\rm{S}}}} \right) & { \;\Rightarrow\; }P\,\left( {{\Phi _{\rm{S}}}\left( N \right) = {\phi _{{\kern 1pt} {\rm{S}}}}} \right) \cr & = \left( {\matrix{ N \cr {N{\phi _{\rm{S}}}} \cr } } \right)\Phi _{\rm{S}}^{N{\phi _{{\kern 1pt} {\rm{S}}}}}{\left( {1 - {\Phi _{\rm{S}}}} \right)^{N - N{\phi _{{\kern 1pt} {\rm{S}}}}}}.} $$

For large N (the regime used exclusively below), the central limit theorem states that this binomial distribution is accurately approximated by a normal distribution ${{\cal N}}\lpar {\mu, \;\sigma^2} \rpar$ with mean μ = Φ_S and variance σ ²:

(4)$$\eqalign{\Phi _{\rm S}\lpar N \rpar & \sim{{\cal N}}\lpar {\Phi_{\rm S},\Phi_{\rm S}\lpar {1-\Phi_{\rm S}} \rpar /N} \rpar \cr & {\Rightarrow} \;p\left( {\Phi_{\rm S}\lpar N \rpar = \displaystyle{X \over N}} \right)\;{\approx}\; \lpar {2\pi \sigma^2} \rpar ^{-1/2}\exp \left[ {-\displaystyle{{{\lpar {\Phi_{\rm S}-X/N} \rpar }^2} \over {2\sigma^2}}} \right].} $$

The width of this distribution is $\sigma = \sqrt {\Phi _{\rm S}\lpar {1-\Phi_{\rm S}} \rpar } /\sqrt N$. For p not close to 0 or 1, $\sqrt {p\lpar {1-p} \rpar }$ does not depend strongly on p and the majority of the uncertainty in Φ_S(N) is determined by the number of reactions, and scales as $1/\sqrt N$. Equation (4) is verified in the Supplementary Information, Section I using a kinetic Monte Carlo algorithm (Voter, Reference Voter, Sickafus, Kotomin and Uberuaga2007).

Information theory

Now that we can calculate the uncertainty in the singlet yield for a finite number of reactions, we can ask how efficiently a cryptochrome-based radical pair compass sensor would perform under low-light conditions. The magnetoreceptors are located in the retinas of both eyes, probably in photoreceptor cells (Mouritsen et al., Reference Mouritsen, Janssen-Bienhold, Liedvogel, Feenders, Stalleicken, Dirks and Weiler2004; Liedvogel and Mouritsen, Reference Liedvogel and Mouritsen2010; Niessner et al., Reference Niessner, Denzau, Gross, Peichl, Bischof, Fleissner, Wiltschko and Wiltschko2011, Reference Niessner, Gross, Denzau, Peichl, Fleissner, Wiltschko and Wiltschko2016; Bolte et al., Reference Bolte, Bleibaum, Einwich, Günther, Liedvogel, Heyers, Depping, Wöhlbrand, Rabus, Janssen-Bienhold and Mouritsen2016; Hore and Mouritsen, Reference Hore and Mouritsen2016; Günther et al., Reference Günther, Einwich, Sjulstok, Feederle, Bolte, Koch, Solov'yov and Mouritsen2018). In the following, each cell is assumed to contain multiple, identically aligned cryptochrome molecules (Lau et al., Reference Lau, Rodgers and Hore2012). The output of a cell would be an average of the singlet yields of all reaction cycles of its cryptochromes during a specified integration time. For a given orientation of the bird's head in the geomagnetic field (the ‘head direction’, θ), the curvature of the retina means that cryptochromes contained in cells at different positions within the retina will experience different magnetic field vectors, and so have different average singlet yields, Φ_S. Thus, information on the head direction is encoded in ${\bf \Phi} _{\rm S}\lpar {\Theta; N} \rpar = \lcub {\Phi_{\rm S}^i \lpar {\Theta; N} \rpar } \rcub$, the set of singlet yields from all of the cells (labelled i) spread across the retina (Ritz et al., Reference Ritz, Adem and Schulten2000). In the following we use upper case to denote random variables, e.g. Θ, and lower case (θ) to denote sample values taken from the set of possible elements of those random variables.

Unfortunately, we know so little about how the yield of the cryptochrome signalling state might be transduced, processed and interpreted (Dodson et al., Reference Dodson, Hore and Wallace2013) that it is very difficult to reach conclusions about the precision of the compass without making assumptions of dubious validity. To circumvent this difficulty, we use a Rate Distortion Theory approach (Lestas et al., Reference Lestas, Vinnicombe and Paulsson2010; Hilfinger et al., Reference Hilfinger, Norman, Vinnicombe and Paulsson2016) to calculate a best-case compass precision for any given radical pair.

According to the title of Shannon's seminal paper (Shannon, Reference Shannon1948), the field of information theory (Jaynes, Reference Jaynes1957; Ihara, Reference Ihara1993; Cover and Thomas, Reference Cover and Thomas2012), provides ‘a mathematical theory of communication’ which encompasses the fundamental topics of signal processing and data compression. Rate distortion theory is specifically the analysis of lossy data compression. The simplest model involves an encoder–decoder pair in which the input data is compressed (generally to afford efficient communication) and then reconstructed by the decoder to give an output. The finite amount of information in the input data means that there is some limit on the extent to which the data can be compressed in order that it can be perfectly reconstructed by the decoder. Otherwise the encoding involves some loss of information that cannot be undone by the decoder, and results in the output being an imperfect copy of the input. In this case, rate distortion theory allows one to calculate the expected error or difference between the input and the output.

Returning to the radical pair model, we can think of the external field direction (relative to the bird's head direction) as the source of data, which is encoded in the singlet yield of cryptochromes in receptor cells in the retina. This, in turn, is decoded by some unknown signal processing pathway to give an estimate of the head direction. If the available singlet yield signal is noisy due to a small photon flux, then the ‘compression’ of the magnetic information will entail some error in the bird's estimation of its head direction, which we can quantify using several general inequalities from rate distortion theory. A schematic representation of this process is given in Fig. 3.

Fig. 3.

A schematic representation of the information theory approach. (a) The source of directional information is the geomagnetic field. θ is the bird's head direction with respect to geomagnetic north. (b) The head direction determines the magnetic field vector experienced by cryptochrome molecules in cells distributed across the retina. This information is encoded in the singlet reaction yield of light-induced radical pairs within the cryptochromes. (c) Each receptor cell i experiences a different field direction and so has a different average singlet yield, Φ_Sⁱ(Θ; N), represented here as a visual modulation pattern (see later). (d) Φ_S(Θ; N), the singlet yield signal from the whole of the retina, is decoded in the bird's eye and brain using some unknown data processing algorithm. (e) This results in an estimate of the head direction. The mutual information between different steps in the process is indicated.

The key quantities in information theory are the differential entropy, h(X), which measures the uncertainty in the value of a random variable, X, based on its domain of possible values, and the conditional entropy h(X|Y) which measures this uncertainty when the domain of X is restricted by knowledge of Y. These quantities are defined in terms of the probability density functions p(x) and p(x, y) as (Shannon, Reference Shannon1948)

(5)$$h\lpar X \rpar = -\int {\,p\lpar x \rpar \ln p\lpar x \rpar \,{\rm d}x}, $$

(6)$$h\lpar {X\vert Y } \rpar = -\int {\int {\,p\lpar {x,y} \rpar \ln p\lpar {x\vert y } \rpar {\rm d}x{\rm d}y}}. $$

Given the mean (μ) and variance (ε²) of the probability density of X, an important result is that the entropy is maximal when the distribution of X is Gaussian (Cover and Thomas, Reference Cover and Thomas2012), i.e. when $X\sim{\cal N}\left( {\mu, { \epsilon} {^2}} \right) = {\left( {2\pi \epsilon{^2}} \right)^{ - 1/2}}\exp [ { - {{\left( {x - \mu } \right)}^2}/2 \epsilon {^2}} \bbig]$. Using Eq. (5), the entropy of X therefore satisfies

(7)$$h\left( X \right) \le h\left( {X\sim{\cal N}\left( {\mu, \epsilon {^2}} \right)} \right) = {1 \over 2}\ln \left( {2\pi \,\hbox{e} \epsilon {^2}} \right).$$

h(X) and h(X|Y) can be combined to compute the ‘mutual information’ about variable X available from a stochastic source of data Y (Shannon, Reference Shannon1948):

(8)$$I\lpar {X;Y} \rpar \equiv h\lpar X \rpar -h\lpar {X\vert Y } \rpar.$$

The narrower the distribution of X when conditioned by Y, the lower the conditional entropy h(X|Y) and so the larger the mutual information (Cover and Thomas, Reference Cover and Thomas2012). If X is entirely determined by Y, the information is maximal and is determined by the entropy of X. If, instead, X and Y are independent, h(X|Y) = h(X) and the mutual information is zero. Equation (8) is the basis of our derivation of a lower bound on the error in the head direction derived from the radical pair sensors.

Returning to our application, we are interested in how precisely the bird can estimate its head direction, Θ, which determines the orientation of the geomagnetic field experienced by the magnetoreceptors in the retina. We assume that the system represents a Markov chain $\Theta \; \to \; {\bf \Phi} _{\rm S}\to \,\hat{\Theta}$ in which the bird's estimate of its head direction, $\hat{\Theta}$, is based solely on the encoding of the geomagnetic field direction in the singlet yield data coming from magnetoreceptor cells distributed across the retina, Φ_S(Θ; N). Our aim is to quantify the mutual information between the two ends of the signalling pathway, $I\lpar {\Theta; \hat{\Theta}} \rpar$, represented by panels (a) and (e) in Fig. 3.

If the number of reactions N = 0 there is clearly no information available about the external magnetic field. Θ and $\hat{\Theta}$ are therefore independent and the mutual information is zero. As N increases, the amount of information available from the singlet yield signal increases as the signal-to-noise ratio improves (the noise scales as $1/\sqrt N$). In the limit N → ∞ (assuming noiseless signal transduction), the information about the magnetic field is perfectly encoded in Φ_S(Θ; N), so that the conditional entropy $h\lpar {\Theta \vert {\hat{\Theta}} } \rpar$ is minimised and $I\lpar {\Theta; \hat{\Theta}} \rpar$ is maximal, meaning that the head direction can be determined exactly.

However, as we do not know how Φ_S(Θ; N) is quantified, transmitted and decoded in the eye and the brain, we cannot quantify $h\lpar {\Theta \vert {\hat{\Theta}} } \rpar$ and therefore cannot calculate the mutual information $I\lpar {\Theta; \hat{\Theta}} \rpar$. To sidestep this difficulty, we use a bound on this mutual information applicable for Markov chains – the data processing inequality (Nielsen and Chuang, Reference Nielsen and Chuang2010; Beaudry and Renner, Reference Beaudry and Renner2011),

(9)$$I\lpar {\Theta; \hat{\Theta}} \rpar \le I\lpar {\Theta; {\bf \Phi}_{\rm S}\lpar {\Theta; N} \rpar } \rpar,$$

which states that the mutual information between Θ and $\hat{\Theta}$ is bounded by I(Θ; Φ_S(Θ; N)). Because the estimate of the head direction is downstream of the singlet yield signal, and because $\hat{\Theta}$ depends on Θ only via the singlet yield profile, any degradation in the input cannot be undone by decoding Φ_S(Θ; N). The significance of this is that, as we can compute the mutual information between Θ and Φ_S(Θ; N), we can compute an upper bound on the information between Θ and $\hat{\Theta}$, for any arbitrary decoding of the singlet yield signal to compute $\hat{\Theta}$. A simpler approach in which the Cramér–Rao lower bound is obtained from the Fisher information is inappropriate here because assumptions about the signal processing would be required to relate Φ_S(Θ; N) to $\hat{\Theta}$ (Lehmann and Casella, Reference Lehmann and Casella1998; Cramér, Reference Cramér1999).

In the following section, we use Eq. (9) together with well-known results from information theory (Jaynes, Reference Jaynes1957; Ihara, Reference Ihara1993; Cover and Thomas, Reference Cover and Thomas2012) to compute a lower bound estimate of the mean-squared difference between Θ and $\hat{\Theta}$ in order to evaluate how the precision of the compass sensor depends on N and the properties of the radical pair.

Lower bound error in the head direction

We start with the mutual information between Θ and $\hat{\Theta}$ (Eq. (8)):

(10)$$I\lpar {\Theta; \hat{\Theta}} \rpar \equiv h\lpar \Theta \rpar -h\lpar {\Theta \vert {\hat{\Theta}} } \rpar.$$

The bird can face in any direction with equal probability, i.e. the head direction is uniformly distributed, so that Eq. (5) gives

(11)$$h\lpar \Theta \rpar = -\int_{-\pi} ^\pi {\left( {\displaystyle{1 \over {2\pi}}} \right)\ln \left( {\displaystyle{1 \over {2\pi}}} \right){\rm d}\theta} = \ln \lpar {2\pi} \rpar.$$

The mutual information, Eq. (10), can be related to $h\lpar {\Theta -\hat{\Theta}} \rpar$ using

(12)$$h\lpar {\Theta \vert {\hat{\Theta}} } \rpar = h\lpar {\Theta -\hat{\Theta} \vert {\hat{\Theta}} } \rpar \le h\lpar {\Theta -\hat{\Theta}} \rpar.$$

The equality in Eq. (12) arises because we are in effect subtracting a constant in the conditional entropy, and Shannon differential entropy is invariant under translation (Cover and Thomas, Reference Cover and Thomas2012). The inequality is Shannon's lower bound which comes from the well-known concept that conditioning cannot increase entropy (Shannon, Reference Shannon1959).

Using Eq. (7), we can evaluate a bound on $h\lpar {\Theta -\hat{\Theta}} \rpar$ in terms of the variance ε² of the distribution of $\Theta -\hat{\Theta}$:

(13)$${\rm \epsilon} ^2 \ge \displaystyle{1 \over {2\pi e}}\exp \lpar {2h\lpar {\Theta -\hat{\Theta}} \rpar } \rpar.$$

Since the expectation value $E\lsqb {\Theta -\hat{\Theta}} \rsqb = 0$, ε² is the mean-squared difference between Θ and $\hat{\Theta}$:

(14)$$ {\rm \epsilon} ^2 = E \left[ \left\{ \left( \Theta -\hat{\Theta} \right) -E \left[ {\Theta -\hat{\Theta}} \right] \right\}^2 \right] = E \left[ {{\left( {\Theta -\hat{\Theta}} \right) }^2} \right].$$

Combining Eqs. (9) and (13) gives an expression for the lower bound, ${\rm \epsilon} _{{\rm LB}}^{}$, the root-mean-squared error in the head direction estimate based on the encoding of the true head direction in the singlet yield signal across the retina:

(15)$${\rm \epsilon} ^2 \ge {\rm \epsilon} _{{\rm LB}}^2 = \displaystyle{{2\pi} \over \hbox{e}}\exp \big [ {-2I\lpar {\Theta; {\bf \Phi}_{\rm S}\lpar {\Theta; N} \rpar } \rpar } \big ].$$

The real power of this method is that Eq. (15) gives the lower bound of ε in the limit that 100% of the information available in the singlet yield signal is retained. This is valid for any arbitrarily sophisticated data processing algorithm used to interpret the singlet yield signal, and so successfully avoids having to make any guesses about signal transduction and processing in vivo, instead assuming the optimal case. Considering the complexity of data processing known in, for instance, the visual system (Hubel, Reference Hubel1995), the ability to account for a highly developed but unknown decoding mechanism is an important one. The evaluation of ${\rm \epsilon} _{{\rm LB}}^{}$ is detailed in the Supplementary Information, Section II.

Magnetoreceptor cells

The next stage is to define the orientations of the magnetoreceptor molecules in the retina with respect to the geomagnetic field. We model a bird's eye as a sphere with the retina spanning a hemispherical section of its surface and the pupil positioned on the optical axis opposite the centre of the retina. Although a crude approximation to reality, this geometry captures the essential physics. The cells that contain the magnetoreceptor molecules are assumed to be identical, uniformly arrayed across the retina, and cylindrically shaped (in the manner of the photoreceptor cells) with their symmetry axes pointing towards the centre of the sphere. Within each cell, all magnetoreceptor molecules, and therefore all radical pairs formed within them, are taken to have identical orientations. For simplicity, the optical axis of the eye is aligned parallel to the horizontal component of the geomagnetic field (θ = 0). The geomagnetic inclination (with respect to the horizontal) is taken as ϕ = 66°. In this model, therefore, the orientation of a radical pair with respect to the geomagnetic field vector is completely determined by the position in the retina of the cell in which it resides. We specify this position by co-latitude (κ) and azimuth (ν) angles and, for simplicity, assume a geometry within each cell such that the angle, ζ, between the magnetic field vector and the main axis of the radical pair is independent of the co-latitude and given by (see Supplementary Information, Section III)

(16)$$\cos \zeta = \sin \lpar \nu \rpar \sin \lpar \phi \rpar.$$

Evaluation of the minimum head direction error, Eq. (15), requires the diagonalisation of a covariance matrix (see Supplementary Information, Section II) whose dimension is the total number of magnetoreceptor cells in the retina (potentially > 10⁶). To make the problem tractable, we reduce the dimension of the system by grouping the cells contained in spherical caps of half-angle ξ = 0.5° to form a smaller number of ‘grouped magnetoreceptor cells’. The cells within each group have such similar positions and orientations that they will have almost identical average singlet yields, Φ_S (Worster et al., Reference Worster, Mouritsen and Hore2017). Each group is treated as a single receptor cell that absorbs the total number of photons that would impinge on the constituent cells during the integration time. The retina is thus described by a square grid of these 1°-spaced grouped receptor cells. In the visual system, similar groupings of photoreceptor cells form under dim light conditions to improve contrast at the expense of resolution (Warrant, Reference Warrant1999). As the increase in the number of photons per receptor compensates the reduction in the number of receptor cells, we consider this simplification to be valid. To verify this, we calculated ε_LB using ξ = 1.0°, representing a four-fold decrease in the number of grouped receptor cells each of which received four times as many photons. The calculated lower bound errors agreed to within a few percent, justifying the approximation.

The double-cone photoreceptor cells have been proposed as the most likely location of the magnetoreceptors (Worster et al., Reference Worster, Mouritsen and Hore2017; Günther et al., Reference Günther, Einwich, Sjulstok, Feederle, Bolte, Koch, Solov'yov and Mouritsen2018). They are arranged on an approximately hexagonal lattice in the retina with a nearest neighbour separation of 7 µm (Kram et al., Reference Kram, Mantey and Corbo2010), giving an average of one double cone per 40 µm² of retina. Taking r = 4 mm as the typical radius of the eye of a small songbird (Thomas et al., Reference Thomas, Szekely, Cuthill, Harper, Newson, Frayling and Wallis2002), the area of a spherical cap with half-angle ξ = 0.5° is 2πr ²(1 − cosξ) = 3800 µm². There are therefore ~100 double cones in each grouped receptor cell and ~30 000 grouped receptor cells in the model retina.

We now need to estimate the number of photons, N, that enter each grouped receptor cell in a given time interval. Birds can orient using their magnetic compass on cloudless, moonless nights at light intensities down to 3 × 10⁻⁴ lux corresponding to an average of ~0.3 photons incident on each 1 µm² of retina per second (Cochran et al., Reference Cochran, Mouritsen and Wikelski2004; Hore and Mouritsen, Reference Hore and Mouritsen2016). A typical photoreceptor cell has diameter ≈2 µm and therefore occupies ~3 µm², meaning that each cell receives roughly 1 photon per second under these conditions. Since there are ~100 cells per grouped receptor cell, we have N ≈ 100 photons per grouped receptor cell per second.

The signal-to-noise ratio of the sensor can be improved by increasing the period during which photons are absorbed and integrating the resulting signals (Warrant, Reference Warrant1999). Following Cochran et al. (Reference Cochran, Mouritsen and Wikelski2004), it seems unlikely that this integration period would be much longer than a few minutes, meaning that the number of photons per grouped receptor cell is in the range 10² ⩽ N ⩽ 10⁴ when the light intensity is ~3 × 10⁻⁴ lux.

A model radical pair

The singlet yield signals produced by grouped receptor cells spread across the whole of the retina can be pictured as ‘visual modulation patterns’ (VMP) (Ritz et al., Reference Ritz, Adem and Schulten2000) in which the singlet yield for each grouped receptor cell is represented here by a colour-coded pixel in a red-to-blue heat map. To obtain a two-dimensional pattern, the retina is projected onto a plane tangent to the centre of the retina using the pupil as the point of perspective. This stereographic projection maps the position of a grouped receptor cell onto a point in this plane defined by Cartesian coordinates (x, y):

(17)$$\eqalign{{\rm tan}\lpar {\kappa /2} \rpar =\; & \sqrt {x^2 + y^2}, \cr {\rm tan}\lpar \nu \rpar = \;& x/y.} $$

The co-latitude and azimuth are thus represented, respectively, by radial and angular features in the VMP. The centre of the retina (κ = 0) projects to the centre (x = y = 0) of the VMP (Lau et al., Reference Lau, Rodgers and Hore2012; Worster et al., Reference Worster, Mouritsen and Hore2017).

Before applying the information theory approach to realistic radical pairs, we begin with a simple parameterised model of the dependence of the average singlet yield on the direction, ζ, of the geomagnetic field with respect to the radical pair (see Supplementary Information, Section III),

(18)$$\Phi _{\rm S} = \displaystyle{1 \over 2}-A\exp \left[ {-\displaystyle{{{\lpar {\zeta -90^\circ} \rpar }^2} \over {2\sigma^2}}} \right].$$

This form of Φ_S is an inverted Gaussian of amplitude A, taken as A = 1.27 × 10⁻³, and width σ, centred at ζ = 90°, superimposed on a constant background of amplitude 0.5. The value of A was chosen to match that found in simulations of a 16-spin model of a cryptochrome-based radical pair with a lifetime, τ = 1/k = 1 µs (Hiscock et al., Reference Hiscock, Worster, Kattnig, Steers, Jin, Manolopoulos, Mouritsen and Hore2016). In this simple example, the alignment of the radical pairs in the receptor cells is such that the signal is independent of the co-latitude of the cell (Eq. (16)) and therefore has no radial variation in the VMP (Fig. 4d).

Fig. 4.

Calculation of the error in the estimate of the head direction for a model radical pair magnetic compass sensor. (a) Visual modulation patterns for different numbers, N, of photons per grouped receptor cell for a singlet yield anisotropy with σ = 26° where each pixel represents $\Phi _{\rm S}^i \lpar {\Theta, N} \rpar$ for a single grouped receptor. (b) An example of the processing of Φ_S(Θ; N) for N = 10⁴. The orange trace is the diameter-averaged signal from the VMP plotted as a function of the azimuth. The blue trace is the result of fitting these data to a Gaussian model of the singlet yield anisotropy. (c) Comparison of the root-mean-squared error for the averaging and fitting methods with the information theory lower bound, ε_LB as a function of the number of photons per grouped receptor cell. The inaccessible region where ε < ε_LB is shaded red. (d) VMPs in the N → ∞ limit for different values of σ. (e) Lower bounds calculated as a function of the width of the Φ_S signal, σ.

Figure 4a gives an example of how the singlet yield signal across the retina varies for different values of N, when σ = 26°. As anticipated, the noise arising from the variance of $\Phi _{\rm S}^i \lpar {\Theta, N} \rpar$ tends to obscure the signal for the smaller numbers of photons per receptor cell, N. Each 100-fold increase in N improves the signal-to-noise ratio by a factor of ~10 (Eq. (4)) and, by eye, the signal only emerges from the noise when N ⩾ 10⁶.

Before considering the information theory lower bound, we assess its achievability using two simple data processing methods to compute the root-mean-squared error in the derived head direction. In the first, we average the noisy signal along the diameters of the VMP and obtain a crude estimate of the head direction, $\hat{\Theta}$, as the value of the azimuth that corresponds to the minimum of this processed signal, e.g. Fig. 4b (where $\hat{\theta} \approx 95$°). This procedure was repeated for many realisations of the singlet yield signal, Φ_S(Θ; N), so that the mean-squared error, ${\rm \epsilon} ^2\;{\rm = \;} E{\rm [(}\Theta -\hat{\Theta} {\rm )}^2{\rm ]}$ can be calculated. The second, more sophisticated, data processing method entails fitting a Gaussian model function to the diameter-averaged signals to obtain an improved estimate of the head direction from the value of the azimuth at the minimum of the fitted model function (Fig. 4b). Again, this is repeated many times, and the mean-squared error in the head direction computed. These calculations have been constructed such that both processing methods give the correct result when N → ∞.

Figure 4c compares the information theory lower bound, ε_LB, with the root-mean-squared error ε derived from the above two methods. In all three cases, the precision improves as N increases, as expected. The curve-fitting method offers a relatively good proof of principle that the information theory bound is not much smaller than that achievable by data processing. Indeed, some more sophisticated processing method could achieve a better agreement, but the importance of the information theory result is that it provides a bound. The more primitive interpretation technique, based on selecting the minimum value of the diameter-averaged signals gives a far larger error.

The same model of the singlet yield (Eq. (18)) can be used to assess how the width, σ, of the signal affects the precision of the compass (Fig. 4d). This is of interest in relation to recent simulations of cryptochrome-based radical pairs which predict that for lifetimes exceeding ~1 µs, Φ_S contains a sharp feature, referred to as a spike or ‘needle’ (Hiscock et al., Reference Hiscock, Worster, Kattnig, Steers, Jin, Manolopoulos, Mouritsen and Hore2016), that could deliver a much more precise compass bearing than previously thought possible. Figure 4e shows ε_LB as a function of σ for N = 100 and 1000. When the number of photons is limited, there is an optimum width that minimises the error in the estimated head direction. The more photons there are, the narrower is the optimal width of Φ_S. It appears that for small N, more spatial averaging is required to extract directional information, and very narrow signals are less well suited to this coarse-graining.

Cryptochrome-based radical pairs

At least three radical pairs have been discussed as potential magnetic sensors in cryptochromes in vivo. [FAD^•− TrpH^•+] is formed by light-induced electron transfer from a tryptophan (TrpH) amino acid residue within the protein to the fully oxidised FAD (flavin adenine dinucleotide) chromophore and gives rise to magnetic field effects on long-lived states of the purified protein (Maeda et al., Reference Maeda, Robinson, Henbest, Hogben, Biskup, Ahmad, Schleicher, Weber, Timmel and Hore2012; Kattnig et al., Reference Kattnig, Evans, Déjean, Dodson, Wallace, Mackenzie, Timmel and Hore2016; Sheppard et al., Reference Sheppard, Li, Henbest, Neil, Maeda, Storey, Schleicher, Biskup, Rodriguez, Weber, Hore, Timmel and Mackenzie2017). Alternatively, a radical pair could be formed during the re-oxidation of the fully (photochemically) reduced state by electron transfer from FADH⁻ to an unknown electron acceptor (Niessner et al., Reference Niessner, Denzau, Stapput, Ahmad, Peichl, Wiltschko and Wiltschko2013, Reference Niessner, Denzau, Peichl, Wiltschko and Wiltschko2014). We consider two possibilities, both of which have been proposed previously: [FAD^•− Z^•] in which Z^• is a radical devoid of hyperfine interactions (Hogben et al., Reference Hogben, Efimova, Wagner-Rundell, Timmel and Hore2009; Ritz et al., Reference Ritz, Wiltschko, Hore, Rodgers, Stapput, Thalau, Timmel and Wiltschko2009; Solov'yov and Schulten, Reference Solov'yov and Schulten2009; Müller and Ahmad, Reference Müller and Ahmad2011), and [FAD^•− Asc^•−], formed when ascorbic acid is the electron donor (Lee et al., Reference Lee, Lau, Hogben, Biskup, Kattnig and Hore2014; Nielsen et al., Reference Nielsen, Kattnig, Sjulstok, Hore and Solov'yov2017). Simulations suggest that both could produce much larger anisotropic magnetic field effects than [FAD^•− TrpH^•+] and could therefore form the basis of a more sensitive direction sensor. With the hyperfine interactions used here (Supplementary Information, Section IV), only [FAD^•− TrpH^•+] is predicted to exhibit the narrow spike mentioned above (Hiscock et al., Reference Hiscock, Worster, Kattnig, Steers, Jin, Manolopoulos, Mouritsen and Hore2016). In the following, we abbreviate these three radical-pair states of cryptochrome as FAD-Trp, FAD-Z and FAD-Asc. To compare their performance under low-light conditions, we use the same intracellular orientation of the cryptochromes as in Fig. 4, any choice being arbitrary given our ignorance of the situation in vivo (see Supplementary Information, Section III). The dependence of Φ_S on the direction of a 50 µT magnetic field was calculated as described previously using the hyperfine interaction tensors listed in the Supplementary Information, Section IV.

Figure 5a shows the lower bound error for FAD-Trp with lifetimes of 1, 2 and 5 µs. As the lifetime increases, the singlet yield signal shows both an increase in anisotropy and the appearance of the compass ‘needle’ (Hiscock et al., Reference Hiscock, Worster, Kattnig, Steers, Jin, Manolopoulos, Mouritsen and Hore2016) for τ ⩾ 5 µs. It is therefore not surprising that ε_LB is smaller for the longer lived radical pairs. When τ for FAD-Trp is extended to 100 µs, the spike becomes much sharper and ε_LB ≈ 0 for N ⩾ 100 (not shown). By contrast, FAD-Z, because of its much larger magnetic anisotropy, has ε_LB ≈ 0 for N ⩾ 100 and τ ⩾ 1 µs (also not shown). For a given lifetime, FAD-Z offers a much more precise estimate of the head direction than does FAD-Trp.

Fig. 5.

The lower bound error in the head direction for the FAD-Trp radical pair state of cryptochrome as a function of the number of photons per grouped receptor cell. (a) A truncated 16-spin model with no exchange or dipolar coupling for lifetimes τ = 1, 2 and 5 µs. (b) The full 29-spin model with τ = 1 µs. The shaded regions indicate the expected number of photons per grouped receptor cell (10² ⩽ N ⩽ 10⁴) when the light intensity is ~3 × 10⁻⁴ lux.

The results shown in Fig. 5a were obtained using the previously studied 16-spin model of FAD-Trp (details are given in Supplementary Information, Section IV) (Lee et al., Reference Lee, Lau, Hogben, Biskup, Kattnig and Hore2014), neglecting the exchange and dipolar interactions between the radicals. Both approximations greatly reduce the computational difficulty of the spin dynamics calculations but lead to overestimates of the anisotropy of Φ_S (Lewis, Reference Lewis2018).

A calculation that includes all 27 hyperfine-coupled nuclear spins in FAD-Trp together with realistic exchange and dipolar interactions is not feasible using an exact quantum mechanical approach. However, Lewis (Reference Lewis2018) has shown that approximate calculations (with an estimated 5% statistical error) can be performed on spin systems of this size using asemiclassical algorithm (Manolopoulos and Hore, Reference Manolopoulos and Hore2013). As shown in Table 1, there is a decrease in the anisotropy of the signal of over an order of magnitude in going from 14 nuclear spins to 27. As expected, this decrease is mirrored in the calculated lower bound error (Fig. 5b) which is large even for N ⩾ 10⁵ photons per grouped receptor when τ = 1 µs. This suggests that there may be a serious problem with a FAD-Trp-based compass sensor under low-light conditions.

Table 1.

Lower bound head direction errors

^a Magnetic parameters are given in Supplementary Information, Section IV. N = 1000.

^b Number of nuclear spins included in the FAD radical and the counter radical, respectively.

^c Anisotropy of the singlet yield, defined as ΔΦ_S = max[Φ_S] − min[Φ_S].

Finally, we turn to FAD-Asc and model a rapidly tumbling Asc^•− radical as the partner of the FAD^•− radical (details are given in Supplementary Information, Section IV). Compared to the 16-spin model of FAD-Trp, the singlet yield anisotropy is larger by a factor of about 50 when τ = 1 µs. This again leads to a very small head direction error (Table 1), and suggests that an organic counter-radical with small hyperfine interactions could prove to be a more sensitive compass sensor in vivo. However, both FAD-Asc and FAD-Z are expected to show a similar decrease in anisotropy as noted above for FAD-Trp pair on inclusion of the radical–radical interactions.