#  Induced Quantized Spin Current in Vacuum

Chong-Sun Chu  Chun-Hei Leung  Department of Physics, National Tsing-Hua
University, Hsinchu 30013, Taiwan  Center for Theory and Computation, National
Tsing-Hua University, Hsinchu 30013, Taiwan

###### Abstract

We uncover a fundamental effect of the QED vacuum in an external
electromagnetic (EM) field. We show that the quantized vacuum of electrons is
spin polarized by the EM field and manifests as a vacuum spin current. An
experiment is proposed to measure the spin torque exerted by the spin current
by measuring the twisted angle of the director axis of a nematic liquid
crystal.

Introduction. — The quantum nature of vacuum is a fascinating place to look
for novel physical phenomena. For example, the large scale fluctuation of the
Universe is supposed to emerge from the quantum ripples of the near de Sitter
vacuum in the inflationary universe [1]. In high energy physics, the
metastability of the Higgs vacuum [2] has nontrivial cosmological consequences
and could provide a new observational window to particle physics well beyond
what collider experiments can achieve.

Recently, it was found that the vacuum fluctuations of boundary QED system in
the presence of an external magnetic field result in a magnetization current
near its boundary [3, 4]. This nontrivial electromagnetic response of the
vacuum is due to the electric charges carried by the virtual electrons and
positrons of the theory. As these quantum fluctuations also carry spin, it is
natural to ask if and how the quantum fluctuation of the spin degrees of
freedom of the vacuum would manifest in observation.

In this Letter, we show that due to the renormalization effect of the spin-
orbit coupling of the electrons, the quantized vacuum can become spin
polarized in the presence of an applied EM field and results in a spin
current; see Eq. (25). The vector spin current is orthogonal to the applied
electric field and the observed spin polarization. We also discuss briefly how
it may be observed experimentally using a setup of nematic liquid crystal.

Definition of spin current. — The transportation of the spin of electrons is
described by a spin current. Classically, the spin current is given by the
3-tensor $J^{ij}=nv^{i}s^{j}$, where $v^{i}$ is the velocity of the electron,
$s^{j}$ is the spin, $n$ is the particle number density, and $i,j=1,2,3$
denotes the spatial directions. Quantum mechanically, the spin polarization of
electrons is represented by the Pauli matrices $\sigma^{i}$ in the rest frame
of the electron and a definition of the spin current:

$J^{\mu i}=\frac{1}{2}(v^{\mu}\sigma^{i}+\sigma^{i}v^{\mu}),\quad
J^{0i}=\sigma^{i}$ (1)

in terms of a symmetrization of the velocity and the spin has been proposed by
Rashba [5]. Although this works well for a nonrelativistic system, the
generalization to the relativistic case requires the use of a relativistic
spin operator. Such an operator was constructed long ago by Bargmann and
Wigner [6] for the free Dirac theory, and it reads

$T_{\rm BW}^{i}=\beta\Sigma^{i}+\gamma^{5}p^{i}/m,\quad T_{\rm
BW}^{0}={\bf\Sigma}\cdot{\bf p}/m,$ (2)

where $\beta=\gamma^{0}$, $\alpha^{i}=\beta\gamma^{i}$,
$\gamma^{5}=i\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}$ is the chirality matrix
and $\Sigma^{i}=(i/2)\varepsilon^{ijk}\gamma^{j}\gamma^{k}$. We will adopt the
particle physics convention of metric signature $(+,-,-,-)$ in this Letter.
The Bergmann-Wigner spin operator $T^{\mu}$ is a Lorentz 4-vector [7], it acts
as the generator of the little group of the Poincare group. That this is so is
because it is indeed equivalent to the Pauli-Lubanski operator

$W^{\mu}=-\frac{1}{2m}\varepsilon^{\mu\alpha\beta\nu}J_{\alpha\beta}p_{\nu}.$
(3)

This can be seen by noticing that the spin angular momentum generator is given
by $J_{\alpha\beta}=(i/2)[\gamma_{\alpha},\gamma_{\beta}]$ in the spinor
representation. It is then easy to see that Eq. (3) is equal to the spin
operator Eq. (2) on shell.

In this Letter, we are interested in the renormalization phenomena of the QED
vacuum in the presence of an external electromagnetic field $A^{\mu}$. As the
velocity is represented by $v^{\mu}=\beta\gamma^{\mu}$ in the Dirac theory, we
propose the following definition of spin current in quantum field theory:

$J_{S}^{\mu\alpha}(x):=\bar{\psi}(x)S^{\mu\alpha}\psi(x),$ (4)

where

$S^{\mu\alpha}:=\frac{1}{2}[\gamma^{\mu}T^{\alpha}(x)+T^{\alpha}(x)\gamma^{\mu}]$
(5)

and $T^{\alpha}$ is the covariant Bergmann-Wigner spin operator in a
background electromagnetic field:

$T^{i}=\beta\Sigma^{i}+\gamma^{5}\frac{\pi^{i}}{m},\quad
T^{0}=\frac{1}{m}{\bf\Sigma}\cdot{\bf\pi}.$ (6)

Here $\pi^{i}=p^{i}-eA^{i}$ is the covariant momentum and $e>0$ is the
magnitude of the electric charge. Similar definitions [8, 9] have been
considered before. Our definition is better justified since we have adopted a
symmetrization prescription such that Eq. (4) reduces to the one (Eq. (1)) of
Rashba in the nonrelativistic limit. A further justification of our definition
(4) can be obtained by noticing that the current Eq. (4) satisfies the
conservation law

$\partial_{\mu}J^{\mu\alpha}=\frac{e}{m}{\cal S}_{\nu}F^{\nu\alpha},$ (7)

with a source term. Here ${\cal S}^{\nu}$ is the density defined by

${\cal S}^{i}:=\psi^{\dagger}\Sigma^{i}\psi,\quad{\cal
S}^{0}:=\psi^{\dagger}\gamma^{5}\psi.$ (8)

We note in passing that $J^{\mu\alpha}$ is not a Noether current, but Eq. (7)
is obtained from the fermion equation of motion. For the spatial directions
$\alpha=i$, Eq. (7) gives explicitly

$\frac{\partial\rho^{i}}{\partial
t}+\nabla\cdot{\pmb{J}}^{i}=\frac{e}{m}\Big{(}({\bf{\cal
S}}\times\pmb{B})^{i}+{\cal S}^{0}E^{i}\Big{)}.$ (9)

Here $\rho^{i}:=J^{0i}$ and $\pmb{J}^{i}:={\bf e}_{k}J^{ki}$. It is clear that
$\rho^{i}$ gives the spin polarization density. As a result, $\pmb{J}^{i}$
does admit the correct interpretation as a current density for spin
polarization in the $i$th direction. We also note that
$-\nabla\cdot{\pmb{J}}^{i}$ and the EM terms on the right-hand side of Eq. (9)
can be interpreted as the spin torque from the matters and the external EM
fields. Classically in a vacuum, ${\cal S}^{\mu}$ vanishes and the spin
current respects the continuity equation

$\frac{\partial\rho^{i}}{\partial t}+\nabla\cdot\pmb{J}^{i}=0.$ (10)

We will show in this Letter that, in the presence of a background
electromagnetic field, the spin current $\pmb{J}^{i}$ and the spin density
$\rho^{i}$ becomes nonzero due to the polarization effect of the
electromagnetic field on the quantum fluctuations of the vacuum. Nevertheless
the conservation law (10) is still satisfied at the quantum mechanical level.

We end this section with a couple of remarks. (i) Properly speaking, the
current Eq. (4) measures the flow of spin polarization and should be called a
spin polarization current. We will follow the common practice of the
literature, e.g., Ref. [5], and refer to it as the spin current. To get the
spin momentum current, one needs to multiply Eq. (4) with $\hbar/2$ of the
spin angular momentum of each fermion. (ii) We note that in the
nonrelativistic limit, the magnetic dipole coupling and the spin-orbit
coupling term in the Dirac Hamiltonian

$H_{E,B}=-\frac{e}{2m}{\bf\Sigma}\cdot\pmb{B}-\frac{e}{4m^{2}}{\bf\Sigma}\cdot{\bf
E}\times{\bf p}$ (11)

can be written in the form

$H_{E,B}=-\frac{e}{8m}\varepsilon_{\mu\nu\lambda\sigma}J^{\mu\nu}F^{\lambda\sigma},$
(12)

where $J^{\mu\nu}$ here is the spin current (1) in the nonrelativistic limit.
The form (12) shows clearly that the spin current, at least in the
nonrelativistic limit, couples to the EM field strengths instead of the EM
potentials like the electric current $J^{\mu}=-e\bar{\psi}\gamma^{\mu}\psi$.
(iii) Equation (12) also shows that a certain amount of energy is needed to
generate the stated vacuum spin current and this is supplied by the external
power source that maintains the background EM field configuration. (iv)
Finally, we remark that for a quantum field theory in curved space, there is a
current $J_{\mu}^{ab}$ that couples to the spin connection
$\omega_{\mu}^{ab}$. Here $a,b=0,1,2,3$ refer to the frame indices. This
current couples to gravity and is also sometimes referred to as a spin current
in the respective community. However we emphasize that this is different from
the spin current we introduced in this Letter.

Vacuum expectation of spin current. — In quantum field theory, the spin
current (4) is a composite operator which needs to be renormalized. We are
interested in the vacuum expectation value (vev) of the spin current in a
background electromagnetic field $A_{\mu}$. This can be computed in
perturbation theory as

$\langle J^{\mu\alpha}(x)\rangle_{A}=-ie\int d^{4}y\langle
J^{\mu\alpha}(x)J^{\beta}(y)\rangle A_{\beta}(y)+O(A^{2}),$ (13)

or, it can be written in the momentum space as

$\langle
J^{\mu\alpha}(q)\rangle_{A}=eT^{\mu\alpha\beta}(q)A_{\beta}(q)+O(A^{2}),$ (14)

where $T^{\mu\alpha\beta}(q)$ is the Green’s function

$T^{\mu\alpha\beta}(q):=-i\int d^{4}xe^{iqx}\,\langle
J^{\mu\alpha}(x)J^{\beta}(0)\rangle.$ (15)

At 1-loop, $T^{\mu\alpha\beta}(q)$ is given by

$T^{\mu\alpha\beta}(q)=-i\int\frac{d^{4}p}{(2\pi)^{4}}(-1)\mbox{tr}\left(S^{\mu\alpha}\frac{i}{{p/}\;-m}\gamma^{\beta}\frac{i}{{p/}\;+q/\hskip
1.49994pt-m}\right)$ (16)

where $S^{\mu\alpha}$ is given by Eq. (5). The trace of the gamma matrices can
be simplified and we obtain

$T^{\mu 0\beta}=0,\qquad T^{\mu
i\beta}=-4m\varepsilon^{i\mu\delta\beta}q_{\delta}\,I(q),\quad i=1,2,3,$ (17)

where $I(q)$ is the momentum function defined by

$I(q):=\int\frac{d^{4}p}{(2\pi)^{4}}\frac{1}{p^{2}-m^{2}}\frac{1}{(p+q)^{2}-m^{2}}.$
(18)

Note that $I(q)$ and hence the Green’s function is logarithmic divergent. This
is due to the singular product of local quantum fields in the expression (4)
for the spin current operator. In order to give a proper definition of this
composite operator, one need to isolate the divergent terms with a
regularization scheme and subtract it away with counterterms in the Lagrangian
[10]. Regularizing Eq. (18), we obtain

$I(q)=\frac{i}{16\pi^{2}}\Big{(}\log\frac{\Lambda^{2}}{\mu^{2}}-1+\log\frac{\mu^{2}}{m^{2}}-h\Big{(}\frac{q^{2}}{m^{2}}\Big{)}+\cdots\Big{)}$
(19)

for a momentum cutoff regularization, and

$I(q)=\frac{i}{16\pi^{2}}\Big{(}\frac{2}{\epsilon}-\gamma_{E}+\log\frac{4\pi\mu^{2}}{m^{2}}-h\Big{(}\frac{q^{2}}{m^{2}}\Big{)}+\cdots\Big{)}$
(20)

for dimensional regularization to $d=4-\epsilon$ dimension. Here $\cdots$
denotes terms of order $O(1/\Lambda^{2})$ or $O(\epsilon)$, and $h(x)$ is the
function

$h(x):=\int_{0}^{1}d\xi\,\ln[1+\xi(1-\xi)x]$ (21)

with $h(0)=0$. We can now subtract away the divergence with a counterterm and
obtain

$\langle J_{R}^{\mu
i}(x)\rangle_{A}=-\frac{emc}{8\pi^{2}\hbar}\varepsilon^{i\mu\delta\beta}\Big{(}\log\frac{m^{2}}{\mu^{2}}+a_{J}+h\Big{(}\frac{\square}{m^{2}}\Big{)}\Big{)}F_{\delta\beta}(x)$
(22)

and $\langle J_{R}^{\mu 0}(x)\rangle_{A}=0$. Here $\mu$ is the standard RG
scale of the QFT and $a_{J}$ is an independent arbitrary constant that is due
to the arbitrariness in the choice of the finite part of the counterterm for
the composite operator $J^{\mu i}$.

To uniquely fix the finite part (i.e., fixing $a_{J}$) and hence the
definition of the renormalized spin current, a normalization condition is
required. Physically, if the electron mass is sent to infinity, then the
quantum loop effects are completely suppressed and the renormalized spin
current mush vanish. In a theory with a cutoff $\Lambda$, the fermion loop
effects are now suppressed as $m$ approaches this scale. As a result, we have
the decoupling condition

$\lim_{m\to\Lambda}\langle J_{R}^{\mu i}\rangle_{A}=O(\frac{1}{\Lambda}).$
(23)

This serves as a natural normalization condition for the renormalized spin
current. However extra care is needed for QED where there is a Landau pole and
the UV cutoff cannot be taken to exceed that. In fact, the positiveness of the
beta function $\beta(\alpha)=2\alpha^{2}/3\pi$ ($\alpha=e^{2}/4\pi\hbar c$ is
the fine structure constant) give rises to the RG flow for the coupling

$\frac{1}{\alpha(\mu)}-\frac{1}{\alpha(\mu_{0})}=-\frac{2}{3\pi}\log\frac{\mu}{\mu_{0}}.$
(24)

This implies the presence of a cutoff scale $\Lambda_{L}=me^{3\pi/2\alpha}$,
the Landau scale, where the bare coupling $\alpha(\Lambda_{L})$ becomes
infinite. The Landau scale represents the highest possible UV cutoff one may
utilize in the renormalization program of QED. Taking $\Lambda=\Lambda_{L}$ in
Eq. (23), the renormalization constant $a_{J}=-\log\Lambda_{L}^{2}/\mu^{2}$ is
fixed up to $O(1/\Lambda_{L})$ terms, which we will ignore. Restoring the
units of $c/\hbar^{2}$ [11], we finally obtain

$\langle J_{R}^{\mu
i}(x)\rangle_{A}=\frac{emc}{8\pi^{2}\hbar^{2}}\Big{(}\frac{3\pi}{\alpha}-h\Big{(}\frac{\square}{m^{2}}\Big{)}\Big{)}\varepsilon^{i\mu\delta\beta}F_{\delta\beta}.$
(25)

Note that the current Eq. (25) is conserved due to the Bianchi identity of the
EM field. For a classical electromagnetic background in vacuum, it is $\square
F_{\alpha\beta}=0$ and we have

$\langle J_{R}^{\mu
i}(x)\rangle_{A}=\frac{3emc}{8\pi\hbar^{2}\alpha}\varepsilon^{i\mu\delta\beta}F_{\delta\beta}$
(26)

or, in terms of components explicitly

$\displaystyle\langle J^{ji}_{R}\rangle_{A}$ $\displaystyle=$
$\displaystyle-\frac{3emc}{4\pi\hbar^{2}\alpha}\varepsilon^{jik}E^{k},$ (27a)
$\displaystyle\langle J^{0i}_{R}\rangle_{A}$ $\displaystyle=$
$\displaystyle\frac{3emc}{4\pi\hbar^{2}\alpha}B^{i}$ (27b)

The results (25), (26), (27) are the main results of this Letter. For the rest
of the Letter, we will be focusing on the case (26) of a classical background.
However for generality we present the result (25) to cover the situation where
the background EM fields go beyond the Maxwell description; for example, if
quantum nonlinear corrections of QED is included, or if the classical EM
fields are coupled to other background fields (e.g., an axion background).

A couple of remarks are in order. (i) Note that our result is independent of
the adopted regularization scheme. Apart from the cutoff regularization and
the dimensional regularization, one can also do the heat-kernel regularization
and obtain the same intermediate result (22), and the final result (25) after
imposing the normalization condition. (ii) It is interesting that the
normalization condition (23) fixes $a_{J}$ in terms of the renormalization
scale $\mu$. In the end, $\mu$ remains free in the QFT and it is remarkable
that the renormalized spin current is independent of it. (iii) We remark that
the renormalization condition (23) is essentially nonperturbative in nature as
the resulting prediction (25) scales inversely with $\alpha$ and is not
smoothly connected with the free theory result.

It is instructive to compare our effect with the famous Schwinger effect [12,
13], which refers to the production of electron-positron pairs under the
influence of an applied electric field. The Schwinger effect is
nonperturbative and requires a strong electric field stronger than the
critical field strength $E\gtrsim E_{\rm critical}=m^{2}c^{3}/(e\hbar)$ in
order to produce an observable amount of particle pairs. Unlike the Schwinger
effect, the spin current predicted in this Letter is a nonperturbative
consequence of the nontrivial spin polarization of the vacuum. There is no
real production of particles involved. Our result is also different from the
analysis of Ref. [9] where the Schwinger-effect-produced electric current is
acted on by a second electric field to produce a spin current. In this case,
the spin current is generated in a similar manner as in an ordinary material
sample except that the source electric current has a nonperturbative origin
and so very small in magnitude.

In the above we have considered pure QED. In a more realistic setting where
QED is embedded as the low energy part of a consistent high energy theory,
e.g., a grand unified theory or string theory when gravity is included, the
Landau pole would be replaced by the corresponding GUT or Planck scale. In
this case, the 1-loop results (25) will have the factor $3\pi/\alpha$ replaced
by $3\pi(1/\alpha-1/\alpha_{G})$ where $\alpha_{G}$ is the fine structure
coupling at the unification scale. It is interesting that the vacuum
expectation value of the spin current actually provides a probe to the UV
physics. The value of $1/\alpha_{G}$ is model dependent. For example, for the
MSSM GUT, one has [14] $1/\alpha_{G}\simeq 24.3$ and $M_{G}\simeq 2\times
10^{16}$ GeV. In any case, $1/\alpha_{G}$ is expected to be small compared to
the observed $1/\alpha$ at the electron mass scale. In the following, we will
continue to analyze the result (25) for pure QED, but keeping in mind the
overall magnitude of the spin current may be different from Eq. (25) by a
small fraction.

Physical picture. — It may appear strange that the switching on of an electric
field in vacuum can produce an observable spin current (27). However the
physical origin of the spin current can be easily understood in turns of the
spin-orbit coupling of the quantum fluctuation of the vacuum. To see this, let
us consider the nonrelativistic expansion of Dirac’s Hamiltonian up to
$O(p^{4})$ where a spin-orbit coupling term arises,
$H_{\text{SO}}=-[(e\hbar)/(4m^{2}c^{2})]{\bf\Sigma}\cdot{\bf E}\times{\bf p}$.
It is well known that the spin-orbit coupling term allows impurity in a
material to scatter the electrons in a spin-dependent way (skew scattering)
and generates a spin current [15, 16]. The spin-orbit coupling also give rises
to a side jump [17, 18] described by the additional velocity

${\bf v}_{\text{SO}}=\frac{\partial H_{\text{SO}}}{\partial{\bf
p}}=-\lambda\,{\bf\Sigma}\times{\bf
E},\quad\lambda:=\frac{e\hbar}{4m^{2}c^{2}},$ (28)

which produces a spin-dependent shift to the trajectory of the electrons and
contributes to the spin current.

In vacuum, there are no charge carriers to start with, but virtual pairs of
electrons and positrons can be created and live for a short time before they
annihilate back. In standard QFT without background field, such vacuum
polarization processes give rise to the running of couplings and scaling of
dimensions in the theory. However, new observable effects may result in the
presence of a background field. Without an electric field, the created
electron-positron pair will move away from each other with opposite velocities
$\pm u_{0}^{j}$ due to momentum conservation. When an $E$ field is turned on,
the particles acquire an additional velocity (28) in accordance with their
spin states. The velocities of the virtual $e^{\pm}$ particles are given by
$u^{j}=\pm u_{0}^{j}+v_{\pm}^{j}$, where
$v_{\pm}^{j}=-\lambda\varepsilon^{jkl}\sigma_{\pm}^{k}E^{l}$ and
$\sigma^{i}_{\pm}$ are the spin operators for the $e^{+}$ and $e^{-}$
respectively. Take a measurement of the spin in the $i$th direction, the
expectation value of the $j$th component of the spin current $\langle
J^{ji}\rangle=\langle n\sigma_{+}^{i}(u_{0}^{j}+v_{+}^{j})\rangle+\langle
n\sigma_{-}^{i}(-u_{0}^{j}+v_{-}^{j})\rangle$ gives

$\langle J^{ji}\rangle=-2n\lambda\varepsilon^{jik}E_{k},$ (29)

where $n$ is the number density of the virtual pairs. In deriving Eq. (29), we
have used $\langle\sigma^{i}_{\pm}\rangle=0$ and the correlations
$\langle\sigma_{\pm}^{i}\sigma_{\pm}^{j}\rangle=\delta^{ij}{\bf 1}$ for the
vacuum. As a result, the free field velocity parts in Eq. (29) make no
contribution, while the anomalous velocity part of the virtual $e^{\pm}$
contribute equally to the spin current. With the estimate that there is one
virtual pair within each volume of Compton wavelength, $n\sim(mc/\hbar)^{3}$,
we recover the result (27a) up to a numerical coefficient (including the sign)
of order 10. It is clear from this quantum mechanical argument that the
induced spin current is a result of the nontrivial spin-spin correlation
imprinted on the vacuum due to the applied electric field. This is similar to
the induced current [3, 4] and the induced Fermi condensate [19, 20] which are
due to nontrivial magnetization and condensation of the quantized vacuum as a
result of external fields.

Proposed experiment. — The induced spin current may be observed by measuring
the torque exerted by the spin current on a probe placed in the vacuum.
Consider an infinitesimal volume element $\delta V=\delta x\delta y\delta z$
in the interior of a probe placed under the influence of an external EM field.
The spin momentum torque acting on $\delta V$ is
$\tau^{i}=-\hbar/2\int_{\delta
V}d^{3}x\,\nabla\cdot\langle{\pmb{J}}^{i}\rangle$. Using the result (27) for a
classical EM field in vacuum, we have

$\tau^{i}=\frac{3emc}{8\pi\hbar\alpha}\overline{\dot{B}^{i}}\delta V,$ (30)

where $\overline{\dot{B}^{i}}:=\delta V^{-1}\int_{\delta V}\partial
B^{i}/\partial t$ is the average rate of change of the magnetic field over the
volume $\delta V$ We will be using the SI units from now on and hence a factor
of $c$ appears in Eq. (30). For a $B$ field pointing in the $z$ direction
described by a wave of the form

$B_{z}=B_{0}g\Big{(}t-\frac{x}{c}\Big{)},$ (31)

where $g$ is as in Fig. 1,

Figure 1: Waveform of the $B$ field.

we obtain

$\tau^{z}=\frac{3emc}{8\pi\hbar\alpha}B_{0}f\delta V,$ (32)

where $f:=1/T_{1}-1/T_{2}$. For an external field with $B_{0}=100$ G, $f=1000$
Hz, the quantum spin torque per unit volume is

$\tau^{z}/\delta V=6.5\times 10^{-5}\,{\rm Nm}{}^{-2}.$ (33)

Not all of the vacuum torque is transferred to the probe. Physically the
angular momentum density $\mathcal{J}^{i}_{0}=(\hbar/2)\rho^{i}$ acquired by
the vacuum generates a vacuum magnetization
$\pmb{M}_{0}=(e/m)\mathbfcal{J}_{0}$ and this corresponds to quantum addition
to the $B$ field, $\Delta\pmb{B}=\mu_{0}\pmb{M}_{0}$. For a material probe
with magnetic susceptibility $\chi$, the $B$ field generates a probe
magnetization $\pmb{M}=\chi\pmb{M}_{0}/(1+\chi)$. These are just the Barnett
effect and the Einstein-de Haas effect for the interplay between angular
momentum and magnetization [21]. This implies the probe receives a torque
$\tau^{i}_{p}$:

$\tau_{p}^{i}=\eta\tau^{i},$ (34)

where $\eta=4\pi\chi/(1+\chi)\approx\chi$ as $|\chi|\ll 1$. The spin torque on
the probe is thus tiny. Nevertheless, the twisting effect may be observable by
using a liquid crystal which is known to be exceptional in sensitivity for
torque measurement. In fact, recently the Casimir torque exerted on the
surface of a liquid crystal [22] has been observed successfully [23]. We
propose here a similar setup to observe the spin torque arising from the spin
current.

Figure 2:  The director axis of a liquid crystal probe is twisted by the
induced spin torque exerted by the vacuum.

Consider a layer of nematic liquid crystal, e.g., 4-cyano-40-pentylbiphenyl
(5CB), placed in a vacuum with one side ($z=0$) anchored at a glass interface
and the rest of the liquid crystal ($0<z\leq L_{z}$) sits in the vacuum
cavity, influenced by the external EM field; See Fig. 2. We place the liquid
crystal such that its director axis is always in the $xy$ plane, i.e.,
$\pmb{n}(z)=[\cos\theta(z),\sin\theta(z),0]$ where $\theta(z)$ describes the
orientation of the liquid crystal molecules with respect to the $x$ axis.
Because of the action of the spin torque (32), the molecules try to orient
themselves correspondingly. An equilibrium configuration is attained when this
torque is balanced by the restoring elastic torque of the liquid crystal.

Suppose there is no bend or splay of the liquid crystal, then only the twist
contribute to distorting energy density $u_{d}=(k/2)[(\partial\theta/\partial
z)]^{2}$, where $k$ is the twist elastic constant. For example, $k=3.6$ pN for
the crystal 5CB. The elastic energy stored in $\delta V$ is thus

$E_{d}[\theta(z)]=A\int_{z}^{z+\delta
z}dz\frac{k}{2}\Big{(}\frac{\partial\theta}{\partial z}\Big{)}^{2},$ (35)

where $A=\delta x\delta y$ is the area of the element. This give rises to a
restoring torque

$\tau_{d}(z)=\frac{\delta
E_{d}}{\delta\theta(z)}=-k\frac{\partial^{2}\theta}{\partial z^{2}}\delta V.$
(36)

This is in analogy with the Newton’s second law $F=m\ddot{x}$ for the inertia.
The equilibrium is attained when Eq. (34) is balanced out by the restoring
torque [24]. Note that the volume factor cancels out. This gives

$\theta(z)=\frac{\beta}{2}z^{2},\quad\beta:=\frac{3\pi\chi}{2\alpha}\frac{mf_{B}f}{k\lambda_{C}},$
(37)

where $f_{B}:=eB_{0}/(2\pi m)$ is the cyclotron frequency and $\lambda_{C}$ is
the Compton wavelength. In deriving Eq. (37), we have taken the boundary
conditions $\theta=\partial_{z}\theta=0$ at the glass contact $z=0$.

To enhance the detection of the twisted angle, a thicker liquid crystal layer
and a stronger magnetic field is preferred. However, both of these are limited
by the properties of the liquid crystal. In order for the liquid crystal to be
able to register the torque, the response time of the liquid crystal should be
smaller than $T_{1}$ and $T_{2}$. This means $f$ is upper bounded by the
response frequency $f_{r}$: $f<f_{r}$. However, the response frequency is
limited by the thickness as generally the response time of the liquid crystal
increases with the thickness, linearly or quadratically depending on the
voltage [25]. As for the EM field, the field strengths cannot be too strong as
otherwise the crystal may be driven into a Freederick transition [25]. Nematic
liquid crystals are typically diamagnetic with $\chi$ of the order of
$10^{-5}$ [25]. Commercial nematic liquid crystal has a $f_{r}$ in the range
of 100 Hz and thickness of the order of 10 $\mu$m. As an estimate, consider a
setup with a 100 G $B$ field and a nematic liquid crystal with an elastic
constant $k$ = 3.6 pN, thickness $L_{z}$ = 1 mm, a response frequency of 1
kHz. The total twist angle accumulated over the thickness of the probe is

$\Delta\theta=0.005^{\circ}\cdot\frac{\chi}{10^{-5}}\cdot\frac{3.6\text{
pN}}{k}\cdot\frac{B_{0}}{100\text{ G}}\cdot\frac{f}{1\text{
kHz}}\cdot\left(\frac{L_{z}}{1\text{ mm}}\right)^{2}.$ (38)

This looks feasible. It will be interesting to perform an experiment to make
observation of the spin current predicted in this Letter. The quantum spin
current could have a wide range of applications, from novel theoretical
properties of physical system, e.g. dark energy, to practical effects on the
workings of micro-machined device.

###### Acknowledgements.

We acknowledge support of this work by the Grant No. 107-2119-M-007-014-MY3 of
the Ministry of Science and Technology of Taiwan.

## References

  * [1] A. R. Liddle and D. H. Lyth, Cosmological Inflation and Large-Scale Structure (Cambridge University Press, Cambridge, England, 2000).
  * [2] F. Bezrukov, J. Rubio and M. Shaposhnikov, Living beyond the edge: Higgs inflation and vacuum metastability, Phys. Rev. D 92, 083512 (2015).
  * [3] C. S. Chu and R. X. Miao, Weyl Anomaly Induced Current in Boundary Quantum Field Theories, Phys. Rev. Lett. 121, 251602 (2018).
  * [4] C. S. Chu and R. X. Miao, Anomalous Transport in Holographic Boundary Conformal Field Theories, J. High Energy Phys. 07 (2018) 005.
  * [5] E. I. Rashba, Spin currents in thermodynamic equilibrium: The challenge of discerning transport currents, Phys. Rev. B 68, 241315(R) (2003).
  * [6] V. Bargmann and E. P. Wigner, Group Theoretical Discussion of Relativistic Wave Equations, Proc. Natl. Acad. Sci. U.S.A. 34, 211 (1948).
  * [7] See, for example, D. M. Fradkin and R. H. Good, Jr., Electron Polarization Operators, Rev. Mod. Phys. 33, 343 (1961) for a review of the construction and properties of the Bergmann-Wigner spin operator.
  * [8] A. Vernes, B. L. Gyorffy, and P. Weinberger, Spin currents, spin-transfer torque, and spin-Hall effects in relativistic quantum mechanics, Phys. Rev. B 76, 012408 (2007).
  * [9] X. G. Huang, M. Matsuo and H. Taya, Spontaneous generation of spin current from the vacuum by strong electric fields, Prog. Theor. Exp. Phys. (2019), 113B02.
  * [10] J.C. Collins, Renormalization, (Cambridge University 448 Press, Cambridge, England, 1986).
  * [11] $c/\hbar^{2}=c\times 1/\hbar^{4}\times\hbar^{2}$ where $c$ comes from the velocity $v^{\mu}=c\alpha^{\mu}$, $1/\hbar^{4}$ comes from the momentum integration $d^{4}p/(2\pi\hbar)^{4}$ and $\hbar^{2}$ comes from the propagator $i\hbar/({p/}\;-m)$.
  * [12] J. Schwinger, On Gauge Invariance and Vacuum Polarization, Phys. Rev. 82, 664 (1951).
  * [13] J. Schwinger, The Theory of Quantized Fields. V, Phys. Rev. 93, 615 (1954); The Theory of Quantized Fields. VI, Phys. Rev. 94, 1362 (1954).
  * [14] See Sec. 94 “Grand Unified Theories” of P. A. Zyla et al. (Particle Data Group), Review of Particle Physics, Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * [15] J. Smit, The spontaneous hall effect in ferromagnetics I, Physica (Amsterdam) 21, 877 (1955);
  * [16] J. Smit, The spontaneous hall effect in ferromagnetics I, Physica (Amsterdam) 24, 39 (1958).
  * [17] L. Berger, Side-Jump Mechanism for the Hall Effect of Ferromagnets, Phys. Rev. B 2, 4559 (1970);
  * [18] L. Berger, Application of the Side-Jump Model to the Hall Effect and Nernst Effect in Ferromagnets, Phys. Rev. B 5, 1862 (1972).
  * [19] C. S. Chu and R. X. Miao, Fermion condensation induced by the Weyl anomaly, Phys. Rev. D 102, 046011 (2020).
  * [20] C. S. Chu and R. X. Miao, Weyl Anomaly induced Fermi Condensation and holography, J. High Energy Phys. 08 (2020) 134.
  * [21] L. D. Landau, L. P. Pitaevskii, E. M. Lifshitz, Electrodynamics of Continuous Media, 2nd ed. (Butterworth-Heinemann, Washington, DC, 1984).
  * [22] D. A. T. Somers and J. N. Munday, Rotation of a liquid crystal by the Casimir torque, Phys. Rev. A 91, 032520 (2015).
  * [23] D. A. T. Somers, J. L. Garrett, K. J. Palm, and J. N. Munday, Measurement of the Casimir torque, Nature (London) 564, 386 (2018).
  * [24] Our analysis is different from that of Ref. [22]. In Ref. [22], there is no external torque acting on the bulk of the liquid crystal and the equilibrium configuration is obtained by balancing the Casimir torque with a restoring elastic torque defined at the surface of the liquid crystal. As there is no torque in the bulk, the equation of motion (36) then give rises to a $\theta$ linear in $z$, which is also obtained in Ref. [22] with a different argument.
  * [25] I.C. Khoo, Liquid Crystals, 2nd ed. (John Wiley & Sons, New York, 2007).


# Pre-oligomerisation stochastic dynamics of prions driven by water molecules

Mairembam Kelvin Singh1<EMAIL_ADDRESS>R.K. Brojen Singh2
<EMAIL_ADDRESS>Moirangthem Shubhakanta Singh1<EMAIL_ADDRESS>1Department of Physics, Manipur University, Canchipur-795003, Imphal, Manipur,
India.
2School of Computational $\&$ Integrative Sciences, Jawaharlal Nehru
University, New Delhi-110067, India.

###### Abstract

Prions are proteinaceous infectious particles that cause neurodegenerative
diseases in humans and animals. The complex nature of prions, with respect to
their conformations and aggregations, has been an important area of research
for quite some time. Here, we develop a model of prion dynamics prior to the
formation of oligomers and subsequent development of prion diseases within a
stochastic framework, based on the analytical Master Equation and Gillespie’s
Stochastic Simulation Algorithm (SSA). The results that we obtain shows that
solvent water molecules act as driving agents in the dynamics of prion
aggregation. Further, it is found that aggregated and non-aggregated proteins
tend to co-exist in an equilibrium state, depending upon the reaction rate
constants. These results may provide a theoretical and qualitative contexts of
possible therapeutic strategies in the treatment of prion diseases.

Keywords: Prions; Protein aggregation; Master Equation; Stochastic Simulation.

##  IIntroduction

Prions are proteinaceous infectious particles that cause neurodegenerative
diseases in humans and animals. Notable diseases caused by prions are the
Bovine spongiform encephalopathy (BSE), scrapie of sheep and Creutzfeldt-Jakob
disease (CJD) of humans. Prions lack nucleic acid and the prion diseases
modify the normal cellular prion protein abbreviated as PrP${}^{\text{C}}$
into the modified form abbreviated as PrP${}^{\text{Sc}}$ Prusiner . This
conversion takes place when they have similar amino acid sequences Prusiner ;
Cohen . PrP${}^{\text{C}}$ contains mainly $\alpha$-helix isoform and is
protease-sensitive while PrP${}^{\text{Sc}}$ is mainly composed of
$\beta$-sheet isoform and is protease-resistant Prusiner ; Srivastava ;
Horwich .

The spread of prion diseases from one species to another is a stochastic
process, where the dynamics of the prions population is governed by
probabilistic equation(s). This process also features the formation of
intermediates Cohen  and a long incubation period from the infected time to
the onset of symptoms Prusiner . Studies have also stated that the diseases
are caused by protein misfolding and aggregation, when monomeric protein
chains diffuse within the same molecule Srivastava ; Lapidus . A unique
characteristic of prion diseases is that the aggregated prions promote further
aggregation of normal proteins, thus making it possible for transmission of
the diseases Horwich .

Therapeutic control strategies like using a particular molecule to stabilise
the structure of PrP${}^{\text{C}}$ or modifying the action of the
intermediates formed during the system reactions, were already suggested
Prusiner . These molecules must be able to prevent aggregation to control the
transmission of diseases Lapidus . However, one of the difficulties that
arises during the control strategies is that prions occur in different
“strains”, producing different patterns and complex behaviors Horwich ;
Lemarre . The lack of nucleic acid also makes the prions harder to remove by
ultraviolet irradiation Mobley . So, studying the complex behavior of prion
aggregation is crucial for developing a sustainable control strategy Lapidus ;
Chong .

When an ensemble of monomeric unfolded prion conformations is subjected to a
large number of solvent water molecules, oligomerisation or the formation of
oligomers takes place. Two aggregated prions can also form a dimerised protein
(dimer) through molecular diffusion, followed by the production of oligomers,
which may be toxic to cells Lapidus ; Stephens ; Chong . These oligomers,
then, eventually form amyloid fibrils Chong ; Eisenberg . The water molecules
play a big role in controlling the dynamics of the prions. The hydration shell
present in prion proteins controls the movement and function of the proteins
and the removal of which can drive protein folding Stephens ; Verma ; Jamal ;
Uversky .

Various mathematical models and computational simulations to study the
aggregation of prions have been done in recent years. In Matthaeus , brain
infection by prion diseases has been studied by developing a mathematical
model based on the reaction-diffusion equations considering the interaction of
prions with transport on networks. Qualitative studies have shown the ability
of the prions to co-exist in a stable state influenced by initial conditions
applied to the mathematical model developed Lemarre . From a stochastic point
of view, a stochastic cellular automata model has been developed in Mobley
based on the probabilistic nature of interaction and diffusion. Liu and Kang
Liu  studied the protein aggregation by developing the chemical master
equation and applying the lognormal moment closure method with the Stochastic
Simulation Algorithm (SSA) by Gillespie Gillespie  confirming its accuracy.

In this paper, we also studied the prion aggregation dynamics using the Master
Equation approach and solve it by the generating function technique McQuarrie
; Gardiner  to obtain the Poissonian probability distribution of dimers,
formed during the process of aggregation and development of oligomers. We also
used Gillespie’s SSA Gillespie  to run various simulations based on different
reaction rates and initial conditions of the solvent water molecules and non-
aggregated prion proteins.

We found that the asymptotic behavior of dimers follow a Poisson process and
the mean population of dimers also saturate to a value equal to the ratio of
the mean production rate to the mean degradation rate. This saturation is also
confirmed by the SSA simulation. Further, we found that certain reaction rates
control the co-stability and co-existence Lemarre  of aggregated and non-
aggregated prions, without producing significant amount of oligomers in the
process.

The paper is organised as follows. Section II describes the model reactions
depicting the various processes involved during the aggregation. Section III
briefly describes the theoretical context of the analytical and computational
methods used. Section IV shows the various analytical and computational
simulation results that we obtain and its subsequent analysis for the
significance in providing possible qualitative control strategies.

##  IIModel Description

We consider an ensemble of well-mixed unfolded prion monomers subjected to a
large number of solvent water molecules $(H)$. Let $X$ denote the non-
aggregated prions and $X^{*}$ denote the aggregated prions. If two $X^{*}$
interact through molecular diffusion, a dimer $D$ is formed. This is followed
by the formation of an oligomer $O$ or production of $X$ and $X^{*}$. These
various processes Srivastava ; Stephens  are then represented in the form of a
system of reactions as shown below:

$\displaystyle\begin{split}X+H&\stackrel{{\scriptstyle\rm
k_{1}}}{{\longrightarrow}}X^{*}+H\\\ X^{*}+H&\stackrel{{\scriptstyle\rm
k_{2}}}{{\longrightarrow}}X+H\\\ 2X^{*}&\stackrel{{\scriptstyle\rm
k_{3}}}{{\longrightarrow}}D\\\ D&\stackrel{{\scriptstyle\rm
k_{4}}}{{\longrightarrow}}X+X^{*}\\\ D&\stackrel{{\scriptstyle\rm
k_{5}}}{{\longrightarrow}}O+H\\\ O&\stackrel{{\scriptstyle\rm
k_{6}}}{{\longrightarrow}}D\end{split}$ (1)

The first two reactions denote the inter-conversions of aggregated $(X^{*})$
and non-aggregated $(X)$ prions through the help of the water molecules $(H)$
with reaction rate constants $k_{1}$ and $k_{2}$. The third reaction denotes
the formation of a dimer $(D)$ from two aggregated prions $(X^{*})$ with
reaction rate constant $k_{3}$. The fourth reaction denotes the conversion of
a dimer $(D)$ into an aggregated $(X^{*})$ and non-aggregated $(X)$ prion with
reaction rate constant $k_{4}$. The fifth reaction denotes the formation of an
oligomer $(O)$ and water $(H)$ from a dimer $(D)$ with reaction rate constant
$k_{5}$. The last reaction denotes the conversion of an oligomer $(O)$ to a
dimer $(D)$ with reaction rate constant $k_{6}$.

Throughout our study, we assume that the system is a well-mixed and well-
stirred system, with the volume taken as $V=1$ so that the stochastic and
deterministic rate constants are equal.

##  IIIMethods

The system of reactions in (1) is studied and analysed through two approaches,
the analytical master equation approach McQuarrie ; Gardiner  and the
computational simulation using Gillespie’s SSA Gillespie . These approaches
are briefly described in the next subsections.

###  III.1 Analytical Method

The dynamics of a particle in a stochastic system can be described by the
Chapman-Kolmogorov (CK) equation Kolmogorov ; Gardiner  where the state
transitions obey Markov process. The CK equation in integral form for a
stochastic particle transitioning from a state $(y^{\prime},t^{\prime})$ to a
state $(y,t)$ is given by,

$\displaystyle\frac{\partial}{\partial t}P(y,t|y^{\prime},t^{\prime})$
$\displaystyle=$ $\displaystyle-\sum_{i}\frac{\partial}{\partial
u_{i}}\left[A_{i}(u,t)P(u,t|y^{\prime},t^{\prime})\right]$ (2)
$\displaystyle+\frac{1}{2}\sum_{i,j}\frac{\partial^{2}}{\partial u_{i}\partial
u_{j}}\left[B_{ij}(u,t)P(u,t|y^{\prime},t^{\prime})\right]$
$\displaystyle+\int dy[W(u|y,t)P(y,t|y^{\prime},t^{\prime})$
$\displaystyle-W(y|u,t)P(u,t|y^{\prime},t^{\prime})]$

where $\left\\{W\right\\}$ are Wiener functions of transition probabilities of
the states, $A_{i}$ and $B_{ij}$ are the functions depending on the properties
of the system and $P(y,t|y^{\prime},t^{\prime})$ is the jump probability. The
master equation , which is the discrete form of the CK equation, derived from
the integral part of the CK equation (2), is given by Kolmogorov ; Gardiner ;
Kampen ,

$\displaystyle\frac{\partial}{\partial t}P(y,t|y^{\prime},t^{\prime})$
$\displaystyle=$ $\displaystyle\int dy[W(u|y,t)P(y,t|y^{\prime},t^{\prime})$
(3) $\displaystyle-W(y|u,t)P(u,t|y^{\prime},t^{\prime})]$

The master equation can be considered as the time evolution of the probability
of a particle in a system obeying Markov process Singh . We can solve this
master equation (3) by the generating function (GF) technique McQuarrie . The
$n$-dimensional GF of the probability distribution function
$P(y_{1},y_{2},...;t)$ is given by,

$\displaystyle
G(s_{1},s_{2},...,s_{n})=\sum_{y_{1},y_{2},...,y_{n}}\prod_{i=1}^{n}s_{i}^{y_{i}}P(y_{1},y_{2},...,y_{n};t)$
(4)

Multiplying the master equation (3) by
$\displaystyle\sum_{y_{1},y_{2},...,y_{n}}\prod_{i=1}^{n}s_{i}^{y_{i}}$ and
using (4), we can transform the master equation (3) to a spatio-temporal PDE
of the GF with the boundary condition, $\displaystyle
G(s_{1},s_{2},...,s_{n};0)=\prod_{i=1}^{n}s_{i}^{N_{i}^{[0]}}$, where
$\\{N_{i}^{[0]}\\};i=1,2,...,n$ is the set of initial populations of the
respective variables $y_{1},y_{2},...y_{n}$. After putting back the solution
of GF to (4) and equating the coefficients, we can obtain the solution
$P(y_{1},y_{2},...,y_{n};t)$ of the master equation (3) Singh .

###  III.2 Computational Simulation

For computational simulation, we use the stochastic simulation algorithm (SSA)
by Gillespie Gillespie , which is based on the assumption that the time
evolution of a stochastic particle is a discrete process obeying Markov
property. The population of the particle increases or decreases by integral
amounts according to the reactions occurring in the system.

Consider a system of $n$ molecular species reacting in $m$ different ways. The
joint reaction probability density function for a random reaction $r$ and
reaction time $T$ Gillespie  is then given as,

$\displaystyle
P(r,T)=\left\\{\begin{array}[]{ll}a_{r}\exp(-a_{0}T),&~{}\text{if}~{}0\leq
T<\infty\\\ &\text{and}~{}r=1,...,m\\\ 0,&\text{otherwise}\end{array}\right.$
(8)

The propensity function of the $r$th reaction is given by,

$\displaystyle a_{r}=h_{r}k_{r}~{}~{}\forall~{}r\in
m;~{}~{}~{}~{}a_{0}=\sum_{r=1}^{m}a_{r}=\sum_{r=1}^{m}h_{r}k_{r}$ (9)

As $P(r)$ and $P(T)$ are independent of each other, we can rewrite (8) as
Gillespie ,

$\displaystyle\begin{split}P(r,T)&=P(r)P(T)\\\
\text{where}~{}~{}P(T)&=a_{0}\exp(-a_{0}T)\\\
\text{and}~{}~{}P(r)&=\frac{a_{r}}{a_{0}}\end{split}$ (10)

From (10), we can generate two uniform random numbers $R_{1}$ and $R_{2}$
using a unit interval uniform random number generator for $P(r)$ and $P(T)$
respectively for simulation,

$\displaystyle\begin{split}T&=\left(\frac{1}{a_{0}}\right)\ln\left(\frac{1}{R_{1}}\right)\\\
\sum_{\nu=1}^{r-1}a_{\nu}&<R_{2}a_{0}\leq\sum_{\nu=1}^{r}a_{\nu}\end{split}$
(11)

The algorithm of the SSA is as follows:

  1. 1.

Initialise the values of the $m$ reaction constants $k_{1},...,k_{m}$ and the
$n$ initial population of the species involved. The time variable $t$ and the
reaction counter $n^{\prime}$ are both set to zero and the unit interval
uniform random number generator is also imported.

  2. 2.

Calculate the $m$ propensities $a_{1}=h_{1}k_{1},...,a_{m}=h_{m}k_{m}$ and
$a_{0}$ using (9).

  3. 3.

According to the unit interval uniform random number generator imported
earlier, two random numbers $R_{1}$ and $R_{2}$ are generated, $r$ and $T$ are
calculated, using (11).

  4. 4.

Increase $t$ by $T$, and change the population of the species according to the
particular reaction.

  5. 5.

Go back to step 2, recalculate $a_{\nu}$, $a_{0}$ and then proceed.

  6. 6.

Stop the loop when either $t$ or $n^{\prime}$ reaches the desired value, or if
$a_{0}$ becomes zero.

Figure 1: Plots of the PDF against the dimer population for different values
of time $t=10,50,100$ units. (A)-(C): These are plots corresponding to the
variation of $k_{3}$. As $k_{3}$ increases, the peak value decreases and the
distribution shifts towards the larger dimer population. Also, as time
increases, the distribution moves towards the larger dimer population but
saturates to a fixed pattern and the peak value corresponding to each $k_{3}$
decreases. (D)-(F): These are plots corresponding to the variation of $k_{6}$.
We can see a similar behavior as in the case of variation of $k_{3}$. Figure
2: (A) Fano factor of prion aggregation dynamics as a function of time. It
saturates to the value 1, implying that the dynamics is a Poisson process. (B)
Mean population of dimer as a function of time. The mean population saturates
to a value equal to $\beta$.

##  IVResults and Discussion

Based on the two methods used, the results are divided into two parts, an
analytical one and a simulation result. The analytical result gives the
qualitative analysis of the aggregation dynamics in terms of the probability
distribution function of the dimers $(D)$ formed during the prion aggregation.
On the other hand, the simulation study considers other species such as the
aggregated prions $(X^{*})$, non-aggregated prions $(X)$ and oligomers $(O)$
apart from the dimers $(D)$. So, the simulation study gives a more broader
picture of the time evolution of the different molecular species considered in
the system.

###  IV.1 Analytical Results

To derive the master equation of the dimers $(D)$ and study their dynamics, we
can ignore the first two reactions of (1) as they do not involve dimers $(D)$.
Considering the state transitions based on the system of reactions in (1) and
using the principle of detailed balance, the probability distribution function
of the dimer population $(D)$ at a particular time $t+\Delta t$, provided that
the previous state occurred in time $t$, is given by Singh ,

$\displaystyle P(D;t+\Delta t)$ $\displaystyle=$
$\displaystyle\frac{X^{*}\left(X^{*}-1\right)}{2}k_{3}P(D-1;t)\Delta t$
$\displaystyle+\left(D+1\right)k_{4}P(D+1;t)\Delta t$
$\displaystyle+\left(D+1\right)k_{5}P(D+1;t)\Delta t$
$\displaystyle+Ok_{6}P(D-1;t)\Delta t$ $\displaystyle+\left[1-\omega\Delta
t\right]P(D;t)$

where
$\displaystyle\omega=\frac{X^{*}\left(X^{*}-1\right)}{2}k_{3}+Dk_{4}+Dk_{5}+Ok_{6}$.

Applying the limit as $\Delta t\rightarrow 0$, we can obtain the master
equation as,

$\displaystyle\frac{\partial P(D;t)}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{X^{*}\left(X^{*}-1\right)}{2}k_{3}P(D-1;t)$ (12)
$\displaystyle+\left(D+1\right)k_{4}P(D+1;t)$
$\displaystyle+\left(D+1\right)k_{5}P(D+1;t)$ $\displaystyle+Ok_{6}P(D-1;t)$
$\displaystyle-\omega P(D;t)$

Now, using the generating function, $\displaystyle
G(s,t)=\sum_{D}s^{D}P(D;t)$, with the initial condition, $\displaystyle
G(s,0)=s^{M}$, where $M$ is the initial dimer population, we obtain the
solution of the master equation (12) as,

$\displaystyle P(D;t)$ $\displaystyle=$ $\displaystyle
e^{\beta\left(e^{-\alpha t}-1\right)}\sum_{i}^{M}C_{i}\frac{1}{(D-i)!}$ (13)
$\displaystyle e^{-i\alpha t}\beta^{D-i}\left(1-e^{-\alpha t}\right)^{M+D-2i}$

where $\alpha=k_{4}+k_{5}$ and
$\displaystyle\beta=\frac{X^{*}(X^{*}-1)k_{3}+2Ok_{6}}{2\alpha}$, the ratio of
the mean production rate to the degradation rate of dimers.

The probability distribution function (PDF) obtained in (13) is plotted with
respect to the dimer population in Figure 1. Upon increase of dimer production
from aggregated prions, i.e., increase of $k_{3}$, the peak value decreases
and the distribution shifts towards the larger dimer population. Also, as time
increases, the distribution moves towards the larger dimer population but
saturates to a fixed pattern, which can be seen by the close similarities of
Figures 1(B) and 1(C). The peak value corresponding to each $k_{3}$ also
decreases as time increases. We can also see the same behavior if we increase
the conversion of oligomer into dimer, i.e., increase $k_{6}$.

To analyse the distribution pattern, we can simplify the expression (13) with
the approximation $M-i\approx M$ and $D-i\approx D$ as,

$\displaystyle P(D;t)=(1-\gamma)^{M}e^{\gamma}\textit{Pois}(\delta)$ (14)

where $\displaystyle\gamma=e^{-\alpha t}$, $\delta=\beta(1-\gamma)$ and
Pois$(\delta)$ denotes the Poisson distribution function in $\delta$. This
deduction proves that the prion aggregation dynamics is a Poisson process.

####

To study the influence of noise in the prion aggregation dynamics, we
calculated the Fano factor Fano , using the relation,
$\displaystyle\mathcal{F}=\frac{\sigma^{2}}{\langle D\rangle}$, where
$\sigma^{2}$ is the standard deviation and $\langle D\rangle$ is the mean
population of dimers.

The mean of the dimer population is determined from the GF using the relation
$\displaystyle\langle D\rangle=\left.\frac{\partial G(s,t)}{\partial
s}\right|_{s=1}$. It is obtained as,

$\displaystyle\langle D\rangle=(M-\beta)\gamma+\beta$ (15)

The asymptotic values of the mean dimer population can be easily determined
as, $\displaystyle\lim_{t\rightarrow 0}\langle D\rangle=M$, the initial dimer
population and $\displaystyle\lim_{t\rightarrow\infty}\langle D\rangle=\beta$,
the ratio of the mean production rate to the degradation rate of dimers.

The standard deviation $\sigma^{2}$ can be determined using the relation,

$\displaystyle\sigma^{2}=\langle D^{2}\rangle-\langle D\rangle^{2}$ (16)

where $\displaystyle\langle
D^{2}\rangle=\left.\frac{\partial^{2}G(s,t)}{\partial
s^{2}}\right|_{s=1}+\langle D\rangle$. From (15) and (16), the Fano factor is
determined as,

$\displaystyle\mathcal{F}=\frac{\left(M\gamma+\beta\right)\left(1-\gamma\right)}{\left(M-\beta\right)\gamma+\beta}$
(17)

The asymptotic values of the Fano factor is also determined as,
$\displaystyle\lim_{t\rightarrow 0}\mathcal{F}=0$ and
$\displaystyle\lim_{t\rightarrow\infty}\mathcal{F}=1$. These values also show
that the prion aggregation dynamics is a Poisson process Rajdl ; Chanu . Thus,
there is no specific role of noise in the prion aggregation dynamics. The
asymptotic behavior of the Fano factor and the mean dimer population are also
shown in Figure 2.

Figure 3: The aggregated and non-aggregated prions approach towards a stable
value at a very quick time and their populations cluster together showing
signs of co-existence towards an equilibrium state (left column). The dynamics
of dimer population remain similar to the analytical result and the final
population also approaches towards the common stable value of the aggregated
and non-aggregated prions (middle column). There is very little to almost no
oligomerisation (right column).

###  IV.2 Simulation Results

Simulations using SSA are done to the reaction system (1) in different ways to
analyse the aggregation dynamics from different perspectives. For the first
simulation, the initial population of water molecules are varied. The values
of the reaction rate constants used are: $k_{1}=0.1$, $k_{2}=0.1$,
$k_{3}=0.01$, $k_{4}=0.05$, $k_{5}=0.01$ and $k_{6}=0.6$. The initial
populations of the different species involved are set as, water,
$H=(10^{4},10^{5},10^{6})$; non-aggregated proteins, $X_{in}=1000$; aggregated
proteins, $X^{*}_{in}=0$; dimer, $D_{in}=0$ and oligomer, $O_{in}=0$. In this
simulation, we kept $k_{1}>k_{3}$, i.e., the rate of formation of aggregated
prion proteins is more than the formation of dimers. The aggregated and non-
aggregated prions approach towards a stable value at a very quick time and
their populations cluster together showing signs of co-existence towards an
equilibrium state. The dynamics of dimer population remain similar to what we
obtained in the analytical result. The final population of dimers also
approaches towards the common stable value of the aggregated and non-
aggregated prions. There is little to almost no oligomerisation. As $H$
increases by a factor of 10, the total number of reactions for the fixed
amount of time increases by a factor of 10. The total number of reactions for
the respective initial water populations are of the order of
$(10^{5},10^{6},10^{7})$. This means that the population of water molecules
determines the total reaction rate of the system. Most of the reactions are
mainly the first two reactions, which also signify that water controls
oligomerisation process during the aggregation. These results are shown in
Figure 3.

Figure 4: (A) and (B): The populations of aggregated and non-aggregated
proteins no longer approach towards a common stable state. The population of
aggregated prions oscillate around their stable state, the non-aggregated
prions seem to decline towards the direction of aggregated prions indicating
the possibility of aggregated and non-aggregated prions approaching co-
stability at a much later time. (C): The population of aggregated prions
exceeds that of the descending population of non-aggregated prions but
eventually decreases to achieve the co-stability, sooner than the first two
cases.

Next, we keep the population of water molecules fixed at $H=10^{5}$ and
observe the aggregation dynamics at different variations of the rate constants
at the same range of time. For $k_{1}<k_{3}$ and $k_{1}=k_{3}$, the
populations of aggregated and non-aggregated prion proteins no longer approach
towards a common stable state (Figure 4A and 4B). However, the distinguishing
feature is that although the population of aggregated prions oscillate around
their stable state, the non-aggregated prions seem to decline towards the
direction of aggregated prions. This indicates the possibility of aggregated
and non-aggregated prions approach co-stability at a much later time. For the
case of $k_{2}<k_{1}$, i.e., the formation of aggregated prion proteins from
the non-aggregated prion proteins occurring at a faster rate than the reverse
process, the population of aggregated prions exceeds that of the descending
population of non-aggregated prions but eventually decreases to achieve the
co-stability, sooner than the first two cases (Figure 4C). In all three cases,
the total number of reactions is of the order of $10^{5}$.

Figure 5: (A)-(C): The aggregated and non-aggregated prions approach together
towards a stable value and their populations cluster together showing signs of
co-existence towards an equilibrium state. (D)-(F): Oligomer formation starts
and increases as $k_{5}$ increases.

Now, the production rate of oligomers $(k_{5})$ is increased with
$k_{1}>k_{3}$ and $H=10^{5}$. The dynamics of aggregated and non-aggregated
prion proteins remain similar and as expected, oligomer formation starts and
increases as $k_{5}$ increases (Figure 5). The total number of reaction is of
the order of $10^{6}$.

Figure 6: (A)-(C): The aggregated and non-aggregated prions approach together
towards a stable value and their populations cluster together showing signs of
co-existence towards an equilibrium state. (D)-(F): The population of dimers
saturate at a quicker time as the initial non-aggregated protein population
decreases. For $X_{in}=100$, the dimer population increases discretely and
there are instances where the population remains constant for a certain period
of time.

Finally, the initial population of non-aggregated proteins is decreased from
1000 to 500 and then to 100. The aggregated and non-aggregated prion dynamics
remain similar again. However, the dynamics of dimers begin to change. They
saturate at a quicker time as the initial non-aggregated protein population
decreases. Particularly, for $X_{in}=100$, the dimer population increases in
discrete steps and there are instances where the population remains constant
for a certain period of time (Figure 6). The total number of reactions also
decreases from the order of $10^{6}$ to $10^{5}$ as $X_{in}$ decreases. This
means that, like the water molecules, the initial population of non-aggregated
proteins controls the total reaction rate of the system.

##  VConclusion

From the study of the aggregation dynamics of prion proteins, we are able to
understand some qualitative behavior of the aggregation before oligomerisation
sets in. The stochastic analysis shows the sensitive dependence of aggregation
on the initial populations of the involved species. Specifically, we have
found that the water molecules and the initial population of non-aggregated
prions control the rate of total reactions of the system. The process of
oligomerisation only starts when the rate of oligomer formation is increased.
Inhibiting the rate of oligomerisation is a way to control the consequential
stages of disease transmission. In all the different observations made, the
presence of water controls the entire reaction system. This is evident from
the fact that out of all the total number of reactions, huge number comes from
the inter-conversion of aggregated to non-aggregated prion proteins through
water molecules and vice-versa. The rate of certain reactions also determine
the state of co-stability and co-existence between the aggregated and non-
aggregated prion proteins.

The analytical method is limited to the dynamics of dimers only. However, it
produces an important result that the dimer dynamics, and hence the
aggregation dynamics is a Poisson process, where noise or fluctuations have no
significant impact. The asymptotic results showing Poissonian dynamics also
describes the random nature of the dynamics.

On the other hand, the computational simulations by SSA enable us to
understand the aggregation dynamics, by considering the important parameters
like the reaction rates and the initial conditions. However, few features like
time delay and larger system size are neglected to simplify the study.
Although the aggregation occurs at a very quick time of the order of
nanoseconds, there may be certain time lag between the reactions. For this, we
need to use the delay algorithms Bratsun ; Barrio ; Cai ; Anderson  to make
realistic visualisations. Moreover, we are able to show the qualitative
picture of how the aggregation dynamics behave according to changing
parameters.

Acknowledgements
MKS is a Junior Research Fellow (JRF) under the National Fellowship for
Scheduled Castes Students (NFSC) scheme implemented by the National Scheduled
Castes Finance and Development Corporation (NSFDC) under the Department of
Social Justice & Empowerment, Ministry of Social Justice & Empowerment,
Government of India with UGC ref. no. 201610025460. The simulations were done
through Python 3 programs written using the Spyder environment of Anaconda
distribution. Simulations that require large computational tasks were done in
HPE ProLiant DL380 Gen10 Plus installed in the Department of Physics, Manipur
University.

Author Contributions:
RKBS and MSS conceptualized the work. MKS did the analytical, computational
work and preparation of the figures. MKS, MSS and RKBS wrote the manuscript,
read, analyzed the results, and approved the final manuscript.

Additional Information
Competing interests:
The authors declare no competing interests.

## References

  * (1) Prusiner, S.B. (1998) Prions. Proc. Natl. Acad. Sci. 95(23), 13363-13383.
  * (2) Cohen, F.E., Pan, K.M., Huang, Z., Baldwin, M., Fletterick, R.J. & Prusiner, S.B. (1994) Structural clues to prion replication. Science, 264(5158), 530-531.
  * (3) Srivastava, K.R. & Lapidus, L.J. (2017) Prion protein dynamics before aggregation. Proc. Natl. Acad. Sci., 114(14), 3572-3577.
  * (4) Horwich, A.L. & Weissman, J.S. (1997) Deadly conformations—protein misfolding in prion disease. Cell, 89(4), 499-510.
  * (5) Lapidus, L.J. (2013) Understanding protein aggregation from the view of monomer dynamics. Mol. BioSyst., 9(1), 29-35.
  * (6) Lemarre, P., Pujo-Menjouet, L. & Sindi, S.S. (2019) Generalizing a mathematical model of prion aggregation allows strain coexistence and co-stability by including a novel misfolded species. J. Math. Biol., 78, 465-495.
  * (7) Mobley, D.L., Cox, D.L., Singh, R.R., Kulkarni, R.V. & Slepoy, A. (2003) Simulations of Oligomeric Intermediates in Prion Diseases. Biophys. J., 85(4), 2213-2223.
  * (8) Chong, S.H. & Ham, S. (2015) Distinct role of hydration water in protein misfolding and aggregation revealed by fluctuating thermodynamics analysis. Acc. Chem. Res., 48(4), 956-965.
  * (9) Stephens, A.D. & Schierle, G.S.K. (2019) The role of water in amyloid aggregation kinetics. Curr. Opin. Struct., 58, 115-123.
  * (10) Eisenberg, D. & Jucker, M. (2012) The amyloid state of proteins in human diseases. Cell, 148(6), 1188-1203.
  * (11) Verma, P.K., Rakshit, S., Mitra, R.K. & Pal, S.K. (2011) Role of hydration on the functionality of a proteolytic enzyme $\alpha$-chymotrypsin under crowded environment. Biochimie, 93(9), 1424-1433.
  * (12) Jamal, S., Kumari, A., Singh, A., Goyal, S. & Grover, A. (2017) Conformational ensembles of $\alpha$-synuclein derived peptide with different osmolytes from temperature replica exchange sampling. Front. Neurosci., 11, 684.
  * (13) Uversky, V.N., Gillespie, J.R. & Fink, A. L. (2000) Why are “natively unfolded” proteins unstructured under physiologic conditions?. Proteins Struct. Funct. Genet., 41(3), 415-427.
  * (14) Matthaeus, F. (2009) The spread of prion diseases in the brain—models of reaction and transport on networks. J. Biol. Syst., 17(04), 623-641.
  * (15) Liu, R.N. & Kang, Y.M. (2020) Stochastic master equation for early protein aggregation in the transthyretin amyloid disease. Scientific Reports, 10(1), 12437.
  * (16) Gillespie, D.T. (1977) Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81(25), 2340-2361.
  * (17) McQuarrie, D.A. (1967) Stochastic approach to chemical kinetics. J. Appl. Probab., 4(3), 413-478.
  * (18) Gardiner, C.W. (1985). Handbook of stochastic methods. Vol. 3. Berlin: Springer.
  * (19) Kolmogorov, A.N. (1938) On analytic methods in probability theory. Uspekhi Mat. Nauk, 5, 5-41.
  * (20) Van Kampen, N.G. (1992) Stochastic processes in physics and chemistry (Vol. 1). Elsevier.
  * (21) Singh, M.S., Singh, M.K. & Singh, R.B. (2024). Stochastic approach to study the properties of the complex patterns observed in cytokine and T cell interaction process. Nonlinear Dynamics, 112(3), 2237-2252.
  * (22) Fano, U. (1947) Ionization yield of radiations. II. The fluctuations of the number of ions. Phys. Rev., 72(1), 26.
  * (23) Rajdl, K., Lansky, P. & Kostal, L. (2020) Fano factor: a potentially useful information. Front. comput. neurosci., 14, 569049.
  * (24) Chanu, A.L., Bhadana, J., & Singh, R.B. (2020) Stochastic fluctuations as a driving force to dissipative non-equilibrium states. J. Phys. A Math. Theor., 53(42), 425002.
  * (25) Bratsun, D., Volfson, D., Tsimring, L.S. & Hasty, J. (2005) Delay-induced stochastic oscillations in gene regulation. Proc. Natl. Acad. Sci. 102 14593-14598.
  * (26) Barrio, M., Burrage, K., Leier, A. & Tian, T. (2006) Oscillatory Regulation of Hes1: Discrete Stochastic Delay Modelling and Simulation. PLoS Comput. Biol. 2(9): e117. DOI: 10.1371/journal.pcbi.0020117.
  * (27) Cai, X. (2007) Exact stochastic simulation of coupled chemical reactions with delays. J. Chem. Phys. 126 124108-1.
  * (28) Anderson, D.F. (2007) A modified next reaction method for simulating chemical systems with time dependent propensities and delays. J. Chem. Phys. 127.21.