Sensing Magnetic Fields with Magnetosensitive Ion Channels

Abstract

Magnetic nanoparticles are met across many biological species ranging from magnetosensitive bacteria, fishes, bees, bats, rats, birds, to humans. They can be both of biogenetic origin and due to environmental contamination, being either in paramagnetic or ferromagnetic state. The energy of such naturally occurring single-domain magnetic nanoparticles can reach up to 10–20 room $k_{B}T$ in the magnetic field of the Earth, which naturally led to supposition that they can serve as sensory elements in various animals. This work explores within a stochastic modeling framework a fascinating hypothesis of magnetosensitive ion channels with magnetic nanoparticles serving as sensory elements, especially, how realistic it is given a highly dissipative viscoelastic interior of living cells and typical sizes of nanoparticles possibly involved.

1. Introduction

Influence of weak electromagnetic fields on living species is perceived by many scientists as a controversial subject matter. Nevertheless, there is a huge body of evidence of a substantial impact ([1,2,3,4], see, especially, the book by Binhi [5] and the references therein). One of such manifestations is given by the microwave auditory effect or Allan Frey hearing effect [6,7,8,9,10,11,12,13,14,15,16], an auditory perception of microwave pulses by humans and animals, which earlier has been considered mysterious. Now, the mystery of this effect is completely resolved within a thermoelastic theory [7,8,9,10,11,12,13,14,15,16] of acoustic wave production in closed resonators (e.g., human or animal head) filled with microwave absorbing tissues having a very large water content (think about heating of food in microwave oven, to realize a possible physical reason). Good reviews are available [13,14,15] and the theory and experiment agree convincingly well. The energy absorption per pulse of $16\mathsf{\mu }\mathrm{J}/\mathrm{kg}$ sufficient to produce the microwave hearing effect [9,13] in humans is 36,000 times lower that the maximal limit of 576 J/kg permitted in the IEEE C95.1 radiation safety standard [13], and a corresponding pulse-like elevation of temperature is really tiny, about ${10}^{-6\circ }\mathrm{C}$ per pulse [13,15], however, rapid (about $\mathsf{\mu }$s). This is currently probably the only one of known profound effects of weak electromagnetic fields on living systems which is explained completely. However, a direct influence of GHz and THz waves on neuronal tissues, which is also evidenced experimentally, is still not convincingly explained and this is the subject of ongoing research [17,18]. Epidemiological evidence for follow-up health effects including a spectrum of neuropsychiatric disorders is extensive, see e.g., in [19], and the references therein. Sensing and navigation of various living species such as magnetosensitive bacteria, fishes, turtles, bees, bats, rats, birds, etc. in the weak magnetic field of the Earth (about 50 $\mathsf{\mu }$T) presents another well established effect [5,20,21]. Differently from electric fields, quasi-static magnetic fields are practically not screened by moving ions and counter-ions, and can deeply penetrate into biological tissues [5]. Currently, such a high sensitivity to magnetic fields presents a puzzle with two basic hypothetical theories proposed to resolve it. One suggests that this is a non-thermal quantum effect based on spin-dependent electron transfer reactions [20,22,23,24] in certain proteins related to vision [21,23,24,25] and it is light-dependent. Whereas it might indeed be relevant for certain birds, it certainly cannot be applied to animals navigating in dark such as fishes, or nocturnal mammals such as bats. An alternative theory is based on a wide presence of magnetic nanoparticles (MNPs), such as magnetite (iron oxide, ${\mathrm{Fe}}_3{\mathrm{O}}_4$), in organs and tissues of many biological species [26,27,28,29], including human brain [30,31,32,33,34]. The biomagnetite nanoparticles were first discovered in magnetotactic bacteria [35], where they form chains of magnetosomes [35,36,37]. Whereas magnetite nanoparticles produced by, e.g., diesel engines, have typically a rounded form, the biomagnetite particles have typically a prismatic or elongated parallelepiped form [26,27]. This characteristic feature is even used as a plausible criterion for that magnetite nanoparticles found in a Martian meteorite ALH84001 may possibly evidence for the existence of life on this planet in the past [38]. The presence of biomagnetite is currently accepted by NASA as one of possible traces of life. How big can such particles be? For example, a typical volume size corresponding to the maximum of probability distribution of biomagnetite particles in a bacterial strain named NIc is $V=144\times 106\times 106{\mathrm{nm}}^3$ [37]. Crystals of such size are in a single-domain ferrimagnetic state at room temperatures [26,27] with saturation magnetization $M_s=4.8·{10}^5$ A/m. Hence, such particles have a permanent magnetic dipole moment $\mu =M_{s}V=7.28\times {10}^{-16}\mathrm{A}·{\mathrm{m}}^2$ and their magnetic energy in the magnetic field of the Earth is $E_B=\mu B_e\approx 3.88\times {10}^{-20}\mathrm{J}\sim$9.47 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{T}}_{\mathrm{r}}$, i.e., almost 10 room ${\mathrm{k}}_{\mathrm{B}}\mathrm{T}$, well over a typical thermal energy. Interestingly, about 10% of biomagnetite particles in human brain have sizes corresponding to such a large magnetic energy [30] with overall particles density 4 ng per gram of tissue in the gray matter on average, or about ${10}^6$ particles per gram of tissue. In hippocampus, their density is even larger, 50 ng/g [32]. A recent PNAS paper [33] questioned the biological origin of this magnetite. Indeed, most particles found in the polluted brains of deceased persons who lived in Mexico City and Manchester were of anthropogenic origin [34]. They were of a rounded form with a median diameter of 18 nm. However, their specific density ranged from 0.2 to 12 $\mathsf{\mu }$g/g of dry tissue [33], i.e., it was from 50 to 3000 times larger than one in unpolluted human brains [30]. No wonder that biomagnetite concentration is such brains was negligible with respect to one incurred by the environmental pollution. As a side remark, a brain polluted by magnetite nanoparticles will absorb more incident microwave energy than one free of such particles, or having a natural concentration. In particular, the local temperature around such particles can be increased stronger. A standard argumentation [39,40] that a temperature increase by ${10}^{-4}-{10}^{-5\circ }\mathrm{C}$ is physically and biologically irrelevant is wrong, as evidenced by the theory and practice of microwave auditory effect [13,15]. This is because the temperature increase is very rapid at a pulse-modulation. Magnetic nanoparticles are known to heavily absorb microwaves between 0.5 and 10.0 GHz through ferromagnetic resonance and also produce hypersound via magnetoacoustic effect [41,42]. The physics of this effect is even used to thermally destroy cancer cells, which are known to accumulate magnetic nanoparticles [43], by using sufficiently strong microwaves (magnetic hyperthermia). An excess accumulation of magnetic nanoparticles in biological tissues can generally be a signature of decease.

Currently, magnetic control of cellular processes using biofunctional MNPs presents a hot topic (see a recent review [44] and the references therein). Apart from hyperthermia-based therapy and controlled drug delivery [43], this includes, e.g., genetically targeted control of neuronal system [45], control of calcium influx in cortical neuronal networks [46], and control of inner ear hair cells [47]. In [45], a genetically engineered magneto-sensitive ion channel named Magneto was created. It comprises of a cation ion channel, TRPV4, fused to a paramagnetic protein ferritin. In [46], mechano-sensitive ion channels were controlled via MNPs modulating tension in biological membranes. In [47], MNPs were fused with cilia of hair cells. In all cases, sufficiently strong magnetic fields were required, of the order of 100 mT, i.e., three orders of magnitude larger than $B_e$. This is because the corresponding MNPs were either (super)paramagnetic or smaller than 100 nm in linear size, or both. All three above cited papers are already dealing experimentally in fact with artificial magneto-sensitive ion channels or complex nanomagnetic biostructures involving ion channels. In addition, in Ref. [48], a membrane pore forming activity of magnetic nanoparticles has been shown. Hence, such man-made biological structures were already in fact demonstrated, however, for sufficiently strong magnetic fields, much stronger than biological species normally experience.

The idea that a magnetic nanorod can play a role in biological sensing of weak magnetic fields by birds has first been proposed by Yorke [49]. Kirschvink et al. [27,28,50] suggested that it can be a magnetosensitive ion channel that involves a magnetic nanoparticle as sensory element. Several variants of such channels were further proposed and discussed [51,52]. Whether such theoretical proposals can be feasible or not requires, however, a serious detailed investigation of the dynamics of such models [40], which is necessarily stochastic. The sensory element unavoidably experiences friction and random thermal noise caused by the environment, which are related by the fluctuation–dissipation theorem [53,54,55], at a local thermal equilibrium. Cytosol is a viscoelastic liquid, when it is functional, to a first rough approximation, with the main water component which accounts for up to 80% of the cytosol’s mass content. However, it is densely stuffed with different protein polymers, which dramatically enhances the cytosol viscosity for particles of the linear size of 100 nm range and even smaller. Thus, in Refs. [28,29,50], the effective viscosity felt by magnetic nanoparticles is assumed to be 100 times larger than the one of water. Notice that a starving cell can do a transition to an anabiosis state, where cytosol behaves more like a superviscous solid with virtually infinite viscosity [56]. However, we are more interested in its functionally active liquid-like state. The characteristic time scale of sensor was estimated it Refs. [28,50] assuming that the sensor is monostable and it fluctuates around a fixed point, which is not affected by the magnetic field. In that original model, one assumes that the ion channel opens when a critical angle fluctuation occurs, and the amplitude of this fluctuation is affected by the magnetic field. However, most biological ion channels do exhibit a characteristic bistable dynamics while fluctuating between open and closed states [57]. Binhi and Chernavskii considered an orientational bistable dynamics of a magnetosome tethered to cytoskeleton [58,59], however, not in the context of ion channels, but rather stipulating that the above quantum mechanism can be mediated by a fluctuating magnetic field of a magnetosome. Indeed, it can largely exceed one of the Earth [60,61]. Anisotropic field of a spherical ferromagnetic magnetosome is estimated to be up to 402 mT strong near to its surface [61].

The magnetic sensor dynamics of hypothetical ion channels should also be bistable, rather than monostable. Such a model was proposed recently in Ref. [61]. The bistability therein is induced by a gating spring type instability as earlier suggested in the context of hair cell ion channels [62,63]. The analysis of this model for realistic parameters showed, however, that for a viscous friction that is 100 times larger than one in water the time scale of switchings would be so large that such a channel would not be functional. Moreover, the effective friction caused by cytosol for the particles of the size of 100 nm can be even larger, e.g., 1000 fold larger than one in water [64,65,66,67,68,69]. Cytosol as a complex fluid [70,71] is, however, not a Newtonian but rather viscoelastic liquid [72,73,74,75,76,77,78,79,80,81,82,83,84,85], and on the appropriate time scales (probably up to hours in some cases) it is characterized by a slowly, power law decaying memory kernel with a memory cutoff at large times [70,71,86,87]. Integration of this memory kernel yields an effective friction at very large times, when the memory effects can be neglected. The discussed memory friction yields subdiffusion on the relevant time scales, which has been experimentally measured for various nanoparticles in cytosol [56,73,74,77,78,79,82,88,89,90,91,92,93,94,95,96,97,98,99] including magnetosomes [81]. Even smaller particles of the size of only several nanometers can subdiffuse on the time scale up to one hour [74,83]. The non-Markovian dynamics, which includes such effects, has also been studied in Ref. [61] using a Markovian embedding approach of Refs. [87,100]. Contrary to naive reasoning involving a largely enhanced normal viscous friction, however, in accord with the results of non-Markovian rate theory [101,102,103,104,105], it has been shown that such bistable sensors can be functional and operate on a millisecond to second time scale. This is in line with some earlier studies showing that viscoelastic subdiffusion largely accelerates (and not hinders, contrary to a common but misleading interpretation [84]) transport processes in living cells over the naive macroscopic Markovian treatment with a largely enhanced normal viscosity. Viscoelastic power-law memory friction yields non-exponential distribution of the waiting times such as stretched exponential distribution, which elegantly explains [61,87,100] the physical origin of such distributions in ionic channels [106,107,108] and $1/f$ noise [109], although other approaches also exist [109,110,111,112,113,114,115,116]. In [61], nanosensor rod consists of several nanoparticles. It is coupled by peptide elastic springs to the gating structural elements of several ion channels forming a cluster and can do a large-amplitude orientational motion (about 150 angle degree change), while moving to a metastable state corresponding to the open state of the channels in cluster. In this paper, I will explore the possibility of sensor consisting of the only one sufficiently large nanoparticle and doing a relatively small orientational motion (about 30 degree change) while creating an opening torque on the gates of the channels within a very similar model. It will be shown that such a magnetic sensor is more realistic and it would operate much faster, on biologically relevant time scales, despite a largely enhanced effective friction. Viscoelastic properties of cytosol are very important for this and cannot be disregarded.

2. Model

We consider the model sketched in Figure 1, where a single magnetosome consisting of an elongated nanoparticle of magnetite in single-domain ferrimagnetic state of length L and width $d<L$, dressed in a protein–lipid membrane, can rotate around one edge fixed on a cytoskeleton element arming the cell membrane inside the cell. For $L=150$ nm and $d=107$ nm its magnetic energy in the magnetic field of the Earth reaches $E_M=\mu B_e=4.12·{10}^{-20}$ J or about $10k_B{T}_r$ for $T_r=297$ K depending on its orientation, which is characterized by the angle $\varphi$, and the orientation of the magnetic field given by the angle $\psi$ in plane geometry. Biomagnetite nanoparticles of such a size are commonly met both in some bacteria and in the human brain. This “nanocompass” is coupled by m elastic peptide linkers (depicted in red) to the gates (depicted in blue) of m ionic channels forming a cluster in the membrane (one is depicted). One end of the linker is attached at the distance l from the rotation axis to magnetosome and another one to a molecular latch, which forms a gate opening and closing an ion-conducting pore formed by a membrane channel protein. In the absence of magnetic field or for its unfavorable orientation (e.g., $\psi =0$), the linkers are relaxed (folded) and the channels are closed (right part of Figure 1). For a properly oriented field, the linkers are fully stretched (unfolded) and the channels are open (left part of Figure 1). Even in predominantly closed state channels can stay open time from time due to thermally activated transitions of sensor between its two metastable states depicted in Figure 1. Likewise, in the predominantly open state channels close stochastically in time. The mean time intervals in the states and the mean opening probability of ion channel complex, which determines the ion current controlled by it, strongly depend on the magnetic field. Sensor experiences friction and thermal noise caused by environment, which crucially determine its stochastic dynamics.

The linker, which is an entropic spring provided by a disordered peptide, is modeled by a finite extensible nonlinear elastic (FENE) chain [117]. Its elastic energy depending on the elongation x is given by $U_{\mathrm{FENE}}\left(x\right)=-\frac{1}{2}{k}_L{l}_{\max }^2\ln (1-x^2/l_{\max }^2)$, where $k_L$ is elastic spring constant, and $l_{\max }$ is the maximal extension length of the linker, when it is fully stretched. The rotation of sensor is thus bounded to some angular interval $[0,{\varphi }_{\max }]$, where ${\varphi }_{\max }$ depends on l and $l_{\max }$. We will choose it sufficiently small, like in Figure 1. The channel gate can be in one of two states. The closed state is characterized by the energy ${ϵ}_1$, and the open one has the energy ${ϵ}_2-f_{0}x$, which depends on the linker elongation x, where $f_0$ is a force constant characterizing the strength of coupling (force exerted by the linker on the gate). The gate fluctuates very fast and its dynamics is slaved to a much slower sensor. Statistical mean force exerted by the channel gate on the linker is $f\left(x\right)=-dG\left(x\right)/dx$, where $G\left(x\right)=-k_{B}T\ln Z\left(x\right)$, is a corresponding potential of mean force with $Z\left(x\right)=\exp [-\beta {ϵ}_1]+\exp [-\beta ({ϵ}_2-f_{0}x)]$ being the statistical sum of the gate, and $\beta =1/\left(k_{B}T\right)$ is inverse temperature. The mean force is $f\left(x\right)=f_{0}p\left(x\right)$, where

$$\begin{array}{c}\hfill p\left(x\right)=\frac{1}{1+\exp [-f_0(x-l_0)/\left(k_{B}T\right)]},\end{array}$$

(1)

is probability of the gate to be open and $l_0=({ϵ}_2-{ϵ}_1)/f_0$. To define some $x_0$ as equilibrium point, we, following [63], redefine the mean force by a shift as $f\left(x\right)=f_0[p\left(x\right)-p\left(x_0\right)]$. Furthermore, the linker elongation is approximated as $x\left(\varphi \right)=2l[\sin (\varphi /2)-\sin ({\varphi }_0/2)]$, where ${\varphi }_0$ is an equilibrium angle. The potential of mean force or rather torque acting on the rod in our model is [61]

$$\begin{array}{ccc}\hfill U\left(\varphi \right)& =& -\frac{1}{2}k{l}_{\max }^2\ln \left\{1-{\left[x\left(\varphi \right)/l_{\max }\right]}^2\right\}\hfill \\ & -& k_{B}Tm\ln \left\{1+\exp [f_0(x\left(\varphi \right)-l_0)/\left(k_{B}T\right)]\right\}\hfill \\ & +& m{f}_{0}p\left({\varphi }_0\right)x\left(\varphi \right)-\mu Bcos(\psi -\varphi ),\hfill \end{array}$$

(2)

where $p\left({\varphi }_0\right)=p(x=2l\sin ({\varphi }_0/2))$, $k=m{k}_L$. We shall scale the energy in units of $U_0=k{l}_{\max }^2$, temperature in the units of $U_0/k_B$, distances in units of $l_{\max }$, and forces in units of $f_u=U_0/l_{\max }$. $U_0$ will be fixed to $U_0=10k_B{T}_r\approx 41\mathrm{pN}·\mathrm{nm}\approx 0.25\mathrm{eV}$. For $m=7$ and a linker with stiffness $k_L=0.0429\mathrm{pN}/\mathrm{nm}$ [95], $k\approx 0.3\mathrm{pN}/\mathrm{nm}$, this $U_0$ corresponds to $l_{\max }\approx 11.69$ nm and force units $f_u\approx 3.51$ pN. In this paper, we choose $l=2$, $l_0=0.91$, $f_0=3$, ${\varphi }_0=0.1rad\approx {5.73}^{\circ }$. The corresponding $U\left(\varphi \right)$ and $p\left(\varphi \right)$ are plotted in Figure 2. $U\left(\varphi \right)$ is bistable due to a gating spring instability featuring such models [62]. Namely, if the sensor pulls the linker sufficiently strong, another metastable state emerges. Notice that the probability of the channel to be half-open, $p=0.5$, belongs in this model to the attraction domain of $U\left(\varphi \right)$ belonging to the open state. When it becomes lower in energy than one corresponding to the relaxed linker (closed channel), the channel becomes predominantly open. This occurs, e.g., when the magnetic field is applied at the angle $\psi \approx {114.59}^{\circ }$.

3. Theory and Results

To obtain the averaged probability $\langle p(\psi ,B)\rangle$ of the ion channel cluster to be in the open state depending on the strength and orientation of the magnetic field, one has to average $p\left(x\right(\varphi \left)\right)$ in Equation (1) over $\varphi$ with the statistical weight function $P(\psi ,B,\varphi )=\exp [-U\left(\varphi \right)/\left(k_{B}T\right)]/\mathcal{Z}$, where $\mathcal{Z}={\int }_0^{{\varphi }_{\max }}\exp [-U\left(\varphi \right)/\left(k_{B}T\right)]d\varphi$ is the corresponding statistical sum. The result is shown in Figure 3

The corresponding averaged current through a cluster of ion channel is $\langle I\rangle =m{i}_0\langle p(\psi ,B)\rangle$, where $i_0$ is unitary conductance of a single channel in the cluster. Already single such sensory complex consisting of $m=7$ large conductance ion channels with $i_0\sim$50 pA can be sufficient to depolarize the membrane above a sensitivity threshold and evoke spiking activity in hypothetical magneto-sensitive neurons [61].

For the magneto-sensor complex to be functional, its dynamics is, however, also very important. For example, if its characteristic times would lie in the range of minutes or hours, it, for sure, would not be of any relevance for animals.

3.1. Stochastic Dynamics without Memory

First, we consider stochastic orientational dynamics of the sensor in the viscous medium under the mean torque $f\left(\varphi \right)=-\partial U\left(\partial \right)/\partial \varphi$, viscous friction term $F_v=-{\eta }_0\dot{\varphi }$, with orientational friction coefficient ${\eta }_0,$ and the corresponding white Gaussian thermal noise ${\xi }_0\left(t\right)$ of the environment. The last two are related by the fluctuation–dissipation relation (FDR) named also the second fluctuation–dissipation theorem (FDT) by Kubo [53,54,55], $\langle {\xi }_0\left(t^{\prime }\right){\xi }_0\left(t\right)=2k_{B}T{\eta }_0\delta (t-t^{\prime })$, at the medium’s temperature T. $\delta \left(t\right)$ here is the Dirac’s delta function signaling that this noise has an infinite root-mean squared amplitude—a common singular model in statistical physics. The overdamped stochastic dynamics reads [54,55,118]

$$\begin{array}{c}\hfill {\eta }_0\dot{\varphi }=f\left(\varphi \right)+{\xi }_0\left(t\right),\end{array}$$

(3)

and a characteristic time-scale entering it is ${\tau }_{sc}={\eta }_0/U_0$. Precise estimation of the rotational frictional coefficient for the particle of the form considered is not easy [119]. A simplest estimate can be obtained by replacing it with the sphere of equal volume $V=L{d}^2$. Then, ${\eta }_0\sim$6${\zeta }_{0}V$, where ${\zeta }_0$ is the medium’s viscosity. For water at T = 20 ${}^{\circ }$C with ${\zeta }_0\sim$ 1 $\mathrm{mPa}·\mathrm{s}$, we obtain ${\tau }_{sc}\sim$ 0.17 ms, which is much smaller than for the rod-like sensor in [61]. This is a first reason why it is faster. ${\tau }_{sc}$ is the time scale used in our simulations done using the second-order stochastic Runge–Kutta method, or stochastic Heun method [120], see in Methods. The sample trajectories are shown in Figure 4 both for predominantly closed channels (part a, for $\psi =0$), and for predominantly open channels (part b, $\psi =2$ rad.)

From very long single trajectories we extract residence time distributions in open and closed states using the following procedure. Two thresholds are placed at the potential minima of the sensor metastable states corresponding to $U\left(\varphi \right)$ or, equivalently, two maxima of the probability distribution of the ionic current. The residence time interval in the closed state starts after a first downward crossing of the lower threshold and continues until the first upward crossing of the higher threshold. Likewise, the residence time interval in the open state starts after a first upward crossing of the upper threshold and continues until the first downward crossing of the lower threshold. This allows to map the continuous fluctuation processes in Figure 4 (central part) onto the corresponding two-state, on-off processes, characterized by the survival probabilities of the residence time-intervals in the corresponding states, $P_c\left(\tau \right)$ and $P_o\left(\tau \right)$. In the case of sufficiently large potential barriers, the Markovian continuous state dynamics yields also two-state Markovian dynamics completely characterized by the exponential survival probabilities $P_{c,o}\left(\tau \right)=\exp (-r_{o,c}\tau )$, and the corresponding probability densities ${\psi }_{c,o}\left(\tau \right)=-d{P}_{c,o}\left(\tau \right)/d\tau =r_{o,c}\exp (-r_{o,c}\tau )$, where $r_o=1/\langle {\tau }_c\rangle$, and $r_c=1/\langle {\tau }_o\rangle$ are the opening and closing rates, respectively. In the case of Markovian two-state dynamics, they are inverse of the mean residence times in the closed state, $\langle {\tau }_c\rangle$, and the open state, $\langle {\tau }_o\rangle$, correspondingly. As a more general model, we use a stretched-exponential or Weibull distribution

$$\begin{array}{c}\hfill P_{c,o}\left(\tau \right)=c_{c,o}\exp \left[-{(\tau /{\tau }_{c,o})}^{{\beta }_{c,o}}\right]\end{array}$$

(4)

for the whole survival probability derived from numerics (then, $c_{c,o}=1$) or some parts of it. The corresponding density

$$\begin{array}{c}\hfill {\psi }_{c,o}\left(\tau \right)=\frac{c_{c,o}{\beta }_{c,o}}{{\tau }_{c,o}^{{\beta }_{c,o}}}\frac{1}{{\tau }^{1-{\beta }_{c,o}}}\exp \left[-{(\tau /{\tau }_{c,o})}^{{\beta }_{c,o}}\right]\end{array}$$

(5)

has a decaying power law part $1/{\tau }^{1-{\beta }_{c,o}}$, for $0<{\beta }_{c,o}<1$ which is the reason why this distribution can be confused for a truly power-law dependence on the corresponding plots for a large ${\tau }_{c,o}$. To avoid this pitfall of interpretation, the plots of survival probability $P\left(\tau \right)$ can be preferred.

The just described procedure of extracting two-state process is well known (see, e.g., in [61,100]). In [61], the derived in this way exponential distributions agree very well with the results of Kramers theory for the transition rates [105], which confirms that this procedure is essentially correct. In the present case, the potential barriers are smaller and the curvatures of the potential wells (at the bottom), and the barrier (at the top) are very different. In such a situation, a good agreement with the Kramers theory is not expected. In this respect, notice large fluctuations of the sensor orientation in Figure 4 in the case of the relaxed linker and the closed state of channel. In a sharp contrast, these fluctuations are much smaller when the linker is in its tense state. This does not mean, however, that there cannot be large fluctuations of conductance when the linker is tense. Vice versa, they are much larger when the linker is tense than when it is relaxed. This is because the midpoint of $p\left(\varphi \right)$ dependence belongs to the attraction basin of the metastable state of $U\left(\varphi \right)$ that corresponds to the open state, and the barrier corresponds to $p\left(\varphi \right)\approx 0.2$ (see Figure 2). The amplitude of current fluctuations which corresponds the large-amplitude fluctuations of sensor with the relaxed linker is quite small because of an exponential dependence in Equation (1). The numerics for the case of predominantly closed channels in Figure 5a show that the closed time residence time distribution is nearly exponential $c_c=1$, and ${\tau }_c\approx \langle {\tau }_c\rangle$. However, the open time distribution is not quite exponential. Initially, it is different from exponential, see in inset of Figure 5a, which is the reason why $c_o\approx 1.214>1$ and fitted ${\tau }_o$ does not agree well with $\langle {\tau }_o\rangle$. The corresponding mean opening probability is calculated as $\langle p\rangle =\langle {\tau }_o\rangle /(\langle {\tau }_o\rangle +\langle {\tau }_c\rangle )$, and it is quite small in Figure 5a. For predominantly open channels in Figure 5b, both open and closed time distributions are nearly exponential. The corresponding mean times $\langle {\tau }_o\rangle =4.254\times 0.17\approx 0.723$ ms and $\langle {\tau }_o\rangle =2.338\times 0.17\approx 0.397$ ms are quite small, being typical for very fast ion channels like sodium channels, which are crucial for neuronal excitability. Notice also that the corresponding $\langle p\rangle \approx 0.645$ is somewhat larger than $\langle p\rangle$ calculated for the continuous process in Figure 2.

Separation of Closed and Open States with a Single Threshold

The distribution of current values in the right panel of Figure 4 implies, however, that it might be difficult practically implement the detecting procedure with two thresholds because the lower threshold must be set at a very tiny, ${10}^{-8}-{10}^{-6}$, value of current. What a theorist can easily do in a Gedankenexperiment, an experimentalist can have difficulties to realize in practice. Setting the lower threshold at some other higher value, say $0.2$, can very essentially modify the thus derived statistics [61], and, hence, it is a rather arbitrary procedure. Another common procedure, which some experimentalists implement both in ion channel research [57], and, especially, in deriving statistics of blinking quantum dots [121,122,123,124] is to operate with the only one detection threshold. It is “naturally” placed at the minimum of the current distribution separating two maxima. An immediate objection of a theorist experienced in the rate theory is that this procedure corresponds to placing a separation threshold on the top of potential barrier (for the sensor, in our case) separating two metastable basins of attraction. It might first look natural. However, physical picture of the rate transitions between two metastable basins of attraction says that it is not. In fact, the particle dwells mostly in a potential well and seldomly, with rate $r_0=({\omega }_0/2\pi )\exp [-\Delta U/\left(k_{B}T\right)]$, where ${\omega }_0$ is circular frequency of oscillations near the bottom of the potential well, and $\Delta U$ is the height of the potential barrier, comes on the barrier top. However, it generally can dwell for a while on the top of this potential barrier and re-cross this single threshold many times, before the particle will make finally transition to another potential well. These multiple crossings yield the so-called transmission coefficient which being multiplied with $r_0$ renders the resulting rate much smaller. The whole rate theory is, roughly speaking, about how to calculate this transmission coefficient, which depends on friction, etc. Its maximal value is one (single crossing). This is why for the overdamped dynamics considered, this second procedure yields totally different residence time distributions, which severally distort the correct distributions characterizing two-state dynamics. In particular, the mean residence times derived in such a way will be much smaller than the correct ones. Experimentally, the problem is softened and can be masked by a finite time resolution $\Delta t_{res}$ of a measurement device. This is why not every fast recrossing will be measured. With $\Delta t_{res}$ becoming large many such recrossing events will be missed. However, this does not mean that statistics will become closer to the correct one. Not at all. To model this second procedure we use $\Delta t_{res}=100\delta t$, where $\delta t$ is the time step used in simulations, and the corresponding results are depicted in Figure 5c. First, the mean residence times are indeed much smaller (they became naturally even much more smaller for $\Delta t_{res}=\delta t$). Second, the residence time statistics is severely distort. In particular, about 90% of the closed time intervals follow now a spurious power law, $P_c\left(\tau \right)\propto 1/{\tau }^{\gamma }$, with $\gamma \approx 0.437$, and hence ${\psi }_c\left(\tau \right)\propto 1/{\tau }^{1+\gamma }=1/{\tau }^{1.437}$. Similar power laws are indeed measured for quantum dots using a similar approach with one threshold (however, the physics there is very different and our reasoning cannot be directly applied). Moreover, a spurious stretched exponential tail appears with ${\beta }_c=0.762$ and weight $c_c=0.0834$. The open time distribution has, however, an exponential tail, ${\beta }_o=1$. Our example shows explicitly how dangerous can be this second detection method and why it cannot be trusted. Interestingly, this second procedure barely affects $\langle p\rangle$, cf. $\langle p\rangle =0.642$ in Figure 5b vs. $\langle p\rangle =0.626$ in Figure 5c.

Notice that with a naíve replacement ${\tau }_{sc}\to 100{\tau }_{sc}=17$ ms in a medium with effective viscosity 100× larger than one of water our sensor dynamics would become 100× slower what would essentially deteriorate its functionality. A much slower model channel in [61] would even cease to be of any interest in biological context, if to apply this Markovian reasoning with ${\eta }_{\mathrm{eff}}$ = 100 ${\eta }_0$. This, however, does not happen in fact in viscoelastic media such as cytosol, where a more careful treatment is required, since the viscoelastic memory effects become very essential.

3.2. Stochastic Dynamics in Viscoelastic Environment

In viscoelastic crowded environments such as cytosol, apart from viscous friction caused by the its main water component, also a viscoelastic memory friction is present, $F_{v-el}\left(t\right)={\int }_0^t{\eta }_{\mathrm{mem}}(t-t^{\prime })\dot{\varphi }\left(t^{\prime }\right)d{t}^{\prime }$, where ${\eta }_{\mathrm{mem}}\left(t\right)$ is a memory kernel. It is necessarily complemented by a corresponding correlated thermal unbiased Gaussian random force of the environment, ${\xi }_{\mathrm{mem}}\left(t\right)$. In accordance with the (second) FDT, $\langle {\xi }_{\mathrm{mem}}\left(t\right){\xi }_{\mathrm{mem}}\left(t^{\prime }\right)\rangle =k_{B}T{\eta }_{\mathrm{mem}}\left(\right|t-t^{\prime }\left|\right)$. Sensor dynamics in this case obeys a generalized Langevin equation (GLE) reading

$$\begin{array}{c}\hfill {\eta }_0\dot{\varphi }=f\left(\varphi \right)-{\int }_0^t{\eta }_{\mathrm{mem}}(t-t^{\prime })\dot{\varphi }\left(t^{\prime }\right)d{t}^{\prime }+{\xi }_{\mathrm{mem}}\left(t\right)+{\xi }_0\left(t\right).\end{array}$$

(6)

The simplest Maxwellian model of viscoelasticity corresponds to an exponentially decaying memory kernel, ${\eta }_{\mathrm{mem}}=k_1\exp (-{\nu }_{1}t)$, where $k_1$ is a spring constant and ${\nu }_1$ is a relaxation rate of stress. For ${\nu }_1\to 0$, $F_{v-el}\left(t\right)$ behaves as an elastic force, while in the limit $k_1\to \infty$, ${\nu }_1\to \infty$, ${\eta }_1=k_1/{\eta }_1=const$ it corresponds to the viscous Stokes friction with the friction coefficient ${\eta }_1$. This is how Maxwell derived the phenomenon of viscosity from the phenomenon of elasticity, i.e., by letting the elastic stress to relax in time [87]. Complex viscoelastic liquids and gels are, however, characterized by a power law decaying memory kernel, ${\eta }_{\mathrm{mem}}\left(t\right)={\eta }_{\alpha }{t}^{-\alpha }/\Gamma (1-\alpha )$, with $0<\alpha <1$, as it has first been established by Gemant [87,100,125], and which now presents a common model. In this particular case, $F_{v-el}\left(t\right)$ can be abbreviated as $F_{v-el}\left(t\right)={\eta }_{\alpha }{d}^{\alpha }\varphi \left(t\right)/d{t}^{\alpha }$, which just defines the fractional Caputo derivative of the order $\alpha$. Hence, ${\eta }_{\alpha }$ is customarily named the fractional friction coefficient, and GLE in this particular case is named fractional Langevin equation or FLE. This is, of course, an idealization. In reality, there are always two memory cutoffs present. A large time memory cutoff ${\tau }_h=1/{\nu }_l$ defines the slowest Maxwellian relaxation mode of the environment, and with ${\eta }_{\alpha }\to {\eta }_{\alpha }\exp (-{\nu }_{h}t)$, an effective friction coefficient ${\eta }_{\mathrm{eff}}={\int }_0^{\infty }{\eta }_{\mathrm{mem}}\left(t\right)dt={\eta }_{\alpha }{\tau }_h^{1-\alpha }$ can be introduced. Notice that ${\eta }_{\mathrm{eff}}$ characterizes diffusion on the time scale $t\gg {\tau }_h$. However, ${\tau }_h$ can be well in the range of minutes and even hours. It depends on the system considered. We assume it to be in the range of seconds for our nanosensor. As long as $t<{\tau }_h$, e.g., for the duration of sojourn times of our sensor in the metastable states, it is ${\eta }_{\alpha }$ that determines the stochastic dynamics and not ${\eta }_{\mathrm{eff}}$, which can be effectively infinite, with ${\tau }_h\to \infty$. This is the reason why thinking in terms of some ${\eta }_{\mathrm{eff}}$ can be very misleading for viscoelastic media. This is a macroscopic type approximation, which can fail completely on micro- and nano-scales. With this reservation we use it because far too many researchers continue to think in terms of some ${\eta }_{\mathrm{eff}}$. In the simulations presented below we fixed $\alpha =0.5$ (one of common experimental values for cytosol [74]) and ${\tau }_h={10}^4$ (or about 1.7 s, when ${\tau }_{sc}=0.17$ ms). ${\eta }_{\alpha }$ will be fixed to two values by choosing ${\eta }_{\mathrm{eff}}$ = 100 ${\eta }_0$ and ${\eta }_{\mathrm{eff}}$ = 1000 ${\eta }_0$. In dimensionless units used in numerics, the former (intermediate) fractional friction is about ${\eta }_{\alpha }\approx {\eta }_0$, whereas the latter one (strong) is ${\eta }_{\alpha }\approx 10$ ${\eta }_0$. For intermediate ${\eta }_{\alpha }$, the relaxation within a potential well is mostly exponential with a heavy power law tail, while for the large ${\eta }_{\alpha }$, it is initially stretched exponential and then changes into a power-law decay (Mittag-Leffler relaxation function) [61]. This latter one corresponds to dielectric Cole–Cole response [126] which is typical for biological media [4]. Furthermore, on physical grounds also a short time cutoff ${\tau }_l=1/{\nu }_h$ is always necessarily present. It ensures that spectral density of the noise ${\xi }_{\mathrm{mem}}\left(t\right)$ does not contain frequencies much above ${\nu }_h$. This physically corrects the approximation of continuum medium where such a cutoff is absent because this approximation neglects the atomistic nature of any real condensed medium. In our numerics, we take it to be ${\nu }_h={\nu }_0={10}^4$, and the numerical method is based on approximating the memory kernel by a sum of exponentials, and using a Markovian embedding (9) of the GLE dynamics in Equation (6), see in Methods. This allows for a numerically highly accurate approach, with a well controlled accuracy [87,100]. Typical trajectories for strong fractional friction are shown in Figure 6. For a weak or intermediate fractional friction, they look more like the ones in Figure 4. One striking feature is immediately seen in Figure 6. This is a highly bursting character of sojourns in the open state, where huge many very short excursions happen to the closed state during a long sojourn in the open state. It visually signals a truly non-exponential kinetics.

3.2.1. Intermediate fractional friction, ${\eta }_{\alpha }\approx {\eta }_0$

The first profound influence of fractional viscoelastic friction on the statistics of the residence time distributions can be reveal in Figure 7. Namely, the distribution of closed times becomes stretched exponential with ${\beta }_c\approx 0.92$ in part a, and ${\beta }_c\approx 0.82$ in part b. However, the distribution of open times remains almost exponential, and the mean times in the states and the mean opening probability remain only weakly affected, compare with Figure 5a,b. This is especially striking because we can further arbitrarily increase ${\eta }_{\mathrm{eff}}$, while keeping the same ${\eta }_{\alpha }$. This is easy to do in our numerics just by further increasing the number N of exponentials in Equation (7) and, correspondingly, the Markovian embedding dimension in Equation (9). The results will not be changed because the essential kinetics in Figure 7 is on the time scale which is already smaller than the cutoff time ${\tau }_h$, that further increases with N. This feature must be shocking for all those who continue thinking in terms of ${\eta }_{\mathrm{eff}}$, rather than fractional friction featuring such complex media as cytosol. All in all, for such an intermediate ${\eta }_{\alpha }$, our sensor operates as fast as in water. This is a good news.

3.2.2. Strong Fractional Friction, ${\eta }_{\alpha }\approx$ 10 ${\eta }_0$

For a strong fractional friction, both closed and open times have Weibull distribution with smaller values of the stretched exponential parameter $\beta$ (see Figure 8). Here, the influence of viscoelastic effects is very strong. Initially, for small time, $\beta$ can exceed one, see for closed times (red fitting curve) in part b. However, the mean opening probabilities are also almost unaffected, as compared with the Markovian case, despite the mean residence times increase. This increase is not very strong, by a factor of less than three only. Our sensor remains very fast.

4. Discussion

Design of a magnetosensitive ion channel complex using a single biomagnetite nanoparticle as a magnetic field sensor makes it more realistic. In essence, the model studied in this paper is a variant of the model in Ref. [61]. The differences are in detail. However, these details do matter. First, a single prolonged nanoparticle is used instead of a rod made of $5÷7$ such nanoparticles, and, second, a shorter linker is used to restrict the orientational motion of sensor when it fluctuates between two metastable states in response to a change of the magnetic field orientation. In our present work, it is just about ${30}^{\circ }$, whereas, in [61], it is about ${150}^{\circ }$. This makes the present variant much faster and more realistic in view of possible sterical restrictions within the cell. This comes, however, at a price: the magnetic moment of the sensor must be larger to ensure its magnetic energy to be about $10k_{B}T$ in the magnetic fields of Earth vs. $3÷4k_{B}T$ in Ref. [61]. This is, however, not a problem because such sufficiently large nanoparticles are met both in certain bacteria, and in the human brain, even if they are certainly less common that ones assumed in [61]. Furthermore, our analysis shows that such a sensor would be operating very fast even in viscoelastic cytosol feeling an effective viscosity 1000× larger than one of water. More precisely, this effective viscosity can be even formally infinite, like in solid, because this is not an effective macroscopic friction defined on a very large time scale which determines the stochastic switching dynamics. In fact, the microscopic fractional friction does matter here and it is crucially relevant. If it is large and dominates the dynamics, as in the studied case of ${\eta }_{\mathrm{eff}}$ = 1000 ${\eta }_0$, the “off–on” two-state dynamics becomes very bursting. It is not characterized by nearly exponential distributions of the residence times in two states anymore, but rather by profoundly non-exponential stretched exponential distributions. Such distributions are indeed measured in several biological ion channels, and our theory tentatively explains their principal origin as one rooted in the viscoelasticity of environment. All this clearly correlates with the dielectric Cole–Cole response in such media (Mittag-Leffler viscoelastic relaxation). It is important that the mean residence times (MRTs) as well as all higher moments remain finite; moreover, MRTs were increased in our model study by a factor of less than three only, as compared to ones in water. Arguably, it can be unbelievable and embarrassing for those who continue to think in terms of an ${\eta }_0\to {\eta }_{\mathrm{eff}}$ renormalization within a Markovian dynamics. However, this is a result of the proper treatment of non-Markovian effects, and it presents a very good news with respect to feasibility of such magnetosensitive complexes in living systems.

If such magnetosensitive ion channels do exist, why then they were not found until now? The situation here can be similar to ion channels associated with cilia in hair cells. The existence of those channels is widely assumed and is taken nowadays by many almost as a real fact. However, they were not identified until now as biological protein structures, such as many other well-known ion channels, despite numerous efforts. It is very difficult to identify them because each cilia is assumed to be connected by elastic protein linkers to the gates of many such channels, and it is difficult to confirm this hypothesis experimentally. The ion channels of cilia in hair cells remain elusive, even if in their existence practically nobody doubts. The hypothesis of magnetosensitive ion channel complexes is much less studied. However, it is a reasonable one and it should attract more attention in the future.

5. Methods

The power law memory kernel is approximated between two cutoffs by a sum of exponentials

$$\begin{array}{c}\hfill {\eta }_{\mathrm{mem}}\left(t\right)=\sum _{i=1}^N{k}_i\exp (-{\nu }_{i}t),\end{array}$$

(7)

with a fractal scaling of relaxation rates ${\nu }_i={\nu }_0/b^{i-1}$ and weights $k_i\propto {\nu }_i^{\alpha }$. Namely,

$$\begin{array}{c}\hfill k_i={\nu }_0{\eta }_{\mathrm{eff}}\frac{b^{1-\alpha }-1}{b^{(i-1)\alpha }[b^{N(1-\alpha )}-1]}.\end{array}$$

(8)

Here, ${\nu }_0={\nu }_h=1/{\tau }_l$ is the largest viscoelastic rate of the environment. Using scaling or dilation parameter $b=10$ allows to achieve 4% accuracy of power law approximation between two cutoffs for $\alpha =0.5$ and a sufficiently large N [87,100]. In our numerics, we use ${\nu }_0={10}^4$ and $N=9$, ${\tau }_h=b^{N-1}/{\nu }_0={10}^4$. The fractional friction coefficient is ${\eta }_{\alpha }={\eta }_{\mathrm{eff}}{\tau }_h^{\alpha -1}/g_{\alpha }$, with an inverse proportionality coefficient $g_{\alpha }$, which slightly differs from unity. For $\alpha =0.5$, $b=10$, and $N>5$, $g_{\alpha }\approx 1.07$. This allows for a Markovian embedding

$$\begin{array}{ccc}\hfill {\eta }_0\dot{\varphi }& =& f\left(\varphi \right)-\sum _{i=1}^N{k}_i(\varphi -y_i)+{\xi }_0\left(t\right),\hfill \\ \hfill {\eta }_i\dot{y_i}& =& k_i(\varphi -y_i)+{\xi }_i\left(t\right),\hfill \end{array}$$

(9)

of GLE dynamics (6) in a space of $N+1$ dimension. Here, $y_i$ are non-dimensional linear auxiliary variables, ${\eta }_i=k_i/{\nu }_i$, and ${\xi }_i\left(t\right)$ are mutually independent auxiliary uncorrelated white Gaussian noises,

$$\begin{array}{ccc}\hfill \langle {\xi }_i\left(t\right){\xi }_j\left(t^{\prime }\right)\rangle & =& 2{\delta }_{ij}{k}_{B}T{\eta }_j\delta (t-t^{\prime }).\hfill \end{array}$$

(10)

They are also uncorrelated with ${\xi }_0\left(t\right)$. The initial angles $y_i\left(0\right)$ are sampled from Gaussian distributions centered around $\varphi \left(0\right)$, $\langle y_i\left(0\right)\rangle =\varphi \left(0\right)$ with variances $\langle {[y_i\left(0\right)-\varphi \left(0\right)]}^2\rangle =k_{B}T/k_i$, in order to have a complete equivalence with the corresponding GLE in Equations (6) and (7) in the ensemble sense [87]. The system of stochastic differential Equations (9) is propagated in time (in units of ${\tau }_{sc}={\eta }_0/U_0$) using stochastic Heun method [120]. The time step was fixed to $\delta t=2·{10}^{-6}$. $T=T_r=297$ K, while $U_0=10k_B{T}_r$.

6. Conclusions

To conclude, in this paper, we studied a model of magnetosensitive ion channel complexes in more realistic detail. Our study underpins theoretically a possible existence of such complexes which should attract more attention of researchers trying to resolve the puzzle of magnetosensitivity of many animals to environmental magnetic fields. The author is convinced that ion channels are involved in magnetosensing, even if the concrete structures involving them can be different. More research in this direction is required and welcome.