11institutetext:  Institute of Theoretical Physics, University of Wrocław,
PL–50204 Wrocław, Poland ExtreMe Matter Institute EMMI, GSI, D-64291
Darmstadt, Germany

# Probing QCD chiral cross over transition in heavy ion collisions with
fluctuations

Krzysztof Redlich<EMAIL_ADDRESS>
###### Abstract

We argue that by measuring higher moments of the net proton number
fluctuations in heavy ion collisions (HIC) one can probe the QCD chiral cross
over transition experimentally. We discuss the properties of fluctuations of
the net baryon number in the vicinity of the chiral crossover transition
within the Polyakov loop extended quark-meson model at finite temperature and
baryon density. The calculation includes non-perturbative dynamics implemented
within the functional renormalization group approach. We find a clear signal
for the chiral crossover transition in the fluctuations of the net baryon
number. We address our theoretical findings to experimental data of STAR
Collaboration on energy and centrality dependence of the net proton number
fluctuations and their probability distributions in HIC.

QCD phase transition and phase diagram  Chiral symmetry and charge
fluctuations  Heavy ion collisions

††articletype:  Communication

##  I QCD phase diagram

Although, the important question addressed in QCD on the existence of a true
2nd order phase transition at finite chemical potential (critical point) has
not been answered yet, nevertheless, there is recently an essential progress
in the quantitative description of the QCD phase diagram. The lattice QCD
(LQCD) provided a final value for the chiral cross over transition temperature
at vanishing chemical potential l1 . The LQCD has also provided arguments that
at small $\mu\simeq 0$ the chiral cross over line is the pseudo critical line
of the 2nd order chiral phase transition belonging to the universality class
of 3-dimensional, O(4) symmetric spin models l2 . Based on the universality
arguments the LQCD has also provided the curvature of the chiral transition
line, $T_{c}(\mu_{c})$ l2 . These results confirm that, at least at small
values of the baryon chemical potential $\mu$, the chiral cross over
transition appears in the near vicinity to the chemical freezeout line f1
obtained from the analysis of particle yields measured in HIC.

The numerical coincidence of thermal parameters for the chiral transition and
the freezeout conditions indicates that the hadron resonance gas partition
function, which describes chemical equilibration of particle yields in HIC,
should describe also the QCD thermodynamics up to a near vicinity to the
transition to a quark gluon plasma phase. Indeed, the equation of state
calculated on the lattice as well as other thermodynamical observables in the
hadronic phase, were shown to be very well quantified by the hadron resonance
gas (HRG) partition function h1 .

Already at vanishing chemical potential, i.e. under conditions realized in the
high energy runs at RHIC or LHC, the question arises to what extent a refined
analysis of freeze-out conditions can establish the existence of a chiral
phase transition. In this lecture we will argue that even at $\mu_{B}/T\simeq
0$ the net baryon number fluctuations and their higher moments can be used to
identify the chiral cross over transition experimentally o1 ; o2 ; o3 ; o4 .

Figure 1:  The temperature dependence of kurtosis
$R_{4,2}:=9\chi_{4}^{B}$/$\chi_{2}^{B}$ and the higher order,
$\chi_{6}^{B}$/$\chi_{2}^{B}$ and $\chi_{8}^{B}/\chi_{2}^{B}$, ratios of
cumulants for different $\mu_{q}/T$ calculated in the PQM model within the FRG
approach o2 ; o4 . The $T_{pc}$ is the pseudo-critical temperature obtained in
the model at the physical pion mass.

##  II Charge fluctuations and the chiral cross over transition

Due to remnants of O(4) criticality related with the chiral phase transition
observed in LQCD, the free energy ($f$) near the chiral phase transition
temperature $T_{c}$ may be represented in terms of singular ($f_{s}$) and
regular contributions ($f_{r}$) as

$f(T,\mu_{q},m_{q})=f_{s}(T,\mu_{q},m_{q})+f_{r}(T,\mu_{q},m_{q})\;,$ (1)

where in addition to the temperature $T$ we also introduced explicit
dependence on the light quark chemical potential, $\mu_{q}=\mu_{B}/3$, and the
(degenerate) light quark masses $m_{q}\equiv m_{u}=m_{d}$. The singular part
of the free energy may be written as s1

$f_{s}(T,\mu_{q},h)=h^{1+1/\delta}f_{s}(z)\ ,\ z\equiv t/h^{1/\beta\delta}$
(2)

with $\beta,\ \delta$ are critical exponents of the 3-dimensional, $O(4)$
universality class and
$t\equiv\frac{1}{t_{0}}\left(\frac{T-T_{c}}{T_{c}}+\kappa_{q}\left(\frac{\mu_{q}}{T}\right)^{2}\right)$,
$h\equiv\frac{1}{h_{0}}\frac{m_{q}}{T_{c}}$. Here $T_{c}$ is the phase
transition temperature in the chiral limit and $t_{0}$, $h_{0}$ are non-
universal scale parameters. The proportionality constant $\kappa_{q}\simeq
0.06$, has recently been determined from a scaling analysis in (2+1)-flavor
QCD l2 . The scaling function $f_{s}$ and its derivatives have recently been
calculated using high precision Monte Carlo simulations of the 3-dimensional
O(4) spin model e1 .

We want to focus here on properties of moments of net baryon number
fluctuations, which are obtained from Eq. 1 by taking derivatives with respect
to $\hat{\mu}_{B}=\mu_{B}/T$,

$\chi_{n}^{B}=-\frac{1}{3^{n}}\frac{\partial^{n}f/T^{4}}{\partial\hat{\mu}_{q}^{n}}=-\frac{1}{3^{n}}\frac{\partial^{n}f_{r}/T^{4}}{\partial\hat{\mu}_{q}^{n}}-\frac{1}{3^{n}}\frac{\partial^{n}f_{s}/T^{4}}{\partial\hat{\mu}_{q}^{n}}=\chi_{n,r}^{B}+\chi_{n,s}^{B}\;.$
(3)

In the hadronic phase and away from transition temperature the regular part
$\chi_{n,r}^{B}$ should be well described by the hadron resonance gas
partition function which will be a reference for critical fluctuations coming
from the singular part $\chi_{n,s}^{B}$. Thus, any deviations from the regular
i.e. hadronic gas contribution could be an indication of criticality due to
remnants of the chiral O(4) transition.

Higher order moments will become increasingly sensitive to the singular part
of the free energy. From Eq. 2 it is apparent that these moments show a strong
quark mass dependence in the vicinity of the critical temperature,

$\chi_{n,s}^{B}\sim\begin{cases}-(2\kappa_{q})^{n/2}h^{(2-\alpha-n/2)/\beta\delta}f_{s}^{(n/2)}(z)&,\
{\rm for}\ \mu_{q}/T=0,\ {\rm and}\ n\ {\rm even}\\\
-(2\kappa_{q})^{n}\left(\frac{\mu_{q}}{T}\right)^{n}h^{(2-\alpha-n)/\beta\delta}f_{s}^{(n)}(z)&,\
{\rm for}\ \mu_{q}/T>0\end{cases}\ .$ (4)

where we used $2-\alpha=\beta\delta(1+1/\delta)$. As $\alpha=-0.2131(34)$ is
negative in the 3-dimensional, $O(4)$ universality class, the 4th order
moments of the net baryon number fluctuations do not diverge yet in the chiral
limit at the chiral transition temperature, $z=0$. The first divergent moment
is obtained for $n=6$ if $\mu_{q}/T=0$ and for $n=3$ if $\mu_{q}/T>0$.

Figure 2:  Left-hand figure: the ratio of quadratic fluctuations and mean net
baryon number ($\sigma^{2}/M$), cubic to quadratic ($S\sigma$) and quartic to
quadratic ($\kappa\sigma^{2}$) baryon number fluctuations calculated in the
HRG model on the freeze-out curve o1  and compared to results obtained by the
STAR Collaboration star . The dashed curves show the approximate
$\tanh(\mu_{B}/T)$ result for $\kappa\sigma^{2}$ and $S\sigma$, respectively.
Middle figure: The probability distributions, uncorrected for event-by-event
counting efficiency, for the net proton number for different centralities
taken by STAR Collaboration in Au-Au collisions at $\sqrt{s_{NN}}=200$ GeV
star . The lines are the Skellam distributions calculated within HRG model o3
. Right-hand figure: Mean (M), variance ($\sigma$), skewness (S) and kurtosis
($\kappa$) of the net proton number calculated from the probability
distributions shown in the middle figure.

In the hadron resonance gas (HRG) the regular part of the fluctuations
$\chi_{n,r}^{B}$ can be directly calculated from the thermodynamic pressure
following Eq. 3. The HRG is a mixture of ideal gases of all particles and
resonances, consequently the thermodynamic pressure exhibits a factorization
of $T$ and $\mu_{B}/T$ dependence. Under the Boltzmann approximation the
pressure in the HRG, $P^{HRG}(T,\mu_{B})\simeq f(T)\cosh(\mu_{B}/T)$, where
$f(T)$ contains contributions from all baryons and baryonic resonances. With
such $P^{HRG}$ there are particular properties of ratios of cumulants,

${{\chi_{2n,r}^{B}}\over{\chi_{2,r}^{B}}}|_{HRG}=1~{}~{},~{}~{}{{\chi_{(2n+1),r}^{B}}\over{\chi_{1,r}^{B}}}|_{HRG}=1~{}~{},~{}~{}{{\chi_{(2n+1),r}^{B}}\over{\chi_{2,r}^{B}}}|_{HRG}\simeq\tanh(\mu_{B}/T)~{}~{},~{}~{}{{\chi_{2n,r}^{B}}\over{\chi_{1,r}^{B}}}|_{HRG}\simeq\coth(\mu_{B}/T),$
(5)

which are independent of the number of baryons, their masses, degeneracy
factors or decay widths.

The above structure of cumulants ratios observed in the HRG will be modified
if the singular part is included in Eq. 3. The $\chi_{n,s}^{B}$ are
increasingly sensitive to the order of cumulants. Thus, ratios of cumulants
with different $n$ should be strongly varying functions of $T$ and $\mu_{B}$
when approaching a chiral cross over transition. Consequently, the observed
deviation from that expected in Eq. 5 could be considered as a signature of
the singular part contribution to the overall fluctuations thus, also of the
chiral cross over transition.

The generic structure of ratios of different cumulants near the chiral cross
over transition is shown in Fig. 1. These ratios were calculated in the
Polyakov loop extended quark-meson (PQM) model at the physical pion mass,
applying the functional renormalization group (FRG) method. In FRG approach
one includes quantum and thermal fluctuations which are needed to preserve the
universal scaling behavior of physical quantities expected in the O(4)
universality class. In the low temperature phase the ratios
${{\chi_{2n,r}^{B}}/{\chi_{2,r}^{B}}}=1$ as expected in the hadron resonance
gas from Eq. 5. For $n=2$, the kurtosis
$R_{4,2}=9{{\chi_{4,r}^{B}}/{\chi_{2,r}^{B}}}$ is not affected by the chiral
critical dynamics since there is no contribution of the singular part to the
first four moments as seen in Eq. 4. The observed in Fig. 1 drop in $R_{4,2}$
is due to ”statistical confinement” property of the PQM model stok . Such
behavior of kurtosis was first observed in LQCD calculations and was
interpreted as being a signature of deconfinement in QCD s1 .

The large deviations of $(R_{4,2}/9)$ from unity in hadronic phase, which are
increasing with $\mu/T$, are due to a singular contribution to the
${\chi_{4}^{B}}$ and ${\chi_{2}^{B}}$ ratio (see Eqs. 3 and 4). With
increasing order of cumulants and the value of the chemical potential their
ratios are dominated by the singular part already deeply below the chiral
cross over transition temperature. Such behavior could be observed in HIC if
freezeout appears near the chiral transition.

Recently, the first data on charge fluctuations and higher order cumulants,
identified through the net-proton fluctuations, were obtained by the STAR
Collaboration in Au- Au collisions at several collision energies star . To
explore possible signs of chiral criticality and a cross over transition, the
STAR data on the first four moments are compared in Fig. 2 to HRG o1 ; o3
following Eq. 5. Different ratios of cumulants of the net proton number can be
directly connected with measured mean (M), variance $(\sigma)$, skewness (S)
and kurtosis $(\kappa)$ o1 .

The basic properties of measured fluctuations and ratios of cumulants are
consistent with that expectated in the HRG model o1 ; o3 . This indicates,
that moments of the net proton number are of thermal origin with respect to
the grand canonical ensemble and that they freezeout along the same chemical
freezout line as particle yields, close to the chiral cross over line.
However, already such first comparison of the HRG model with STAR data reveals
that there are deviations o1 ; o3 ; cp1 . This is seen in Fig. 2 on the level
of different ratios of cumulants as well as by comparing directly the measured
probability distributions with Skellam distribution expected in the HRG for
the net proton number o3 . The HRG model, as seen in Fig. 2, results in a
broader distribution than observed in data, particularly at the most central
collisions. This implies deviations of data on different moments from the HRG
model which are increasing with centrality. Shrinking of measured widths of
the net proton number distribution relative to the HRG results, seen in Fig.
2, is to be already expected due to deconfinement properties of QCD o3 . We
have to stress, however, that the experimental net-proton distributions in
Fig. 2 are not corrected for event-by-event proton/anti-proton counting
efficiency. While the HRG lines in this figure used efficiency corrected mean
multiplicities from the same experiment. If uncorrected data are used in the
Skellam distribution, then the agreement of the HRG model and measure
probability distributions is found to be much better cp2 .

From the model calculations in Fig. 1, as well as from Eqs. 3 and 4, it is
clear that contributions of the O(4) singular part to fluctuations increase
with the order of cummulants. Preliminary data of STAR Collaboration on
$\chi_{6}^{B}/\chi_{2}^{B}$ ratio, taken in Au-Au collisions at the top RHIC
energy, show a strong deviation of this ratio from the HRG model expectations
with a rather moderate deviations of the lower order cumulant ratios.
Recently, the STAR Collaboration has also observed that deviations in central
Au-Au collisions of $\chi_{4}^{B}/\chi_{2}^{B}$ ratio from the HRG is non-
monotonic in energy. Such behaviors of ratios of cumulants are to be expected
due to remnants of criticality related with a chiral cross over transition. It
is a further challenge to understand and quantify the observed energy and
centrality dependence of these deviations from the HRG results.

I am grateful for fruitful collaboration and discussions with P. Braun-
Munzinger, B. Briman, F. Karsch and V. Skokov. Partial support by MEN (Poland)
and EMMI (Germany) and discussion with Nu Xu is also acknowledged.

## References

  * (1) Y. Aoki, et al., JHEP 0906, 088 (2009). A. Bazavov et. al. (HotQCD collaboration), arXiv:1111.1710 [hep-lat].
  * (2) S. Ejiri, et al., Phys. Rev. D 80, 094505 (2009). O. Kaczmarek, et al., Phys. Rev. D 83, 014504 (2011)
  * (3) P. Braun-Munzinger, et al., in Quark-Gluon Plasma 3, Eds. R.C. Hwa and X.N. Wang, (World Scientific Publ., 2004). J. Cleymans, K. Redlich, Phys. Rev. C 60, 054908 (1999); Phys. Rev. Lett. 81, 5284 (1998).
  * (4) F. Karsch, et al., Eur. Phys. J. C 29, 549 (2003). C. Ratti, Nucl. Phys. A 855, 253 (2011).
  * (5) F. Karsch, K. Redlich, Phys. Lett. B 695, 136 (2011).
  * (6) B. Friman, et al., Eur. Phys. J. C 71, 1694 (2011).
  * (7) P. Braun-Munzinger, et al., Phys. Rev. C 84, 064911 (2011); arXiv:1111.5063 [hep-ph].
  * (8) V. Skokov, et al., Phys. Rev. C 83, 054904 (2011). J. Phys. G 38,124102 (2011). Phys. Lett. B 708, 179 (12).
  * (9) See also: Lizhu Chen, these proceedings.
  * (10) See also: Xiaofeng Luo and B. Mohanty, these proceedings.
  * (11) S. Ejiri, et al., Phys. Lett. B 633, 275 (2006).
  * (12) J. Engels, F. Karsch, arXiv:1105.0584 [hep-lat].
  * (13) B. Stokic, et al., Phys. Lett. B 673, 192 (2009).
  * (14) M. M. Aggarwal et al. [STAR Collaboration], Phys. Rev. Lett. 105, 022302 (2010). X. Luo, et al., [for the STAR Collaboration], arXiv:1106.2926v1.


# Biochemical pathways of forward/backward steps of motor protein kinesin-1

Beibei Shen and Yunxin Zhang<EMAIL_ADDRESS>Shanghai Key Laboratory for
Contemporary Applied Mathematics, School of Mathematical Sciences, Fudan
University, Shanghai 200433, China.

###### Abstract

Kinesin-1 is an ATP-driven, two-headed motor protein that transports
intracellular cargoes along microtubule. Based on recent experimental
observations, we formulate a mechanochemical model for it, in which
forward/backward/futile cycle of kinesin-1 can be realized through multiple
biochemical pathways. Our results show that both forward and futile cycles
consist of two biochemical pathways, while backward cycle may be realized
through six possible pathways. Backward motion of kinesin-1 is mainly through
backward sliding along microtubule when it is in semi-detach state, and
usually coupled with ATP hydrolysis. Under low external load, about 80% of ATP
is wasted by kinesin-1, and this proportion of waste is almost independent of
ATP concentration. At high ATP concentration or under high external load, both
heads of kinesin are always in ATP or ADP$\cdot$Pi binding state and bind to
microtubule tightly. But at low ATP concentration and low load, kinesin mainly
stays at the one head bound state. In contrast to mean run length, the mean
run time of kinesin-1 along microtubule decreases with ATP concentration.
Unless external load is near the stall force, the motion of kinesin-1 is
almost deterministic, not as predicted by usual Browian ratchet models.

forward/backward/futile biochemical pathways, mean run time/length,
conventional kinesin.

##  I Introduction

Kinesin-1 plays a central role in the intracellular transport of various
vesicles and organelles Howard (2001); Vale (2003); Schliwa (2003); Kolomeisky
(2015). Each kinesin dimer consists of two motor domains (heads), which are
connected to a coiled-coil stalk through a $\sim$14-amino-acid-long sequence
known as the neck linker (NL) Case et al. (2000); Al-Bassam et al. (2003);
Tomishige et al. (2006); Kalchishkova and Böhm (2008); Hariharan and Hancock
(2009). Kinesin proceeds unidirectionally toward the plus ends of microtubule
(MT) by an asymmetric hand-over-hand fashion, hydrolyzing an ATP molecule for
each 8 nm step Asbury et al. (2003); Yildiz et al. (2004); Schnitzer and Block
(1997); Coy et al. (1999); Svoboda et al. (1993); Higuchi et al. (2004);
Guydosh and Block (2009). Motion of kinesin is significantly processive Carter
and Cross (2005); Toprak et al. (2009), taking about 100 steps before
detaching from the MT Howard (2001); Clancy et al. (2011a); Ariga et al.
(2018). At low loads, kinesin-1 always steps forward to the plus ends of MT,
while at high loads, it can also result in backward steps Svoboda and Block
(1994); Visscher et al. (1999a); Block et al. (2003).

In previous studies, kinds of theoretical models have been presented to
understand the mechanism of kinesin motion; see Parrondo and Cisneros (2002);
Kolomeisky and Fisher (2007); Chowdhury (2013); Mugnai et al. (2020) for
details. Generally, the motion of kinesin can be regarded as a Markov process
with multiple biophysical and biochemical states, and in steady-state,
meaningful quantities can be obtained by using Fokker-Planck equation, or
master equation, or Langevin dynamics, or any other complicated hybrid models
Kampen (2007); Gardiner (2010). Among most of these models, a backward step of
kinesin is usually thought to occur by directional reversal of a forward step
and therefore happens with ATP synthesis.

Generally, both forward and backward steps of kinesin may be accomplished
through different biochemical pathways, and backward motion of kinesin may be
resulted simply from pure biophysical slips along MT Toleikis et al. (2020).
According to a mechanochemical model, our study shows that there are two
biochemical pathways for forward step. At low ATP concentration and low
external load, kinesin will spend more time waiting for ATP, and its MT-bound
head might be in the nucleotide-free state. However, at high ATP concentration
or high external load, all MT-bound heads are usually in ATP or ADP$\cdot$Pi
binding state, and ATP binding to the front head of kinesin is earlier than
the release of Pi and the detachment of the trailing head from MT. There are
altogether six possible pathways for backward motion, unless at some extreme
conditions the backward motion of kinesin is through backward slips along MT
when it is in semi-detach state, and usually coupled with ATP hydrolysis.
Directional reversal of forwarding pathways is only non-negligible when the
ATP molecule is scarce while the external load is very high. One surprising
finding is that about 80% of ATP is wasted by kinesin, and this proportion of
waste is almost independent of ATP concentration.

Figure 1:  The mechanochemical model of kinesin-1. (a) The depiction of
detailed configurations of the two heads of kinesin-1. There are six possible
configurations at each binding site of MT. (b) Mechanochemical model used in
this study to describe the periodic motion of kinesin-1, where $n_{j}$ denotes
the state that kinesin is at position $n$ and in configuration $j$, with
$j=1,2,\cdots,6$.

##  II The Mechanochemical model of kinesin-1

Based on experimental observations in Toleikis et al. (2020), this study uses
the mechanochemical model as depicted in Fig. 1. The biochemical cycle of
kinesin is assumed to begin with the one-head-bound, ATP-waiting state, with
the nucleotide-free front head is tightly bound to the MT while the ADP-bound
rear head unbound Block (2007). ATP binding induces NL docking and ATP
hydrolysis completes NL docking Case et al. (2000); Tomishige et al. (2006);
Kalchishkova and Böhm (2008) (see Fig. 1(a) 1$\to$2). In either ATP or
ADP$\cdot$Pi bound states, the head binds tightly to the MT, so these two
states are assumed to be the same and denoted as state 2 in Fig. 1. Due to NL
docking, the unbound head of kinesin swings forward and becomes the front
head. ADP release facilitates front head binding to MT, leading kinesin
changes from state 2 to state 3, generating a two-heads-bound state. At low
ATP concentrations, the release of phosphate ion Pi from the rear head usually
occurs before an ATP molecule reaches the newly bound front head. Then
dissociation of the ADP-bound rear head from the MT leads kinesin to return to
biochemical state 1 again. However, if the ATP concentration reaches
saturating, after ADP releases from the front head, another ATP molecule might
bind to the nucleotide-free head immediately, leading kinesin to enter the
biochemical state 4. Next, the release of Pi from the rear head will lead
kinesin to return to state 2. In either case, 1$\to$2$\to$3$\to$1 or
2$\to$3$\to$4$\to$2, one biochemical cycle of kinesin is completed once the
ADP bound rear head dissociates from the MT. With completion of one
biochemical cycle, kineisn will move forward about 8 nm along the MT, and is
ready to start a new mechanochemical cycle again Schnitzer and Block (1997);
Coy et al. (1999); Higuchi et al. (2004).

Recent experiments discovered that the dwell time for forwarding step is
shorter than that for backward step, and dwell times for backward step and
detachment are almost the same Toleikis et al. (2020). Here dwell time is the
waiting time preceding each step, consisting of the time spent waiting for ATP
to bind and the time is taken to process ATP and complete the coupled
mechanical step. This implies that backward step originates from a different
and later state in the biochemical cycle than the one that generates forward
steps, and backward step and detachment might emanate from the same state.
Therefore, our model assumes that forward step can only occur before Pi
release from the ADP$\cdot$Pi bound head (see state 2 in Fig. 1), whereas
backward step/slip and detachment can only occur after Pi release. The release
of Pi leads to a weakly bound ADP state (state 5). With increasing external
load, forward step of motor becomes difficult to complete before Pi release,
while backward step may occur from state 5. If the rear head releases its ADP
before it is semi-detached from the MT, kinesin in state 5 can return to state
1. Evidently, one ATP is hydrolyzed during cycle 1$\to$2$\to$5$\to$1 but
without output of mechanical work, which we called futile cycle in this study.
As concluded in Clancy et al. (2011b), the inclusion of futile cycles is
experimentally meaningful and theoretically essential to make the model fit
the experimental data well.

Semi-detach state 6 is the weakest bounding state in the biochemical cycle of
kinesin, and any backward slipping necessarily requires passage through this
state. From state 6, kinesin may dissociate completely from MT, and enter into
the detachment state, denoted as Off in Fig. 1. More importantly, as described
in Toleikis et al. (2020), kinesin in state 6 may also slip backward along the
MT. For simplicity, our model assumes that, unless it is completely detached
from the MT, kinesin always slips backward $d=8$ nm along the MT in each
slipping process, then releases the ADP from its semi-detached head and
returns to state 1. This assumption can be relaxed to allow kinesin to slip
backward any distance of integer multiple of 8 nm. But results of fitting to
experimental data show that the load free rate $k_{61}^{0}$ is very small, so
the possibility of backward sliding is small, and this simplification is
reasonable, see Tab. 1. Meanwhile, kinesin in semi-detach state 6 may also
slide forward along the MT. But, since the probability of kinesin in state 6
is very small (see Fig. 3), compared with the total forward probability flux,
the flux produced by forwarding sliding is negligible.

The model used in this study is similar to the one constructed in Toleikis et
al. (2020). One of the key features distinguishing it from earlier ones is the
inclusion of semi-detach state 6, from which the motor can slip backward to
state 1 or dissociate from MT completely. Generally, biochemical cycles of
kinesin in this model can be classified into three categories, (i) forward
cycle coupled with one forward mechanical step, (ii) backward cycle coupled
with one backward mechanical step, and (iii) futile cycle with one ATP
hydrolyzed but without change of position.

As illustrated in Fig. 1(a), biochemical transitions 2$\rightleftharpoons$3,
3$\rightleftharpoons$1, and 4$\rightleftharpoons$2 are coupled with head
attachment/detachment to/from MT binding site, which may result in the change
of mass center of kinesin. So rates of these transitions, i.e., $k_{23}$,
$k_{32}$, $k_{31}$, $k_{13}$, $k_{42}$, $k_{24}$, are assumed to be external
load $F$ dependent. Similar as in Fisher and Kolomeisky (2001); Zhang (2012),
this study assumes that

$\displaystyle k_{23}(F)=k_{23}^{0}e^{\frac{-\delta_{23}Fd}{k_{B}T}},\quad
k_{32}(F)=k_{32}^{0}e^{\frac{\delta_{32}Fd}{k_{B}T}},$ (1) $\displaystyle
k_{31}(F)=k_{31}^{0}e^{\frac{-\delta_{31}Fd}{k_{B}T}},\quad
k_{13}(F)=k_{13}^{0}e^{\frac{\delta_{13}Fd}{k_{B}T}},$ $\displaystyle
k_{42}(F)=k_{42}^{0}e^{\frac{-\delta_{42}Fd}{k_{B}T}},\quad
k_{24}(F)=k_{24}^{0}e^{\frac{\delta_{24}Fd}{k_{B}T}},$

where $k_{ij}^{0}$ are load free transition rates, $F$ is the external load
(positive if it points to the minus end of MT), $k_{B}$ is the Boltzmann
constant, $T$ is the absolute temperature, $d=8$ nm is the step size of
kinesin. $\delta_{ij}\geq 0$ is called load distribution factor that reflects
how the external load affects the rate of transition from state $i$ to state
$j$, which satisfies

$\delta_{23}+\delta_{32}+\delta_{31}+\delta_{13}=1,\quad\delta_{23}+\delta_{32}+\delta_{42}+\delta_{24}=1.$

Since transitions 4$\to$2 and 3$\to$1 describe the same biochemical process,
during which phosphate Pi is released and then the trailing head is detached
from MT, we let $k_{42}^{0}=k_{31}^{0}$ and $k_{24}^{0}=k_{13}^{0}$. However,
the position of energy barrier between states 4 and 2 might be different from
that between states 3 and 1, so load distribution factors
$\delta_{42},\delta_{24}$ are generally different from
$\delta_{31},\delta_{13}$ Fisher and Kolomeisky (2001); Howard (2001).

The detachment rate $k_{\rm{off}}$ of kinesin from MT is assumed to be load
$F$ dependent and formulated as

$k_{\rm{off}}(F)=k_{\rm{off}}^{0}e^{\frac{Fd_{\rm{off}}}{k_{B}T}},$

with $d_{\rm{off}}\geq 0$ a characteristic distance. This expression is
actually equivalent to the one used previously in Müller et al. (2008); Zhang
and Fisher (2010), where the detachment rate is assumed the be the form
$\epsilon_{0}\exp(|F|/F_{d})$ with parameter $F_{d}$ called detachment force.

It is biophysically reasonable to assume the sliding process from states 6 to
1 is load-dependent, and the rate $k_{61}$ is given as

$k_{61}(F)=k_{61}^{0}e^{\frac{Fd}{k_{B}T}},$

since each sliding process is always coupled with a backward translocation of
$d=8$ nm.

In addition, transitions 1$\to$2 and 3$\to$4 describe ATP binding to the
nucleotide-free head of kinesin. So rates $k_{12}$ and $k_{34}$ depend on ATP
concentration ${\rm[ATP]}$ in environment. For simplicity, we let

$k_{12}({\rm ATP})=k_{12}^{0}{\rm[ATP]}=k_{34}^{0}{\rm[ATP]}=k_{34}({\rm
ATP}).$

The inverse transitions 2$\to$1 and 4$\to$3 describe the disassociation of ATP
molecule from kinesin head. So rates $k_{21}$ and $k_{43}$ are independent of
load $F$ and ATP concentration [ATP]. Also for simplicity, we let
$k_{21}=k_{43}$ in our model.

Rates $k_{25},k_{52},k_{51},k_{15}$ describe the release/recruiting of
phosphate Pi or ADP. For simplicity, we always assume the change of
concentrations of ADP and Pi in environment are negligible. So, all these
rates are assumed to be constants, though more sophisticated methods can be
used to handle their dependence on [ATP] Fisher and Kolomeisky (2001).

##  III Theoretical analysis of the model

###  III.1 Mean velocity and diffusion coefficient

According to the model illustrated in Fig. 1(b), the probability
$p_{j}^{n}(t)$ of finding kinesin in biochemical state $j$, at binding site
$n$, and at time $t$, is governed by the following master equations,

$\displaystyle\begin{aligned}
&\partial_{t}p_{1}^{n}=k_{31}p_{3}^{n-1}+k_{21}p_{2}^{n}+k_{51}p_{5}^{n}+k_{61}p_{6}^{n+1}\\\
&\qquad\quad-(k_{12}+k_{13}+k_{15})p_{1}^{n},\\\
&\partial_{t}p_{2}^{n}=k_{12}p_{1}^{n}+k_{32}p_{3}^{n}+k_{42}p_{4}^{n-1}+k_{52}p_{5}^{n}\\\
&\qquad\quad-(k_{21}+k_{23}+k_{24}+k_{25})p_{2}^{n},\\\
&\partial_{t}p_{3}^{n}=k_{23}p_{2}^{n}+k_{43}p_{4}^{n}+k_{13}p_{1}^{n+1}-(k_{32}+k_{31}+k_{34})p_{3}^{n},\\\
&\partial_{t}p_{4}^{n}=k_{34}p_{3}^{n}+k_{24}p_{2}^{n+1}-(k_{43}+k_{42})p_{4}^{n},\\\
&\partial_{t}p_{5}^{n}=k_{15}p_{1}^{n}+k_{25}p_{2}^{n}+k_{65}p_{6}^{n}-(k_{51}+k_{52}+k_{56})p_{5}^{n},\\\
&\partial_{t}p_{6}^{n}=k_{56}p_{5}^{n}-(k_{61}+k_{65})p_{6}^{n}.\end{aligned}$
(2)

Note that, the detachment process will be added in the analysis of mean run-
length/time. The distance between two adjacent binding sites, i.e. the step
size of kinesin, is $d=8$ nm.

Using Eq. (2), the Fourier transform of probability $p_{j}^{n}(t)$, defined by
$P_{j}(q,t):=\sum_{n=-\infty}^{\infty}p_{j}^{n}(t)e^{-iqn}$, satisfies

$\partial_{t}\bm{P}(q,t)=\bm{M}(q)\bm{P}(q,t),$ (3)

where $\bm{P}(q,t)=[P_{1}(q,t),\cdots,P_{6}(q,t)]^{T}$. Matrix $\bm{M}(q)$
depends on all rates $k_{ij}$ labeled in Fig. 1(b) and the variable $q$, see
Supplemental Material for its detailed elements.

Taking Laplace transform of Eq. (3), we get a $6$th order characteristic
polynomial in $s$,

$\left|s\bm{I}-\bm{M}(q)\right|=s^{6}+\cdots+\alpha(q)s^{2}+\beta(q)s+\gamma(q),$
(4)

with $\bm{I}$ the identity matrix. As stated in Koza (1999, 2000); Chemla et
al. (2008), to obtain expressions of the mean velocity $v$, diffusion
coefficient $D$, coefficients of the last three terms of the characteristic
polynomial are enough. It can be shown that, see Supplemental Material for
detailed derivation,

$v=-\frac{id\dot{\gamma}(0)}{\beta(0)},$ (5)
$D=\frac{d^{2}}{2}\left(\frac{\ddot{\gamma}(0)}{\beta(0)}-\frac{2\dot{\gamma}(0)\dot{\beta}(0)}{\beta(0)^{2}}+\frac{2\alpha(0)\dot{\gamma}(0)^{2}}{\beta(0)^{3}}\right),$
(6)

where ‘$\cdot$’ denotes the derivative with respect to $q$. From Eqs. (5) and
(6), the randomness defined as $r=2D/vd$ is

$r=i\left(\frac{\ddot{\gamma}(0)}{\dot{\gamma}(0)}-2\frac{\dot{\beta}(0)}{\beta(0)}+2\frac{\alpha(0)\dot{\gamma}(0)}{\beta(0)^{2}}\right).$
(7)

Obviously, $v,D$ and $r$ are all functions of transition rates $k_{ij}$, and
depend on external load $F$ and ATP concentration [ATP].

###  III.2 Mean run length/time along MT

Let $\rho_{i}=\lim_{t\to\infty}\sum_{n=-\infty}^{\infty}p_{i}^{n}(t)$ be the
steady state probability of finding the motor in biochemical state $i$, then
from Eq. (2), we have

$M\rho=\bm{0},$ (8)

with matrix $M=\bm{M}(0)$, $\rho=[\rho_{1},\cdots,\rho_{6}]^{T}$, and $\bm{0}$
the zero vector. Generally, $\rho$ can be obtained from Eq. (8) combined with
the normalization condition $\sum_{i=1}^{6}\rho_{i}=1$.

Denoting $F_{m}^{n_{j}}(t)$ the probability density of kinesin, which is
initiated from site $n$ and biochemical state $j$ at time $t=0$ and detaches
for the first time from site $m$ of MT at time $t$. It can be shown that
$F_{m}^{n_{j}}(t)$ satisfy the following backward master equations, see Kampen
(2007); Gardiner (2010); Kolomeisky (2015); Zhang and Kolomeisky (2018),

$\displaystyle\partial_{t}F_{m}^{n_{1}}=k_{13}F_{m}^{(n-1)_{3}}+k_{12}F_{m}^{n_{2}}+k_{15}F_{m}^{n_{5}}$
(9) $\displaystyle\qquad\qquad-(k_{12}+k_{13}+k_{15})F_{m}^{n_{1}},$
$\displaystyle\partial_{t}F_{m}^{n_{2}}=k_{21}F_{m}^{n_{1}}+k_{23}F_{m}^{n_{3}}+k_{24}F_{m}^{(n-1)_{3}}+k_{25}F_{m}^{n_{5}}$
$\displaystyle\qquad\qquad-(k_{21}+k_{23}+k_{24}+k_{25})F_{m}^{n_{2}},$
$\displaystyle\partial_{t}F_{m}^{n_{3}}=k_{31}F_{m}^{(n+1)_{1}}+k_{32}F_{m}^{n_{2}}+k_{34}F_{m}^{n_{4}}$
$\displaystyle\qquad\qquad-(k_{31}+k_{32}+k_{34})F_{m}^{n_{3}},$
$\displaystyle\partial_{t}F_{m}^{n_{4}}=k_{42}F_{m}^{(n+1)_{2}}+k_{43}F_{m}^{n_{3}}-(k_{42}+k_{43})F_{m}^{n_{4}},$
$\displaystyle\partial_{t}F_{m}^{n_{5}}=k_{51}F_{m}^{n_{1}}+k_{52}F_{m}^{n_{2}}+k_{56}F_{m}^{n_{6}}$
$\displaystyle\qquad\qquad-(k_{51}+k_{52}+k_{56})F_{m}^{n_{5}},$
$\displaystyle\partial_{t}F_{m}^{n_{6}}=k_{61}F_{m}^{(n-1)_{1}}+k_{65}F_{m}^{n_{5}}+k_{\textrm{off}}\delta_{nm}\delta(t)$
$\displaystyle\qquad\qquad-(k_{61}+k_{65}+k_{\textrm{off}})F_{m}^{n_{6}},$

where $\delta_{nm}=0$ if $n\neq m$, and $\delta_{nn}=1$. $\delta(t)$ in the
last equation means that if kinesin detaches from MT, the first-passage
process is accomplished immediately.

From Eq. (9), the conditional mean run length of kinesin which starts from
biochemical state $i$,

$l_{i}^{(1)}:=\sum_{m}m\int_{0}^{+\infty}F_{m}^{0_{i}}(t)\mathrm{d}t,$

satisfies $K\bm{l}^{(1)}=\bm{b}^{(1)}$, with
$\bm{l}^{(1)}=[l_{1}^{(1)},\cdots,l_{6}^{(1)}]^{T}$,
$\bm{b}^{(1)}=[k_{13},k_{24},-k_{31},-k_{42},0,k_{61}]^{T}$, and matrix $K$ is
similar to the coefficient matrix of Eq. (9), for details see Supplemental
Material. The mean run length of kinesin along MT is then obtained by

$\displaystyle\langle
l\rangle=\sum_{i=1}^{6}\rho_{i}l_{i}^{(1)}=\rho\cdot\bm{l}^{(1)}.$ (10)

In the same way, the second moment of conditional run length

$l_{i}^{(2)}:=\sum_{m}m^{2}P\int_{0}^{+\infty}F_{m}^{0_{i}}(t)\mathrm{d}t,$

satisfies $K\bm{l}^{(2)}=2\bm{b}^{(1)}\circ\bm{l}-\bm{b}^{(2)}$, where
‘$\circ$’ is the Hadamard product,
$\bm{l}^{(2)}=[l_{1}^{(2)},\cdots,l_{6}^{(2)}]^{T}$,
$\bm{l}=[l_{3}^{(1)},l_{4}^{(1)},l_{1}^{(1)},l_{2}^{(1)},0,l_{1}^{(1)}]^{T}$,
$\bm{b}^{(2)}=[k_{13},k_{24},k_{31},k_{42},0,k_{61}]^{T}$ is the absolute
value of $\bm{b}^{(1)}$. So the variance of run length is

$\displaystyle{\rm Var}(l)$ $\displaystyle=$ $\displaystyle\langle
l^{2}\rangle-\langle
l\rangle^{2}=\sum_{i=1}^{6}\rho_{i}l_{i}^{(2)}-\left(\sum_{i=1}^{6}\rho_{i}l_{i}\right)^{2}$
(11) $\displaystyle=$
$\displaystyle\rho\cdot\bm{l}^{(2)}-(\rho\cdot\bm{l}^{(1)})^{2}.$ (12)

The mean run time of kinesin initiated from biochemical state $i$ is

$t_{i}^{(1)}:=\int_{0}^{+\infty}t\sum_{m}F_{m}^{0_{i}}(t)\mathrm{d}t.$

Using Eq. (9), $\bm{t}^{(1)}:=[t_{1}^{(1)},\cdots,t_{6}^{(1)}]^{T}$ satisfies
$K\bm{t}^{(1)}=-\bm{e}$ with $\bm{e}=[1,\cdots,1]^{T}$. So the mean run time
of kinesin along MT is

$\displaystyle\langle
t\rangle=\sum_{i=1}^{6}\rho_{i}t_{i}^{(1)}=\rho\cdot\bm{t}^{(1)}.$ (13)

Similarly, the variance of run time is

$\displaystyle{\rm Var}(t)$ $\displaystyle=$ $\displaystyle\langle
t^{2}\rangle-\langle
t\rangle^{2}=\sum_{i=1}^{6}\rho_{i}t_{i}^{(2)}-\left(\sum_{i=1}^{6}\rho_{i}t_{i}\right)^{2}$
(14) $\displaystyle=$
$\displaystyle\rho\cdot\bm{t}^{(2)}-(\rho\cdot\bm{t}^{(1)})^{2},$ (15)

where $\bm{t}^{(2)}=[t_{1}^{(2)},\cdots,t_{6}^{(2)}]^{T}$ satisfies
$K\bm{t}^{(2)}=-2\bm{t}^{(1)}$, see Gardiner (2010).

Figure 2:  Theoretical predictions (solid lines) and experimental data
(markers) of various biophysical properties of kinesin-1 purified from wild-
type Drosophila. The data in (a,b,d,e) are from Schnitzer et al. (2000), and
in (c,f) are from Visscher et al. (1999b), with $D$ is estimated from
$r=2D/vd$. In (d) the left axis is for 5 $\mu$M ATP while the right axis is
for 2 mM ATP. Theoretical results are obtained from formulations given in Eqs.
(5,6,7,10), with model parameters listed Tab. 1.

##  IV Mechanochemical properties of kinesin-1

###  IV.1 Biophysics of kinesin motion along MT

Using parameter values listed in Tab. 1, our model can reproduce the related
experimental results well, see Fig. 2. In the following, we will always use
parameter values listed in Tab. 1 to discuss the mechanochemical properties of
kinesin-1.

Table 1:  Model parameter values obtained by fitting to experimental data of kinesin-1 purified from wild-type Drosophila measured in Schnitzer et al. (2000); Visscher et al. (1999b), see Figs. 1 and 2.  Parameter | Value | Parameter | Value
---|---|---|---
$k_{12}^{0},k_{34}^{0}$ | 9.827 $\mu$M${}^{-1}\cdot$s-1 | $k_{21},k_{43}$ | 5074.875 s-1
$k_{23}^{0}$ | 1627.099 s-1 | $\delta_{23}$ | 0.055
$k_{32}^{0}$ | 0.006 s-1 | $\delta_{32}$ | 0.416
$k_{31}^{0},k_{42}^{0}$ | 137.582 s-1 | $\delta_{31}$ | 0
$\delta_{42}$ | 0.006 | $k_{13}^{0},k_{24}^{0}$ | 1.666 s-1
$\delta_{13}$ | 0.529 | $\delta_{24}$ | 0.523
$k_{25}$ | 5146.371 s-1 | $k_{52}$ | 77.252 s-1
$k_{15}$ | 4.344 s-1 | $k_{51}$ | 1455.909 s-1
$k_{56}$ | 345.215 s-1 | $k_{65}$ | 2923.252 s-1
$k_{61}^{0}$ | 0.001 s-1 | $k_{\rm off}^{0}$ | 5.286 s-1
$d_{\rm{off}}$ | 0.803 nm |  |

From Tab. 1, $\delta_{23}+\delta_{32}\approx 0.47$, which implies that there
may be two substeps within each 8 nm step of kinesin, with size 3.76 nm and
4.23 nm respectively, as previously observed in Coppin et al. (1996);
Nishiyama et al. (2001). According to Fig. 1(a), these two substeps are caused
by binding of the nucleotide-free head to the front site of MT and detachment
of the ADP bound rear head from MT, respectively. But, for high ATP
concentration, rates $k_{12},k_{34}$ and $k_{23}$ are very large, the substep
coupled with transition $2\rightleftharpoons 3$ is difficult to be observed.

The load distribution factor $\delta_{31}=0$ means the energy barrier between
state 3 and 1 is very close to state 3. Transition rate $k_{31}$ is
independent of external load $F$ while $k_{13}$ increases rapidly with $F$.
The process $4\rightleftharpoons 2$ is similar since $\delta_{42}=0.006$ is
also very small.

Figs. 3(a,d) show that, at low load $F$ and low ATP concentration [ATP],
kinesin mainly stays in biochemical state 1, i.e., with only one nucleotide-
free head bound to MT while the other ADP bound head detached. At high ATP
concentration, kinesin stays mainly in state 4, with both heads in ATP or
ADP$\cdot$Pi bound state and bound to MT tightly. Meanwhile, if the load $F$
is not high, the probability that kinesin in the two ADP bound state 5 is also
relatively high (about 1/5), see Figs. 3(a,b,e). At low [ATP] while high load
$F$, kinesin stays mainly in state 3 with both two heads bound to MT but one
in nucleotide-free state, see Figs. 3(c,d). Provided the load $F$ is high,
kinesin stays mainly in two head bound states, and stays mainly in state 4
when [ATP] is high while in state 3 otherwise, see Figs. 3(c).

Generally, probability $\rho_{1}$ decreases while
$\rho_{2},\rho_{4},\rho_{5},\rho_{6}$ increases with ATP concentration. At
high load $F$, $\rho_{3}$ decreases with [ATP] monotonically, otherwise
$\rho_{3}$ may increase first and then decrease with [ATP]. At low ATP
concentration, kinesin mainly stays in the ATP waiting state 1 or 3. But with
the increase of load $F$ kinesin switches from the one head bound state 1 to
the two heads bound state 3, see Figs. 3(a-d). When ATP concentration is high,
all probabilities $\rho_{i}$ change only slightly with the load $F$, see Fig.
3(e).

Figure 3: (a-c) Probability $\rho_{i}$ of kinesin in the six biochemical
states as depicted in Fig. 1, versus ATP concentration [ATP] at different
load. (d-e) Probability $\rho_{i}$ of kinesin versus load $F$ at different ATP
concentration. (f) Stall force of kinesin as a function of [ATP], where solid
line shows the theoretical predictions using parameters listed in Tab. 1, and
triangles denote experimental data measured in Visscher et al. (1999b).

Fig. 3(f) shows the stall force $F_{s}$, under which the mean velocity of
kinesin vanishes, increases slightly with ATP concentration, and eventually
approaching to a constant between 7 and 8 pN, which agrees well with
experimental measurements Coppin et al. (1997); Visscher et al. (1999b);
Fisher and Kolomeisky (2001).

Figure 4: (a,b) Mean run time $\langle t\rangle$ with the change of load $F$
and [ATP], see Eq. (13). (c) Variance of run length ${\rm Var}(l)$ with the
change of $F$, see Eq. (11). (d) Variance of run time ${\rm Var}(t)$ with the
change of [ATP], see Eq. (14). (e,f) Probability flux of kinesin detachment
from MT, $\textrm{flux}_{\rm off}:=k_{\rm off}\rho_{6}$, with the change of
load $F$ and [ATP]. In (a,e), left axis is for [ATP]$=5$ $\mu$M, and right
arix is for [ATP]$=2$ mM.

Figs. 4(a,b) show mean run time $\langle t\rangle$ of kinesin before its
detachment from MT decreases with ATP concentration monotonically, while first
decreases slightly and then increases rapidly with load $F$. These results
seem surprising. One possible explanation is that, in general, kinesin does
not detach from MT unless it has transported its cargo to corresponding
destination Toba et al. (2006); Toprak et al. (2009). So, with low ATP
concentration, kinesin will stay on MT for more time, since the period of
single biochemical cycle becomes long due to the lack of ATP molecule.
Actually, Fig. 4(b) and Fig. 2(b) show both mean run time $\langle t\rangle$
and mean run length $\langle l\rangle$ tend to limit constants at high ATP
concentration. One possible reason is that, in cells this run time/length is
enough to transport any cargoes to their corresponding destinations. At low
ATP concentration and high external load, $\langle t\rangle$ increases with
$F$, which implies kinesin will try to run more time along MT to complete its
task when its motion is slowed down by resistance, see Figs. 4(a,b).

Generally, both mean value $\langle l\rangle$ and variance ${\rm Var}(l)$ of
kinesin run length along MT decrease with load $F$, and if $F$ is not too
large both $\langle l\rangle$ and ${\rm Var}(l)$ increases with ATP
concentration, see Figs. 2(b,e), 4(c), and S1(d). The variance of run time
${\rm Var}(t)$ decreases with ATP concentration, while decreases first and
then increases with the load $F$, see Figs. 4(d), S1(a). Further calculations
show both coefficient of variation ${\rm CV}_{l}:=\langle l\rangle/\sqrt{{\rm
Var}(l)}$ and ${\rm CV}_{t}:=\langle t\rangle/\sqrt{{\rm Var}(t)}$ almost
equal to 1, see Fig. S1(b,c,e,f). So both the run time of kinesin along MT and
its run length are approximately exponentially distributed.

The mean run time $\langle t\rangle$ of kinesin along MT is actually the time
spent by kinesin in average to detach from MT. So $\langle t\rangle$ is
related to the probability flux of detachment $\textrm{flux}_{\rm off}:=k_{\rm
off}\rho_{6}$. Intuitively, $\textrm{flux}_{\rm off}$ should change with $F$
and [ATP] in a similar manner as $1/\langle t\rangle$. In contrast to Fig.
4(a,b), the plots in Fig. 4(e,f) show detachment flux $\textrm{flux}_{\rm
off}$ increases monotonically with [ATP], while increases first and then
decreases with load $F$.

Figure 5:  Probability flux of forward/backward/futile cycle of kinesin-1, as
well as their proportions, versus ATP concentration [ATP], with load $F=0$ pN,
4 pN and 8 pN, respectively.

###  IV.2 Pathways of forward/backward/futile biochemical cycle of kinesin-1

Given the model depicted in Fig. 1, there is no evident information about the
coupling between biochemical cycle of ATP hydrolysis and mechanical step of
kinesin. As pointed out before, there are generally three categories of the
biochemical cycle, forward/backward/futile cycle. Mean velocity $v$ is the sum
of probability fluxes of different biochemical cycles, weighted by mechanical
step. It is generally difficult to sort out these three categories of
biochemical cycle and then calculate their corresponding fluxes, since net
flux between any two adjacent biochemical states $i$ and $j$,
$\textrm{flux}_{ij}:=k_{ij}\rho_{i}-k_{ji}\rho_{j}$, usually changes with
external load $F$ and ATP concentration [ATP]. Fortunately, thorough numerical
calculations show that, with parameter values listed in Tab. 1 and
biophysically meaningful load $F$ and [ATP], flux variation in our model for
kinesin-1 consists of only 5 different cases, in each case the direction of
net flux between any two adjacent states does not change, see Figs. S2-S4. For
each case, forward/backward/futile cycle can be phenomenologically sorted out
by analyzing all possible biochemical pathways. The results are summarized in
Tab. 2.

Conditions | Cycles | Pathway | ATP | Steps | Total flux
---|---|---|---|---|---
$++++$ (Case 1) | Forward | $\overrightarrow{1231}$ | 1 | 1 | $\textrm{flux}_{23}$
$\overrightarrow{2342}$ | 1 | 1
Backward | $\overrightarrow{12561}$ | 1 | -1 | $\textrm{flux}_{61}$
Futile | $\overrightarrow{1251}$ | 1 | 0 | $\textrm{flux}_{51}$
$-+++$ (Case 2) | Forward | $\overrightarrow{2342}$ | 1 | 1 | $\textrm{flux}_{23}$
Backward | $\overrightarrow{12561}$ | 1 | -1 | $\textrm{flux}_{61}$
$\overrightarrow{1342561}$ | 1 | -1
Futile | $\overrightarrow{1251}$ | 1 | 0 | $\textrm{flux}_{51}$
$\overrightarrow{134251}$ | 1 | 0
$--++$ (Case 3) | Forward | — | — | — | $0$
Backward | $\overrightarrow{12561}$ | 1 | -1 | $\textrm{flux}_{61}+\textrm{flux}_{51}\frac{\textrm{flux}_{32}}{\textrm{flux}_{25}}$
$\overrightarrow{1342561}$ | 1 | -1
$\overrightarrow{13251}$ | 0 | -1
$\overrightarrow{132561}$ | 0 | -2
Futile | $\overrightarrow{1251}$ | 1 | 0 | $\textrm{flux}_{51}\frac{\textrm{flux}_{12}+\textrm{flux}_{42}}{\textrm{flux}_{25}}$
$\overrightarrow{134251}$ | 1 | 0
$---+$ (Case 4) | Forward | — | — | — | $0$
Backward | $\overrightarrow{12561}$ | 1 | -1 |  $\textrm{flux}_{61}+\textrm{flux}_{24}$ +$\textrm{flux}_{51}\frac{\textrm{flux}_{13}}{\textrm{flux}_{25}}$
$\overrightarrow{132561}$ | 0 | -2
$\overrightarrow{13251}$ | 0 | -1
$\overrightarrow{2432}$ | -1 | -1
Futile | $\overrightarrow{1251}$ | 1 | 0 | $\textrm{flux}_{51}\frac{\textrm{flux}_{12}}{\textrm{flux}_{25}}$
$----$ (Case 5) | Forward | — | — | — | $0$
Backward | $\overrightarrow{1321}$ | -1 | -1 | $\textrm{flux}_{32}+\textrm{flux}_{61}$
$\overrightarrow{132561}$ | 0 | -2
$\overrightarrow{2432}$ | -1 | -1
$\overrightarrow{13251}$ | 0 | -1
Futile | — | — | — | $0$
Table 2: Pathway details of forward/backward/futile biochemical cycle of
kinesin-1. Column ‘Conditions’ lists the signs of
$\textrm{flux}_{31},\textrm{flux}_{23},\textrm{flux}_{34}$, and
$\textrm{flux}_{12}$, respectively. Column ‘ATP’ lists the number of ATP
molecule consumed in the corresponding pathway, where ‘$-1$’ means one ATP is
synthesized. Column ‘Steps’ gives the number of mechanical step coupled with
the pathway, with one step 8 nm. Column ‘Total flux’ lists the total flux of
the forward/backward/futile cycle. $\overrightarrow{ij\cdots ki}$ denotes the
biochemical pathway $i\\!\to\\!j\\!\to\\!\cdots\\!\to\\!k\\!\to\\!i$. It is
evident from Fig. 1 that $\textrm{flux}_{42}=\textrm{flux}_{34}$. Flux
$\textrm{flux}_{25}$, $\textrm{flux}_{56}=\textrm{flux}_{61}$, and
$\textrm{flux}_{51}$ are always nonnegative for $0\leq F\leq 9$ pN and 1
$\mu$M $\leq$ [ATP] $\leq$ 10000 $\mu$M, see Fig. S2.

Results show that, at low load $F$ the flux of backward cycle
$\textrm{flux}_{-}$ is negligible, indicating the impossibility of backward
motion Toba et al. (2006); Toprak et al. (2009), see Figs. 5(a,b,d,e) and Fig.
S5. Both forward flux $\textrm{flux}_{+}$ and futile flux $\textrm{flux}_{0}$
increase with ATP concentration [ATP]. At low load $F$, proportions of
$\textrm{flux}_{+}$ and $\textrm{flux}_{0}$ are almost independent of [ATP],
which implies that the utilization ratio of ATP is independent of [ATP] at low
load, see Figs. 5(d,e). At super stall load, kinesin motion is dominated by
backward cycle when [ATP] is small, while dominated by futile cycle when [ATP]
is large. Particularly, for $F=8$ pN, the proportion of backward flux
$\textrm{flux}_{-}$ is almost 1 when [ATP] $<22.23$ $\mu$M, indicating that
kinesin-1 moves backward almost deterministically. For large enough [ATP], the
lower limit value of the proportion of $\textrm{flux}_{-}$ is almost the same
as the upper limit value of the proportion of $\textrm{flux}_{+}$, so the
motion of kinesin is completed stalled, see Figs. 5(c,f). With increasing
[ATP], the proportion of $\textrm{flux}_{0}$ rises rapidly, until 263.47
$\mu$M, at which the sign of $\textrm{flux}_{23}$ changes from negative to
positive, switching from Case 3 to Case 2, producing the flux of forward cycle
$\textrm{flux}_{+}>0$, see Fig. S3(c) and Figs. 5(c,f). On the whole, at high
load and low ATP concentration, kinesin-1 can only move backward. As [ATP]
increases, the futile flux begins to rise rapidly, and then forward flux is
generated. For large enough [ATP], backward flux is balanced by forward flux,
and the motion of kinesin-1 is stalled completely.

Further calculations show that both forward flux $\textrm{flux}_{+}$ and its
proportion decrease, while the proportion of $\textrm{flux}_{-}$ increases,
with load $F$, which implies the utilization ratio of ATP decreases with load
$F$, see Fig. S5. At [ATP] = 5 $\mu$M and [ATP] = 2 mM, the intersection
points of $\textrm{flux}_{+}$ and $\textrm{flux}_{-}$ are found at $F=5.85$ pN
and $F=7.76$ pN, respectively, which are exactly the stall forces $F_{s}$ of
kinesin at [ATP] = 5 $\mu$M and [ATP] = 2 mM, see Fig. 3(f). One surprising
result is that, except the extreme cases with high load $F$ and very low
[ATP], the biochemical process of kinesin-1 is mainly dominated by futile
cycle, and utilization ratio of ATP is only about 20%, see Fig. 5. At small
load $F$, the backward flux is almost zero. While for load $F$ larger than
stall force, the forward flux decreases rapidly to zero. This indicates that,
unless at external load which is near the stall force, the motion of kinesin-1
is almost deterministic, and at small load, kinesin-1 hardly steps backward.
This is different with the conclusion draw from Brownian ratchet models
Astumian (1997); Zhang (2008).

Tab. 2 shows that forward cycle consists of two possible pathways,
$\overrightarrow{1231}$ and $\overrightarrow{2342}$, both of them include
hydrolysis of one ATP and a forward mechanical step of 8 nm. Futile cycle
consists of two possible pathways, $\overrightarrow{1251}$ and
$\overrightarrow{134251}$, both of them include hydrolysis of one ATP but
without mechanical step. Backward cycle consists of six possible pathways,
which are $\overrightarrow{12561}$, $\overrightarrow{1342561}$,
$\overrightarrow{13251}$, $\overrightarrow{132561}$, $\overrightarrow{2432}$
and $\overrightarrow{1321}$. These six pathways can be further classified into
three categories, (i) $\overrightarrow{12561}$ and $\overrightarrow{1342561}$,
which are mainly induced by backward sliding 6$\to$1 through the semi-detach
state, (ii) $\overrightarrow{13251}$, $\overrightarrow{2432}$, and
$\overrightarrow{1321}$, which can be roughly regarded as the reversal of
corresponding forward pathways, and (iii) $\overrightarrow{132561}$, which
includes both the reversal of a forward pathway and the backward sliding
6$\to$1, and therefore coupled with a backward motion of total 16 nm. Pathway
$\overrightarrow{2432}$ and $\overrightarrow{1321}$ include the synthesis of
ATP, while $\overrightarrow{12561}$ includes the hydrolysis of ATP. Note that
pathway $\overrightarrow{ij\cdots ki}$ can be denoted equivalently as
$\overrightarrow{j\cdots kij}$, $\overrightarrow{kij\cdots k}$, etc. For
convenience, we denote the probability flux through pathway
$\overrightarrow{ij\cdots ki}$ by $f_{\overrightarrow{ij\cdots ki}}$.

With small [ATP] and load $F\ll F_{s}$, the forward cycle is mainly realized
through pathway $\overrightarrow{1231}$, see Figs. S6(a,d), S8(a,d), and
S9(a,d). Otherwise, the forward cycle is mainly realized through pathway
$\overrightarrow{2342}$, Figs. S6-S10(a,d). During pathway
$\overrightarrow{2342}$, binding of ATP to the front head of kinesin is
earlier than the release of phosphate Pi from the rear head, while during
pathway $\overrightarrow{1231}$, ATP binds to the front head only after the
detachment of rear head from MT. So our results imply that at low external
load and low ATP concentration, the MT bound head might be in the nucleotide-
free state, while at high ATP concentration or high load, the head bound to MT
is always in ATP or ADP$\cdot$Pi binding state. This is consistent with
results shown in Fig. 3.

Except some extreme cases, the futile cycle is mainly realized through the
pathway $\overrightarrow{1251}$. Moreover, for small load $F$ or small [ATP],
the probability that the futile cycle realized through
$\overrightarrow{134251}$ is almost negligible, see Figs. S6-S10(b,e).

If load $F$ is small, the backward cycle is mainly realized through pathway
$\overrightarrow{12561}$, and other pathways are all negligible, see Figs.
S6-S9(c,f). If both $F$ and [ATP] are large, the backward cycle is mainly
realized through pathway $\overrightarrow{12561}$ and pathway
$\overrightarrow{1342561}$, and also with $\overrightarrow{12561}$ the most
prominent one, see Figs. S7(c,f) and S10(c,f). If $F$ is large but [ATP] is
small, pathways $\overrightarrow{13251}$, $\overrightarrow{2432}$,
$\overrightarrow{1321}$ and $\overrightarrow{132561}$ are all important, while
both pathway $\overrightarrow{12561}$ and pathway $\overrightarrow{1342561}$
are negligible, see Figs. S6(c,f) and S10(e,f). For large $F$ and small [ATP],
prominent pathways of backward cycle are roughly related to the reversal of
corresponding pathways of forward cycle, and consequently might lead to the
synthesis of ATP, see Tab. 2. However, this is difficult to observe
experimentally, since kinesin-1 may have already detached from MT before
achieving such harsh conditions.

Figure 6:  Probability flux of 8 nm backstep, 8 nm backslip and 16 nm
backslip, as well as their corresponding proportions at various [ATP] and load
$F$. (a-c) show the probability flux at [ATP] = 5 $\mu$M, [ATP] = 2 mM, and
$F=8$ pN, respectively, and (d-f) show their corresponding proportions. The
flux of 8 nm backstep is obtained by
$f_{\overrightarrow{13251}}+f_{\overrightarrow{2432}}+f_{\overrightarrow{1321}}$,
the flux of 8 nm backslip is obtained by
$f_{\overrightarrow{12561}}+f_{\overrightarrow{1342561}}$, and the flux of 16
nm backslip is $f_{\overrightarrow{132561}}$.

During pathways $\overrightarrow{12561}$ and $\overrightarrow{1342561}$,
backward motion of kinesin-1 is induced only by the backward sliding of 8 nm
through transition 6$\to$1\. But during pathway $\overrightarrow{132561}$, the
backward motion is accomplished by both one backward sliding of 8 nm and one
backward step of 8 nm. Fig. 6 shows that the probability flux through pathway
$\overrightarrow{132561}$ is very rare and the backward step of 16 nm is
observable only under the very harsh condition that external load $F$ is high
but ATP is scarce. If ATP concentration reach saturating, both probability
flux of 8 nm backstep and 16 nm backslip decrease to zero, only backslip 8nm
existing, see Figs. 6(b,c).

##  V Conclusions and remarks

According to recent experiments, a mechanochemical model of kinesin-1 is
constructed, in which the motor can step forward/backward through multiple
biochemical pathways. Biophysical quantities of kinesin, including its mean
velocity, diffusion constant, and mean run time/length along microtubule, can
be obtained theoretically.

Our study shows that mean run time of kinesin before its detachment from
microtubule decreases with ATP concentration monotonically and finally tends
to an external load-dependent limit value, but decreases first and then
increases rapidly with external load, which are different with properties of
mean run length observed in experiments. However, both distributions of run
time and run length are approximately exponential.

Forward motion of kinesin-1 can be realized through two biochemical pathways.
At high ATP concentration or under high external load, the kinesin head bound
to microtubule will always be in ATP or ADP$\cdot$Pi binding state. But at low
ATP concentration and low load, new ATP molecule will not bind to the
nucleotide-free head before the release of phosphate Pi and the detachment of
the trailing head from microtubule. Kinesin-1 may stay at the one head bound
state for a long time to wait for ATP arrival.

Backward motion of kinesin-1 is mainly caused by backward sliding along
microtubule through semi-detach state, and hardly through directional reversal
of forward pathways as usually believed in previous studies. Large backward
step and directional reversal of forward pathways can only happen in very
harsh environments with scarce ATP molecules but under high external load.
Usually, kinesin-1 has already been pulled down from microtubule before
achieving such conditions.

One surprising finding is that the utilization ratio of ATP is only about 20%,
and about 80% of ATP is consumed in futile cycle, which consists of two
possible biochemical pathways. Unless ATP concentration is extremely low but
external load is high, ATP is required in almost all backward steps of
kinesin-1, though these backward steps are mostly accomplished through the
sliding process in its semi-detach state.

## References

  * Howard (2001) J. Howard, _Mechanics of Motor Proteins and the Cytoskeleton_ (Sinauer Associates and Sunderland, MA, 2001).
  * Vale (2003) R. D. Vale, Cell 112, 467 (2003).
  * Schliwa (2003) M. Schliwa, _Molecular Motors_ (Wiley-Vch, Weinheim, 2003).
  * Kolomeisky (2015) A. B. Kolomeisky, _Motor Proteins and Molecular Motors_ (CRC Press, 2015).
  * Case et al. (2000) R. B. Case, S. Rice, C. L. Hart, B. Ly, and R. D. Vale, Curr. Biol. 380, 157 (2000).
  * Al-Bassam et al. (2003) J. Al-Bassam, Y. Cui, D. Klopfenstein, B. O. Carragher, R. D. Vale, and R. A. Milligan, J. Cell. Biol. 163, 743 (2003).
  * Tomishige et al. (2006) M. Tomishige, N. Stuurman, and R. D. Vale, Nat. Struct. Mol. Biol. 13, 887 (2006).
  * Kalchishkova and Böhm (2008) N. Kalchishkova and K. J. Böhm, J. Mol. Biol. 382, 127 (2008).
  * Hariharan and Hancock (2009) V. Hariharan and W. O. Hancock, Cell. Mol. Bioeng. 2, 177 (2009).
  * Asbury et al. (2003) C. L. Asbury, A. N. Fehr, and S. M. Block, Science 302, 2130 (2003).
  * Yildiz et al. (2004) A. Yildiz, M. Tomishige, R. D. Vale, and P. R. Selvin, Science 303, 676 (2004).
  * Schnitzer and Block (1997) M. J. Schnitzer and S. M. Block, Nature 388, 386 (1997).
  * Coy et al. (1999) D. L. Coy, M. Wagenbach, and J. Howard, J. of Biol. Chem. 274, 3667 (1999).
  * Svoboda et al. (1993) K. Svoboda, C. F. Schmidt, B. J. Schnapp, and S. M. Block, Nature 365, 721 (1993).
  * Higuchi et al. (2004) H. Higuchi, C. E. Bronner, H. W. Park, and S. A. Endow, EMBO J. 23, 2993 (2004).
  * Guydosh and Block (2009) N. R. Guydosh and S. M. Block, Nature 461, 125 (2009).
  * Carter and Cross (2005) N. J. Carter and R. A. Cross, Nature 435, 308 (2005).
  * Toprak et al. (2009) E. Toprak, A. Yildiz, M. T. Hoffman, S. S. Rosenfeld, and P. R. Selvin, Proc. Natl. Acad. Sci. USA 106, 12717 (2009).
  * Clancy et al. (2011a) B. E. Clancy, W. M. Behnke-Parks, J. O. Andreasson, S. S. Rosenfeld, and S. M. Block, Nat. Struct. Mol. Biol. 18, 1020 (2011a).
  * Ariga et al. (2018) T. Ariga, M. Tomishige, and D. Mizuno, Phys. Rev. Lett. 121, 218101 (2018).
  * Svoboda and Block (1994) K. Svoboda and S. M. Block, Cell 77, 773 (1994).
  * Visscher et al. (1999a) K. Visscher, M. J. Schnitzer, and S. M. Block, Nature 400, 184 (1999a).
  * Block et al. (2003) S. M. Block, C. L. Asbury, J. W. Shaevitz, and M. J. Lang, Proc. Natl. Acad. Sci. USA 100, 2351 (2003).
  * Parrondo and Cisneros (2002) J. M. R. Parrondo and B. J. D. Cisneros, Appl. Phys. A 75, 179 (2002).
  * Kolomeisky and Fisher (2007) A. B. Kolomeisky and M. E. Fisher, Ann. Rev. Phys. Chem. 58, 675 (2007).
  * Chowdhury (2013) D. Chowdhury, Physics Reports 529, 1 (2013).
  * Mugnai et al. (2020) M. L. Mugnai, C. Hyeon, M. Hinczewski, and D. Thirumalai, Rev. Mod. Phys. 92, 025001 (2020).
  * Kampen (2007) N. G. V. Kampen, _Stochastic Processes in Physics and Chemistry_ (Elsevier Science & Technology Books, 2007).
  * Gardiner (2010) C. W. Gardiner, _Handbook of Stochastic Methods_ (Springer, 2010).
  * Toleikis et al. (2020) A. Toleikis, N. J. Carter, and R. A. Cross, Biophys. J. 119, 1984 (2020).
  * Block (2007) S. M. Block, Biophys. J. 92, 2986 (2007).
  * Clancy et al. (2011b) B. E. Clancy, W. M. Behnke-Parks, J. O. Andreasson, S. S. Rosenfeld, and S. M. Block, Nat. Struct. Mol. Biol. 18, 1020 (2011b).
  * Fisher and Kolomeisky (2001) M. E. Fisher and A. B. Kolomeisky, Proc. Natl. Acad. Sci. USA 98, 7748 (2001).
  * Zhang (2012) Y. Zhang, PLoS ONE 7, e32717 (2012).
  * Müller et al. (2008) M. J. I. Müller, S. Klumpp, and R. Lipowsky, Proc. Natl. Acad. Sci. USA 105, 4609 (2008).
  * Zhang and Fisher (2010) Y. Zhang and M. E. Fisher, Phys. Rev. E 82, 011923 (2010).
  * Koza (1999) Z. Koza, J. Phys. A 32, 7637 (1999).
  * Koza (2000) Z. Koza, Physica A 285, 176 (2000).
  * Chemla et al. (2008) Y. R. Chemla, J. R. Moffitt, and C. Bustamante, J. Phys. Chem. B 112, 6025 (2008).
  * Zhang and Kolomeisky (2018) Y. Zhang and A. B. Kolomeisky, J. Phys. Chem. B 122, 3272 (2018).
  * Schnitzer et al. (2000) M. J. Schnitzer, K. Visscher, and S. M. Block, Nat. Cell Biol. 2, 718 (2000).
  * Visscher et al. (1999b) K. Visscher, M. J. Schnitzer, and S. M. Block, Nature 400, 184 (1999b).
  * Coppin et al. (1996) C. M. Coppin, J. T. Finer, J. A. Spudich, and R. D. Vale, Proc. Natl. Acad. Sci. USA 93, 1913 (1996).
  * Nishiyama et al. (2001) M. Nishiyama, E. Muto, Y. Inoue, T. Yanagida, and H. Higuchi, Nat. Cell. Biol. 3, 425 (2001).
  * Coppin et al. (1997) C. M. Coppin, D. W. Pierce, L. Hsu, and R. D. Vale, Proc. Natl. Acad. Sci. USA 94, 8539 (1997).
  * Toba et al. (2006) S. Toba, T. M. Watanabe, L. Yamaguchi-Okimoto, Y. Y. Toyoshima, and H. Higuchi, Proc. Natl. Acad. Sci. USA 103, 5741 (2006).
  * Astumian (1997) R. D. Astumian, Science 276, 917 (1997).
  * Zhang (2008) Y. Zhang, Biophys. Chem. 136, 19 (2008).


# JAMMING SOURCE LOCALIZATION USING AUGMENTED PHYSICS-BASED MODEL

###### Abstract

Monitoring interferences to satellite-based navigation systems is of paramount
importance in order to reliably operate critical infrastructures, navigation
systems, and a variety of applications relying on satellite-based positioning.
This paper investigates the use of crowdsourced data to achieve such detection
and monitoring at a central node that receives the data from an arbitrary
number of agents in an area of interest. Under ideal conditions, the pathloss
model is used to compute the Cramér-Rao Bound of accuracy as well as the
corresponding maximum likelihood estimator. However, in real scenarios where
obstructions and reflections are common, the signal propagation is far more
complex than the pathloss model can explain. We propose to augment the
pathloss model with a data-driven component, able to explain the complexities
of the propagation channel. The paper shows a general methodology to jointly
estimate the interference location and the parameters of the augmented model,
showing superior performances in complex scenarios such as those encountered
in urban environments.

Index Terms—  Jamming localization, GNSS, augmented physics-based model,
neural networks.

##  1 Introduction

GNSS is a general term that encompasses satellite-based navigation systems
such as GPS, Galileo, Glonass and Beidou among others [1, 2]. It relies on a
constellation of satellites synchronously emitting known signals, which enable
computation of a receiver position, velocity, and time (PVT) unknowns. GNSS is
arguably the primary technology when it comes to position, navigation, and
timing (PNT) applications, such as intelligent transportation systems [3, 4,
5], critical infrastructures [6], environmental applications [7], or
agriculture applications [8] to name a few. However, the increasing dependence
on GNSS has created a concern about its potential vulnerabilities [9], mostly
related to GNSS vulnerability to simple jamming interferences [10, 11] or more
complex spoofing attacks [12].

It is therefore of crucial importance to develop interference monitoring
systems that can detect and localize the possible threats in an area [13, 14].
In this paper we are particularly interested in congested areas, where
malicious users are more prone to operate, where crowdsourced [15] data can be
leveraged in order to implement such monitoring. In particular, we foresee a
system where agents navigate an area with capabilities of transmitting the
measured signal power at the corresponding GNSS frequency bands. Recall that
GNSS signals are transmitted using a spread-spectrum system such that they are
below the receiver’s noise floor, as a consequence signal values with powers
larger than the noise floor would likely correspond to a jamming signal being
present in the vicinity. If a large number of agents convey this information
to a central unit, the data can be processed in order to localize the source
of jamming interferences and, ultimately, take corrective measures. Several
challenges of this approach are addressed in this paper, namely the need for
an accurate signal propagation model from which one can perform inference
from. The nominal pathloss model has severe disadvantages when it comes to
complex environments, given that accounting for reflections and obstructions
is not simple. To overcome this, we propose to augment the pathloss model with
a data-driven component (e.g. a neural network, or NN for short) that is able
to learn the position dependent parts of the propagation model that are
complex to accurately model. We refer to these models as augmented physics-
based models (APBMs), which were succesfully used in other contexts [16, 17,
18].

Section 2 presents the jammer localization problem, discussing the nominal
pathloss signal propagation model that is used as a measurement equation in
the proposed scheme. Under the nominal model, the Cramér-Rao bound (CRB) of
estimation accuracy is derived as well as the maximum likelihood estimator
(MLE). The proposed APBM augmentation to the pathloss is discussed in Section
3, where it is shown how we propose to jointly estimate the jammer’s location
and the parameters of the augmented model. Experiments are provided in Section
4 for $i)$ an open-sky scenario where the pathloss is an accurate model, which
is used to validate the proposed methods; and $ii)$ a more complex urban
scenario, where the received signal powers are generated using ray tracing
techniques, where the APBM shows its superiority to plain pathloss modeling.
Finally, Section 5 discusses the main conclusions and outlook of this work.

##  2 Pathloss-based estimation of jammer’s position

###  2.1 The pathloss measurement model

We propose a crowdsourced framework where the position of the jammer is
estimated from a sequence of $N$ observations of the jammer power at several
locations. In particular, the $n$-th observation $y_{n}$ (dBW) is related to
the vector of parameters we want to estimate through a generic function
$f(\mathbf{x}_{n};\bm{\theta})$, such that

$y_{n}=f(\mathbf{x}_{n};\bm{\theta})+\xi_{n}$ (1)

where $\mathbf{x}_{n}=({x}_{n}^{(1)},\dots,{x}_{n}^{(D)})^{\top}$ is the
location where the measurement took place,
$\bm{\theta}=(\theta_{1},\dots,\theta_{D})^{\top}$ is the jammer location in
$D$ dimensions (often 3), and $\xi_{n}$ is the measurement noise, which is
assumed to be additive and independent of $\mathbf{x}_{n}$ and $\bm{\theta}$.
This measurement model is particularly relevant in the context of GNSS
interference monitoring since useful signals are received below the noise
floor and large signal powers can be regarded as jamming signals.

A common model for received signal strength (RSS) observations is the pathloss
model [19]:

$f(\mathbf{x}_{n};\bm{\theta})=P_{0}-\gamma
10\log_{10}d(\mathbf{x}_{n},\bm{\theta})$ (2)

where $P_{0}$ (dBW) is the jammer power at the reference distance of 1 m, and
$d(\mathbf{x}_{n},\bm{\theta})$ is the distance between the $n$-th observer at
$\mathbf{x}_{n}$ and the jammer [20]. It is defined as

$d(\mathbf{x}_{n},\bm{\theta})=\|\mathbf{x}_{n}-\bm{\theta}\|=\sqrt{(\mathbf{x}_{n}-\bm{\theta})^{\top}(\mathbf{x}_{n}-\bm{\theta})}$
(3)

with $\|\cdot\|$ being the Euclidean norm.

###  2.2 The Cramér-Rao bound

The CRB of the $i$-th parameter is defined through the Fisher information
matrix (FIM) [21] as

$[\bm{I}(\bm{\theta})]_{ij}=-E\left[\frac{\partial^{2}\ln
p{(\mathbf{y};\bm{\theta})}}{\partial\theta_{i}\partial\theta_{j}}\right]$ (4)

where $p{(\mathbf{y};\bm{\theta})}$ is the _likelihood_ function of the
observed data and $\mathbf{y}=(y_{1},\dots,y_{N})^{\top}$ is the vector of
observations. As a result

$\operatorname{var}([\hat{\bm{\theta}}]_{i})\geq[\bm{I}^{-1}(\bm{\theta})]_{ii}$
(5)

gives the minimum variance attainable by an unbiased estimator
$\hat{\bm{\theta}}$ of the parameters $\bm{\theta}$.

In (2), the unpredictable shadowing effects are modeled by $\xi_{n}$ and are
experienced by measurements which have the same locations, but different
clutter on the propagation path. It has been shown that in nominal conditions
this measurement noise can be modeled with a log-normal distribution [22, 23],
thus resulting in $\xi_{n}\sim\mathcal{N}(0,\sigma)$ when values are expressed
in dB as in (2). As a result, the vector of observations $\mathbf{y}$ is
distributed according to

$\mathbf{y}|\bm{\theta}\sim\mathcal{N}(\bm{\mu}(\bm{\theta}),\mathbf{C}(\bm{\theta}))$
(6)

We can thus resort to a generalized CRB for the general Gaussian case [21].
When (6) holds, equation (4) becomes

$[\bm{I}(\bm{\theta})]_{ij}=\left[\frac{\partial\bm{\mu}(\bm{\theta})}{\partial\theta_{i}}\right]^{\top}\mathbf{C}^{-1}(\bm{\theta})\left[\frac{\partial\bm{\mu}(\bm{\theta})}{\partial\theta_{j}}\right]\\\
+\frac{1}{2}\operatorname{tr}\left[\mathbf{C}^{-1}(\bm{\theta})\frac{\partial\mathbf{C}(\bm{\theta})}{\partial\theta_{i}}\mathbf{C}^{-1}(\bm{\theta})\frac{\partial\mathbf{C}(\bm{\theta})}{\partial\theta_{j}}\right]\
\;,$ (7)

which under the assumption of independent measurements and constant variance
we get $\mathbf{C}(\bm{\theta})=\sigma I$, can be simplified to

$\displaystyle[\bm{I}(\bm{\theta})]_{ij}$
$\displaystyle=\left[\frac{\partial\bm{\mu}(\bm{\theta})}{\partial\theta_{i}}\right]^{\top}\frac{1}{\sigma^{2}}I\left[\frac{\partial\bm{\mu}(\bm{\theta})}{\partial\theta_{j}}\right]$
(8)
$\displaystyle=\frac{1}{\sigma^{2}}\sum_{n=1}^{N}\frac{\partial[\bm{\mu}(\bm{\theta})]_{n}}{\partial\theta_{i}}\frac{\partial[\bm{\mu}(\bm{\theta})]_{n}}{\partial\theta_{j}}$
(9) $\displaystyle=\frac{1}{\sigma^{2}}\sum_{n=1}^{N}\frac{\partial
f(\mathbf{x}_{n};\bm{\theta})}{\partial\theta_{i}}\frac{\partial
f(\mathbf{x}_{n};\bm{\theta})}{\partial\theta_{j}}\ .$ (10)

Using the particular expression of the pathloss model in (2), the FIM can be
further particularized to

$\displaystyle\frac{\partial
f(\mathbf{x}_{n};\bm{\theta})}{\partial\theta_{i}}$
$\displaystyle=-10\gamma\frac{\partial\log_{10}(d(\mathbf{x}_{n},\bm{\theta}))}{\partial\theta_{i}}$
(11)
$\displaystyle=\frac{-10\gamma}{d(\mathbf{x}_{n},\bm{\theta})\ln{10}}\frac{(2\theta_{i}-2{x}_{n}^{(i)})}{2d(\mathbf{x}_{n},\bm{\theta})}$
(12)
$\displaystyle=\frac{-10\gamma}{\ln{10}}\frac{(\theta_{i}-{x}_{n}^{(i)})}{d^{2}(\mathbf{x}_{n},\bm{\theta})}\
.$ (13)

The elements of the FIM then become

$[\bm{I}(\bm{\theta})]_{ij}=\frac{100\gamma^{2}}{\sigma^{2}(\ln{(10)})^{2}}\sum_{n=1}^{N}\frac{(\theta_{i}-{x}_{n}^{(i)})(\theta_{j}-{x}_{n}^{(j)})}{d^{4}(\mathbf{x}_{n},\bm{\theta})}$
(14)

resulting in

$\bm{I}(\bm{\theta})=\frac{100\gamma^{2}}{\sigma^{2}(\ln{(10)})^{2}}\sum_{n=1}^{N}\frac{1}{d^{4}(\mathbf{x}_{n},\bm{\theta})}(\bm{\theta}-\mathbf{x}_{n})^{\top}(\bm{\theta}-\mathbf{x}_{n})\
.$ (15)

The CRB for the $i$-th element of $\hat{\bm{\theta}}$ can be obtained by
replacing (15) into (5). For the 2-dimensional case we can derive a simple
analytic CRB for the estimator of the jammer’s position. Thus, the inverse of
$\bm{I}(\bm{\theta})$ can be easily computed leading to

$\operatorname{var}({\hat{\theta}_{1}})\geq\frac{\sigma^{2}(\ln{(10)})^{2}}{100\gamma^{2}}\frac{b}{ab-c^{2}}$
(16)
$\operatorname{var}({\hat{\theta}_{2}})\geq\frac{\sigma^{2}(\ln{(10)})^{2}}{100\gamma^{2}}\frac{a}{ab-c^{2}}\,\cdot$
(17)

where we defined

$a=\sum_{n=1}^{N}\frac{\
(\theta_{1}-{x}_{n}^{(1)})^{2}}{d^{4}(\mathbf{x}_{n},\bm{\theta})}$ (18)
$b=\sum_{n=1}^{N}\frac{\
(\theta_{2}-{x}_{n}^{(2)})^{2}}{d^{4}(\mathbf{x}_{n},\bm{\theta})}$ (19)
$c=\sum_{n=1}^{N}\frac{\
(\theta_{1}-{x}_{n}^{(1)})(\theta_{2}-{x}_{n}^{(2)})}{d^{4}(\mathbf{x}_{n},\bm{\theta})}\
.$ (20)

###  2.3 Maximum Likelihood Estimator

In this section we present the MLE estimator for the jammer’s position
$\bm{\theta}$, considering an arbitrary RSS function
$f(\mathbf{x}_{n};\bm{\theta})$. For this, let us define
$\mathbf{X}=\\{\mathbf{x}_{1},\dots,\mathbf{x}_{N}\\}$ and
$\mathcal{D}=\\{y_{n},\mathbf{x}_{n}\\}_{n=1}^{N}$ as a dataset composed of
$N$ i.i.d. pairs of loci $\mathbf{x}_{n}$ and RSS measurements $y_{n}$.
Assuming (6) to define the statistical characteristics of the model in (1),
the likelihood distribution can be written as:

$\displaystyle p(\mathbf{y}|$
$\displaystyle\mathbf{X},\bm{\theta})=\prod_{n=1}^{N}p(y_{n}|\mathbf{x}_{n},\bm{\theta})$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{N/2}}\exp{\left\\{-\frac{1}{2\sigma^{2}}\sum_{n=1}^{N}(y_{n}-f(\mathbf{x}_{n};\bm{\theta}))^{2}\right\\}}\
.$ (21)

The log-likelihood function is therefore

$\ln{p(\mathbf{y}|\mathbf{X},\bm{\theta})}=-\frac{N}{2}\ln{(2\pi\sigma^{2})}-\frac{1}{2\sigma^{2}}\sum_{n=1}^{N}(y_{n}-f(\mathbf{x}_{n};\bm{\theta}))^{2}$
(22)

The estimator for the jammer’s location, $\hat{\bm{\theta}}$, can be found by
maximizing the $\log$-likelihood as

$\hat{\bm{\theta}}_{\mathrm{MLE}}=\mathop{\arg\max}_{\bm{\theta}}\ln
p(\mathbf{y}|\mathbf{X},\bm{\theta})\;.$ (23)

One drawback of the pathloss model is that
$f(\mathbf{x}_{n};\bm{\theta})\to\infty$ when
$d(\mathbf{x}_{n},\bm{\theta})\to 0$. This generate singularities in the
likelihood function that need to be addressed when solving the optimization in
(23). To circumvent this issue we modified (2) to

$\bar{f}(\mathbf{x}_{n};\bm{\theta})=P_{0}-\gamma
10\log_{10}\\{\max(d(\mathbf{x}_{n},\bm{\theta}),d_{F})\\}$ (24)

where $d_{F}$ is the far-field distance, a limit that depends on the
characteristics of the transmitting antenna and above which the far-field
assumption that motivates the pathloss formulation in (23) is considered to
hold [24]. The effect of this choice is twofold. First, infinite values are
removed from the observation function and singularities are avoided. Secondly,
the “holes” visible in Figure 1 (left panel) are filled before the gradient
becomes too steep, towards the singularity points. In this way the
maximization of the log-likelihood can be pursued efficiently through
gradient-based methods. The resulting log-likelihood is shown in Figure 1
(right panel). This strategy is also effective when minimizing the cost
function described next.

Fig. 1: Log-likelihood with respect to different $\bm{\theta}$, measured under
an ideal pathloss propagation scenario with $f(\mathbf{x}_{n};\bm{\theta})$
(left panel) and $\bar{f}(\mathbf{x}_{n};\bm{\theta})$ (right panel).

##  3 Augmented physics-based estimation of jammer’s position

Although the pathloss model in (1) is well suited to describe the jamming RSS
field over open-spaces, it fails to capture the complexity of urban or indoor
scenarios where multiple reflections and attenuations might occur.
Additionally, our practical experience indicated that naive, in the sense that
no prior information regarding the jamming field is incorporated in the model,
NNs can only represent the jamming field accurately enough if humongous
amounts of data is readily available. For this reason, in this section we
propose an augmented physics-based (i.e., pathloss-based) neural network model
where the NN acts as a nonlinear correction term around the pathloss function.
Recently, APBMs have been used successfully in different scenarios [25, 16,
17, 18]. Here, we propose to augment the modified pathloss model with NNs
leading to an enhanced model:

$h(\mathbf{x}_{n};\bm{\theta},\bm{\phi})=\bar{f}(\mathbf{x}_{n};\bm{\theta})+g(\mathbf{x}_{n};\bm{\phi})$
(25)

where $\bm{\phi}\in\mathbb{R}^{M}$ is the vector of neural network’s
parameters. The measurement model in (1) can then be re-written as

$y_{n}=h(\mathbf{x}_{n};\bm{\theta},\bm{\phi})+\xi_{n}\;,$ (26)

and used to perform inference as described later in this section.

To limit the NN from overpowering the physics-based model, the authors in [17,
18] considered different regularizations employed to force the NN to act as a
correction term around the physics-based model, maintaining interpretability
of the overall model. Similarly, we consider a $\ell_{2}$ regularization to
control the NN contribution. Thus, we define the regularized cost function as

$\mathcal{C}(\mathcal{D},\bm{\theta},\bm{\phi})=\sum_{n=1}^{N}\|\mathbf{y}_{n}-h(\mathbf{x}_{n};\bm{\theta},\bm{\phi})\|^{2}_{2}+\beta\|\bm{\phi}\|^{2}_{2}$
(27)

where $\beta\in\mathbb{R}_{+}$ is a scalar controlling the regularization over
the NN’s parameters, such that when $\beta=0$ the NN dominates the APBM and
when $\beta\rightarrow\infty$ the NN contribution is reduced. Finally, we
define the optimization problem with respect to the jammer location
$\bm{\theta}$ and the NN’s parameters $\bm{\phi}$ as a empirical risk
minimization problem as:

$(\hat{\bm{\theta}},\hat{\bm{\phi}})=\mathop{\arg\min}_{\bm{\theta},\bm{\phi}}\mathbb{E}_{\hat{p}(\mathcal{D})}\bigg{\\{}\mathcal{C}(\mathcal{D},\bm{\theta},\bm{\phi})\bigg{\\}}$
(28)

where the expectation operator is taken with respect to the empirical data
distribution $\hat{p}(\mathcal{D})$.

##  4 Experiments

In this section we present our experimental results in two simulated
scenarios: an ideal scenario, where the jamming signal propagates according to
a nominal pathloss (PL) model, defined by $\gamma=2$ and described by (2); and
an urban scenario, subject to intense multipath and shadowing effects. The
urban scenario is modeled exploiting a ray tracing approach [26, 27],
computing multiple propagation paths using 3-D environment geometry and
electromagnetic analysis, including free-space loss and reflection losses up
to 4 reflections.

Fig. 2: Power measurements under nominal path loss propagation.

Fig. 3: Power measurements under ray tracing model propagation.

Particularly, we compare the localization performance of $(i)$ the MLE in (23)
under the PL model assumption; $(ii)$ the proposed APBM ($\beta=1$)
implemented as in (28); $(iii)$ the APBM ($\beta=1$) where the parameter
$P_{0}$ is also estimated jointly along $\bm{\theta}$ and $\bm{\phi}$, which
is of practical importance; $(iv)$ a PL-only solution for (28) where
$h(\mathbf{x}_{n};\bm{\theta},\bm{\phi})=\bar{f}(\mathbf{x}_{n};\bm{\theta})$;
and $(v)$ a NN-only solution for (28) where
$h(\mathbf{x}_{n};\bm{\theta},\bm{\phi})=g(\mathbf{x}_{n};\bm{\phi})$.
Additionally, we plot the CRB (under the PL model) for the sake of
benchmarking.

Both scenarios involve a static jammer, transmitting with $P_{0}=10$ dB. We
collected $10000$ observations over an area of $1$ km2, yielding, on average,
1 observation every $10\times 10$ m square. Since the signal power decreases
very fast with respect to the distance from the source, weak power
observations are not informative and buried in the noise floor. Only the $15$
most powerful observations were fed to the estimators. We repeated the data
collection for different measurement noise variance and therefore different
interference-to-noise ratio (INR), defined as
$\text{INR}=10\log_{10}P_{0}/\sigma^{2}\ .$ For each INR level,
$N_{\textrm{MC}}=100$ Monte Carlo realizations were simulated to collect
statistics about the estimators performances. To compare the different
approaches we consider the root-mean-square error (RMSE) computed over each
dimension of $\bm{\theta}$ as

$\text{RMSE}_{\theta_{i}}=\sqrt{\frac{1}{N_{\textrm{MC}}}\sum_{n=1}^{N_{\textrm{MC}}}(\theta_{i}-\hat{\theta}_{i,n})^{2}}$
(29)

The NN model is a feed-forward network with 2 hidden layers of 200 and 100
neurons and a hyperbolic tangent activation function. The PL model in $(ii)$
and $(iv)$ is initialized with the same information available to the MLE, i.e.
$P_{0}$ and $\gamma$. The models are trained over 200 epochs with a learning
rate of $0.4$, exploiting the Adam optimizer [28] and batch learning.

Fig. 4: RMSE in the PL propagation scenario.

Fig. 5: RMSE in the ray tracing propagation scenario.

In the PL scenario, addressed in Figure 4, the jamming signal propagates in
nominal conditions, according to (2), and $(i)$ and $(iv)$ are based on a
model that matches exactly this propagation rule. Not surprisingly, as the INR
increases, they are all approaching the CRB, which is also based on perfect
modeling. The use of the PL model component alone enables to closely follow
the MLE, even at low INR values, where the fixed implementation settings might
prevent such method to perform at its best. Indeed, it is worth stressing that
the common settings described at the beginning of this section are kept
constant for all the scenarios and methods. Moreover, the solution adopted to
circumvent the pathloss formula singularities (see Section 2.3) might affect
the learning process in a systematic way for a given set of measurements and
jammer locations. On the contrary, despite the simple propagation rule, the
NN-only solution is unable to effectively learn the model and it is not able
to improve further above 20 dB of INR. When the separate components are joined
again to realize the APBM we observe a low RMSE, very close to the MLE, also
around 10 dB where the complete system $(ii)$ is outperforming the PL-only
model $(iv)$. In these conditions, the NN component is probably helping the
APBM to overcome the learning flaws of the PL model. Conversely at very low
INR the high amount of measurement noise has a counter-productive effect on
the joint learning. The APBM can also work being $P_{0}$-blind $(iii)$, with
no knowledge about the jammer’s power and indeed converge to a correct
estimation (Figure 4). In this case the NN is effectively able to learn a
constant term, which compensate for the missing jammer information, providing
a fundamental aid. The low RMSE values attained by the MLE and PL model in
Figure 4 are possible thanks to a perfect knowledge of the propagation
environment. When it comes to a realistic scenario such as the one in Figure
3, few valid assumptions on the pathloss can be made. This results in the
jammer location estimation performances shown in Figure 5. The MLE and PL-only
solution are far form the attainable minimum, instead the APBM and NN-only
model are performing better. How close are their respective RMSE values
strictly depends on how far the actual propagation model is from an ideal
pathloss. The more mismatched is the PL model the closer the APBM performance
is to the NN, as in this case. It is worth noting that despite the large
difference between the propagation scenarios, the APBM can seamlessly adapt to
both without changes. Also, in this realistic scenario, the ability of the
APBM to perform a successful $P_{0}$-blind estimation is even more relevant,
since this information is not usually available.

##  5 Conclusion

In this work we proposed a jamming source localization strategy using
crowdsourced measurements. We investigated the use of an augmented physics-
based model (APBM) approach that models the complex propagation channel with a
physics-based pathloss model and a data-driven NN component, jointly
estimating the position of the jamming source and the NN model parameters. The
proposed method exhibited flexibility to characterize both ideal propagation
conditions and complex urban scenarios, attaining a performance comparable to
a perfect model-aware MLE in the first case, and better than a stand-alone NN
in the latter. Moreover, such performance can be reached also with no
information about the jamming source, thus without knowledge about the
transmission power, a parameter not available in realistic scenarios. Further
works will investigate the adaptability of the method in dynamic scenarios,
involving a moving and intermittent jammer, considering also continual
learning to address the changes in data distribution. To address privacy
concerns, we intend to connect this work with federated learning schemes.

## References

  * [1] D. Dardari, P. Closas, and P. M. Djurić, “Indoor tracking: Theory, methods, and technologies,” _IEEE Transactions on Vehicular Technology_ , vol. 64, no. 4, pp. 1263–1278, Apr. 2015.
  * [2] Y. J. Morton, F. van Diggelen, J. J. Spilker Jr, B. W. Parkinson, S. Lo, and G. Gao, _Position, navigation, and timing technologies in the 21st century: integrated satellite navigation, sensor systems, and civil applications_. John Wiley & Sons, 2021\.
  * [3] A. Yastrebova, M. Höyhtyä, S. Boumard, E. S. Lohan, and A. Ometov, “Positioning in the arctic region: State-of-the-art and future perspectives,” _IEEE Access_ , vol. 9, pp. 53 964–53 978, 2021.
  * [4] N. Williams, P. B. Darian, G. Wu, P. Closas, and M. Barth, “Impact of Positioning Uncertainty on Connected and Automated Vehicle Applications,” _SAE International Journal of Connected and Automated Vehicles_ , vol. 6, no. 12-06-02-0010, 2022.
  * [5] A. Minetto, A. Nardin, and F. Dovis, “GNSS-only collaborative positioning among connected vehicles,” in _Proceedings of the 1st ACM MobiHoc Workshop on Technologies, MOdels, and Protocols for Cooperative Connected Cars_ , ser. TOP-Cars ’19. Association for Computing Machinery, 2019, p. 37–42.
  * [6] R. T. Ioannides, T. Pany, and G. Gibbons, “Known vulnerabilities of global navigation satellite systems, status, and potential mitigation techniques,” _Proceedings of the IEEE_ , vol. 104, no. 6, pp. 1174–1194, 2016.
  * [7] R. Imam, M. Pini, G. Marucco, F. Dominici, and F. Dovis, “Data from GNSS-Based Passive Radar to Support Flood Monitoring Operations,” in _2019 International Conference on Localization and GNSS (ICL-GNSS)_ , 2019, pp. 1–7.
  * [8] W. Wang, N. Okati, I. Tanash, T. Riihonen, and E.-S. Lohan, “Location-based beamforming architecture for efficient farming applications with drones,” in _2019 International Conference on Localization and GNSS (ICL-GNSS)_ , 2019, pp. 1–6.
  * [9] M. G. Amin, P. Closas, A. Broumandan, and J. L. Volakis, “Vulnerabilities, threats, and authentication in satellite-based navigation systems [scanning the issue],” _Proceedings of the IEEE_ , vol. 104, no. 6, pp. 1169–1173, 2016.
  * [10] D. Borio, F. Dovis, H. Kuusniemi, and L. L. Presti, “Impact and detection of GNSS jammers on consumer grade satellite navigation receivers,” _Proceedings of the IEEE_ , vol. 104, no. 6, pp. 1233–1245, 2016.
  * [11] R. Morales-Ferre, P. Richter, E. Falletti, A. de la Fuente, and E. S. Lohan, “A survey on coping with intentional interference in satellite navigation for manned and unmanned aircraft,” _IEEE Communications Surveys Tutorials_ , vol. 22, no. 1, pp. 249–291, 2020.
  * [12] H. Sathaye, G. LaMountain, P. Closas, and A. Ranganathan, “SemperFi: A Spoofer Eliminating GPS Receiver for UAVs,” in _Proceedings of Network and Distributed System Security Symposium (NDSS 2022)_ , 2023.
  * [13] F. Dovis, _GNSS interference threats and countermeasures_. Artech House, 2015.
  * [14] S. Thombre, M. Z. H. Bhuiyan, P. Eliardsson, B. Gabrielsson, M. Pattinson, M. Dumville, D. Fryganiotis, S. Hill, V. Manikundalam, M. Pölöskey _et al._ , “GNSS threat monitoring and reporting: Past, present, and a proposed future,” _The Journal of Navigation_ , vol. 71, no. 3, pp. 513–529, 2018.
  * [15] E. Arias-de Reyna, P. Closas, D. Dardari, and P. M. Djuric, “Crowd-based learning of spatial fields for the internet of things: From harvesting of data to inference,” _IEEE Signal Processing Magazine_ , vol. 35, no. 5, pp. 130–139, 2018.
  * [16] A. Demirkaya, T. Imbiriba, K. Lockwood, S. Rampersad, E. Alhajjar, G. Guidoboni, Z. Danziger, and D. Erdoğmuş, “Cubature Kalman Filter Based Training of Hybrid Differential Equation Recurrent Neural Network Physiological Dynamic Models,” in _2021 43rd Annual International Conference of the IEEE Engineering in Medicine Biology Society (EMBC)_ , 2021, pp. 763–766.
  * [17] H. Li, R. A. Borsoi, T. Imbiriba, P. Closas, J. C. M. Bermudez, and D. Erdoğmuş, “Model-based deep autoencoder networks for nonlinear hyperspectral unmixing,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 19, pp. 1–5, 2022.
  * [18] T. Imbiriba, A. Demirkaya, J. Duník, O. Straka, D. Erdogmus, and P. Closas, “Hybrid neural network augmented physics-based models for nonlinear filtering,” in _2022 25th International Conference on Information Fusion (FUSION)_ , 2022, pp. 1–6.
  * [19] P. Wu, T. Imbiriba, G. LaMountain, J. Vilà-Valls, and P. Closas, “Wifi fingerprinting and tracking using neural networks,” in _Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2019)_ , 2019, pp. 2314–2324.
  * [20] T. S. Rappaport, _Wireless Communications: Principles and Practice, 2nd Edition_. Upper Saddle River, New Jersey: Prentice Hall, 2002, pp. 102–105.
  * [21] S. M. Kay, _Fundementals of Statistical Signal Processing: Estimation Theory_. Upper Saddle River, New Jersey: Prentice Hall, 1993.
  * [22] R. Bernhardt, “Macroscopic diversity in frequency reuse radio systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 5, no. 5, pp. 862–870, 1987.
  * [23] D. C. Cox, R. R. Murray, and A. W. Norris, “800-MHz attenuation measured in and around suburban houses,” _AT & T Bell Laboratories Technical Journal_, vol. 63, no. 6, pp. 921–954, 1984.
  * [24] T. S. Rappaport, _Wireless communications: principles and practice_. Prentice-Hall, 1996.
  * [25] I. Santos, J. J. Murillo-Fuentes, and P. M. Djurić, “Recursive Estimation of Dynamic RSS Fields Based on Crowdsourcing and Gaussian Processes,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 5, pp. 1152–1162, 2019.
  * [26] K. Schaubach, N. Davis, and T. Rappaport, “A ray tracing method for predicting path loss and delay spread in microcellular environments,” in _[1992 Proceedings] Vehicular Technology Society 42nd VTS Conference - Frontiers of Technology_ , 1992, pp. 932–935 vol.2.
  * [27] Z. Yun and M. F. Iskander, “Ray tracing for radio propagation modeling: Principles and applications,” _IEEE Access_ , vol. 3, pp. 1089–1100, 2015\.
  * [28] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014. [Online]. Available: https://arxiv.org/abs/1412.6980