# Role of Brownian motion and Néel relaxations
in Mössbauer spectra of magnetic liquids

A.Ya.Dzyublik<EMAIL_ADDRESS>V.Yu.Spivak  Institute for Nuclear Research,
National Academy of Sciences of Ukraine, avenue Nauki 47, 03680 Kyiv, Ukraine

###### Abstract

The absorption cross section of Mössbauer radiation in magnetic liquids is
calculated, taking into consideration both translational and rotational
Brownian motion of magnetic nanoparticles. Stochastic reversals of their
magnetization are also regarded in the absence of external magnetic field. The
role of Brownian motion in ferrofluids is considered in the framework of the
diffusion theory, while for the magnetorheological fluids with large
nanoparticles it is done in the framework of the Langevin’s approach. For
rotation we derived the equation analogous to Langevin’s equation and gave the
corresponding correlation function. In both cases the equations for rotation
are solved in the approximation of small rotations during lifetime of the
excited state of Mössbauer nuclei. The influence of magnetization relaxations
is studied with the aid of the Blume-Tjon model.

###### pacs:

06.20.-f, 42.55.-f, 32.90.+a, 03.5.Nk

##  I Introduction

Suspensions of magnetic nanoparticles (MNPs) attract great attention due to
their numerous applications in technique, medicine and biology [1-17]. It is
provided by large magnetic moment of MNPs, which allows to manipulate them by
moderate magnetic fields. Depending on the dimensions of MNPs, they can be
divided into magnetorheological fluids formed by MNPs with the diameter of the
order of $1\mu$m and ferrofluids with dimensions of MNPs $\sim 10$ nm (see,
e.g., Brazil ). Viscosity of the magnetorheological liquids, being subjected
to the magnetic field, enormously increases, so that they may even transform
into a solid body. This property gives possibility to use such suspensions in
dampers, brakes and clutches Wang .

Ferrofluids are widely used in computers, loudspeakers, semiconductors, motion
controllers, sensors, ink-jet printers, seals, bearings, stepper motors, etc
[1-8].

In medicine ferrofluids are employed in hyperthermia Mody , for drug
delivering to local ill regions of the body and as contrast agents in magnetic
resonance imaging (MRI) med . Recently, much more progressive method of
magnetic particle imaging (MPI) was developed for visualizing MNPs in humans
and animals [11-15]. The advantage of MPI is that it is more fast,
quantitative, and sensitive than MRI.

Note that MNPs are always coated with a polymer shell to prevent their
agglomeration in a solution. The commercial ferrofluids are predominantly
based on magnetite Fe304 particles.

Usually MNPs have single easy-magnetization axis $\zeta$ and their
magnetization ${\bf M}$ tends to be oriented along it or in the opposite
direction, keeping the constant value $|{\bf M}|$. The anisotropy potential
energy of such particles in the absence of external magnetic fields, versus
the angle $\Theta$ between ${\bf M}$ and axis $\zeta$, is represented by two
potential wells at $\Theta=0$ and $\pi$, separated by the potential barrier:

$E=K_{\textrm{\scriptsize{eff}}}V\sin^{2}\Theta,$ (1)

where $K_{\textrm{\scriptsize{eff}}}$ denotes the effective magnetic
anisotropy, $V$ the particle volume. The magnetization, oscillating in one of
the potential wells, from time to time gets sufficient energy to jump over the
barrier into the neighboring well. In the symmetric potential (1) the
magnetization reversals occur with equal rate $w=1/\tau$ in both sides, where
the relaxation time $\tau_{N}$ is determined by Néel’s formula Neel

$\tau_{N}=\tau_{0}\exp{\left(\frac{K_{\textrm{\scriptsize{eff}}}V}{kT}\right)}$
(2)

with the constant factor $\tau_{0}\sim 10^{-9}-10^{-13}s^{-1}$ (see, e.g.,
Landers ).

Effectiveness of the MNPs in different applications strongly depends on such
parameters as the Nèel relaxation time and Brownian rotational relaxation time
of MNPs as well as their relation to temperature and viscosity of the carrier
fluid. The Mössbauer spectroscopy is the most powerful method allowing us to
determine all such characteristics. Soon after discovery the Mössbauer effect
has been applied for investigation of the Brownian motion of nanoparticles in
liquids [18-26]. Foundation of these studies has been laid by Singwi and
Sjölander Singwi , who expressed the absorption cross section of Mössbauer
rays by chaotically moving nanoparticle in terms of the Van Hove auto-
correlation function. Having described the translational motion of the
Brownian particle by diffusion equation, they found that the broadening of
Mössbauer line linearly depends on the ratio of temperature $T$ and viscosity
of the liquid $\eta$. This conclusion was supported experimentally for small
nanoparticles [18-24]. But somewhat later for large particles it was observed
considerable curvature of the line, which was explained theoretically in Refs.
K ; B2 , where the Brownian motion was described by means of the stochastic
Langevin’s equation.

For the first time the role of the Brownian rotation was analyzed by Zatovskiǐ
Zat , who found most strict solution for the absorption cross section. In Ref.
Dz1  this task has been solved in small-angle approximation, taking into
account that during the lifetime of the excited Mössbauer nucleus the root
mean-square angle of the Brownian rotation is usually much less than unity.
Such simple calculations were expanded to the case of ellipsoidal Brownian
particles in Ref. Dz2 . Another more cumbersome approach to the problem was
developed in Ref. Af .

Besides, manifestation of the Néel relaxations along with the Brownian motion
in Mössbauer spectra was studied in Ref. Dz3 , whereas Landers et al. [17]
seem to be the first who observed the Mössbauer spectra in ferrofluids.
Interesting experimental results have been also reported in Refs.[33-36].

In this paper we research in much more detail the impact of Brownian rotation
together with translational diffusion on the shape of Mössbauer spectra. In
particular, for the first time the Brownian rotation of large nanoparticles
will be regarded in the framework of Langevin’s formalism on the same footing
as the translational motion.

For ferrofluids we shall apply a simplest relaxation model, when the
magnetization vector of the particle ${\bf M}(t)$ makes stochastic jumps
between the values ${\bf M}$ and $-{\bf M}$ along the easy axis $\zeta$.
Respectively, the magnetic field at the nucleus, being antiparallel to ${\bf
M}(t)$, takes the values ${\bf h}_{0}(t)={\bf h}_{0}f(t)$ with $f(t)=\pm 1$.
The magnetic field ${\bf h}_{0}$ causes splitting of the nuclear sublevels
giving rise to a Zeeman sextet. For generality we adopt that along the field
${\bf h}_{0}$ there is an electric field gradient, which ensures a quadrupole
splitting of the lines. This model was previously applied to calculations of
Mössbauer spectra by Blume and Tjon Blume .

##  II Basic equations

In order to separate the translational and rotational motion we first
introduce the coordinate frame $\\{x,\;y,\;z\\}$ with the origin in the center
of the particle and axis $z$ along the beam of incident $\gamma$-quanta. In
addition, we introduce the frame $\\{\xi,\;\eta,\;\zeta\\}$ with axis $\zeta$
along the easy-magnetization axis of the particle. Position of the Mössbauer
nucleus 57Fe in the laboratory frame $\\{x_{L},\;y_{L},\;z_{L}\\}$ is
determined by the radius-vector

${\bf X}={\bf R}+{\bf r}+{\bf u},$ (3)

where the vector ${\bf R}$ indicates position of the center of the Brownian
particle, ${\bf r}$ specifies the equilibrium position of the nucleus in the
frame $\\{x,\;y,\;z\\}$, and ${\bf u}$ is the displacement from this cite.

Random reversals of the magnetization ${\bf M}$ and the Brownian motion are
independent processes. Therefore the absorption cross section of
$\gamma$-quanta with the energy $E=\hbar\omega$ and wave vector
${\bm{\kappa}}$ by the Mössbauer nucleus 57Fe, embedded in the Brownian
particle, may be written as Dz3

$\displaystyle\sigma_{a}(\omega)=\frac{\sigma_{0}\Gamma_{a}}{2}e^{-2W}$ (4)
$\displaystyle\times\mbox{Re}\int_{0}^{\infty}\frac{dt}{\hbar}e^{i(\omega-\omega_{a})t-\Gamma_{a}t/2\hbar}G_{B}({\bm{\kappa}},t)G_{N}(t),$

where $\sigma_{0}$ is the resonant value of the absorption cross section of
$\gamma$-quanta by a fixed nucleus in the absence of the hyperfine structure,
$E_{a}=\hbar\omega_{a}$ and $\Gamma_{a}$ are the energy and width of the
resonant level of the absorbing nucleus, $e^{-2W}$ is the Debye-Waller factor,
$G_{B}({\bm{\kappa}},t)$ denotes the Fourier-transform of the classical
autocorrelation function for the Brownian motion, $G_{N}(t)$ the correlation
function for the Néel relaxations of magnetization.

This cross section is to be averaged over the energy distribution of
$\gamma$-quanta emitted by a source without recoil

$w_{e}(E)=\frac{2\pi}{\Gamma_{e}}\frac{1}{(E-E_{e}-s)^{2}+(\Gamma_{e}/2)^{2}},$
(5)

where $s=vE_{e}/c$ denotes the Doppler shift for a source, moving with the
velocity $v$ relative to an absorber. Then experimentally measured cross
section takes the form

$\displaystyle\sigma_{a}(s)=\frac{\sigma_{0}\Gamma_{a}}{2}e^{-2W}$ (6)
$\displaystyle\times\mbox{Re}\int_{0}^{\infty}\frac{dt}{\hbar}e^{ist/\hbar-\Gamma_{0}t/2\hbar}G_{B}({\bm{\kappa}},t)G_{N}(t),$

where $\Gamma_{0}=\Gamma_{e}+\Gamma_{a}$ means the width observed when any
broadening due to Brownian motion or Néel relaxations is absent.

For spherical particles the translational and rotational Brownian motions are
separated, so that

$G_{B}({\bm{\kappa}},t)=G_{t}({\bm{\kappa}},t)G_{r}({\bm{\kappa}},t),$ (7)

where $G_{t}({\bm{\kappa}},t)$ and $G_{r}({\bm{\kappa}},t)$ are the Fourier
transforms of the correlation functions for translational motion and rotation,
respectively.

##  III Correlation functions

In this section we shall give the correlation functions for the translational
and rotational Brownian motion of spherical nanoparticles in a liquid,
provided by corresponding diffusion equations. Besides, the correlator
responsible for the Néel relaxations of the MNPs magnetization, derived in
Ref. Blume , will be reproduced below in somewhat changed form.

### 3.1 Translational Brownian motion

For the translational Brownian motion, described by simple diffusion equation,
the correlation function has the form Singwi

$G^{t}_{s}({\bf R},t)=(4\pi D_{t}|t|)^{-3/2}\exp[-R^{2}/4D_{t}|t|].$ (8)

Here we suppose that at the initial moment $t=0$ the particle is located in
the origin of the laboratory frame. The Fourier-transform of the function (8)
reads

$G_{s}^{t}({\bm{\kappa}},t)=\exp(-\kappa^{2}D_{t}|t|).$ (9)

As a result, the broadening of the Mössbauer line caused by translational
diffusion of a spherical nanoparticle is given by Singwi

$\Delta\Gamma_{t}=2\hbar\kappa^{2}D_{t},$ (10)

where the translational diffusion coefficient

$D_{t}=\frac{kT}{6\pi\eta a_{h}},$ (11)

$\eta$ is the viscosity coefficient of the liquid, $a_{h}$ the hydrodynamic
radius of the nanoparticle being a sum of the core radius $a_{c}$ and a
thickness of its coating $d$.

### 3.2 Brownian rotation

The mean-square angle of rotation $\Delta_{r}^{2}$ of the Brownian particle in
a liquid during time $t$ is determined by Einstein’s formula Ein

$\Delta^{2}_{r}=2D_{r}t,$ (12)

depending on the rotation diffusion coefficient

$D_{r}=\frac{kT}{8\pi\eta a_{h}^{3}},$ (13)

where $a_{h}$ is the hydrodynamic radius of the particle.

Let us estimate the $\Delta^{2}_{r}$ for rotation during the time $t$ of the
order of the lifetime $\tau_{N}=\hbar/\Gamma_{a}=141$ ns for 57Fe. We take the
parameters of the experiment Landers , which correspond to maximal value of
$\Delta^{2}_{r}$: $a_{h}=7$ nm and $\eta=22.5$ cp (viscosity of the 70%
glycerol solution at $T=293$ K). In this case $\Delta^{2}_{r}\approx 4\cdot
10^{-3}$. In all other measurements Landers , corresponding to lower
temperatures and larger particles, $\Delta^{2}_{r}$ is much less. Thus, we can
really treat the Brownian rotation in the small-angle approximation.

The Fourier-transform of the rotational correlation function is calculated
with the aid of the probability density $W({\bf n},{\bf n}_{0};t)$ of the
Brownian rotation from ${\bf n}_{0}$ to ${\bf n}$ during time $t$:

$\displaystyle G_{s}^{r}({\bm{\kappa}},t)=\int d{\bf n}W({\bf n},{\bf
n}_{0};t)e^{i{\bm{\kappa}}({\bf n}-{\bf n}_{0})r},$ (14)

where orientation of the unit vectors ${\bf n}_{0}={\bf r}(0)/r$ and ${\bf
n}={\bf r}(t)/r$ in the frame $\\{x,\;y,\;z\\}$ are determined by the
spherical angles $\theta_{0},\varphi_{0}$ and $\theta,\varphi$, respectively.
The function $W({\bf n},{\bf n}_{0};t)$ is looked for as a solution of the
rotational diffusion equation Leon

$\frac{\partial W}{\partial
t}=D_{r}\left[\frac{\partial^{2}W}{\partial\theta^{2}}+\cot\theta\frac{\partial
W}{\partial\theta}+\frac{1}{\sin^{2}\theta}\frac{\partial^{2}W}{\partial\varphi^{2}}\right]$
(15)

with the initial condition

$W({\bf n},{\bf n}_{0};0)=\delta({\bf n}-{\bf n}_{0}).$ (16)

The probability of all possible events equals unity, therefore the probability
density is normalized as

$\int_{0}^{\pi}\sin\theta d\theta\int_{0}^{2\pi}d\varphi W({\bf n},{\bf
n}_{0};t)=1.$ (17)

Let us introduce one more frame $x^{\prime},\;y^{\prime},\;z^{\prime}$, whose
axis $z^{\prime}$ is directed along ${\bf n}_{0}$. The orientation of the
vector ${\bf n}$ in this frame is determined by the spherical angles
$\vartheta,\phi$. In the approximation of small Brownian rotations, when

$({\bf n}-{\bf n}_{0})^{2}\approx\vartheta^{2}<<1,$ (18)

the basic equation (15) transforms to

$\frac{\partial W}{\partial
t}=D_{r}\left[\frac{\partial^{2}W}{\partial\vartheta^{2}}+\frac{1}{\vartheta}\frac{\partial
W}{\partial\vartheta}+\frac{1}{\vartheta^{2}}\frac{\partial^{2}W}{\partial\phi^{2}}\right].$
(19)

Notice that the same form has the equation, which describes diffusion of the
point-like particle on the plane, expressed in polar coordinates with the
radial coordinate $\vartheta$ and azimute angle $\phi$. Hence, a solution of
Eq. (19), satisfying the initial condition (16), is

$W({\bf n},{\bf n}_{0};t)=\frac{1}{4\pi
D_{r}|t|}\exp\left(-\frac{\vartheta^{2}}{4D_{r}|t|}\right).$ (20)

From here we see that the mean-square rotation angle is really determined by
Eq. (12).

In order to calculate the Fourier-transform of $W({\bf n},{\bf n}_{0};t)$ we
shall express the components $n_{x},\;n_{y},\;n_{z}$ of the of unit vector
${\bf n}$ in spherical angles:

$n_{x}=\sin\theta\cos\varphi,\quad n_{y}=\sin\theta\sin\varphi,\quad
n_{z}=\cos\theta.$ (21)

Simple calculation gives

$({\bf n}-{\bf
n}_{0})^{2}=(\theta-\theta_{0})^{2}+\sin^{2}\theta_{0}(\varphi-\varphi_{0})^{2}.$
(22)

Comparing (18) with (22) we rewrite (20) as

$\displaystyle W({\bf n},{\bf n}_{0};t)=\frac{1}{4\pi D_{r}|t|}\qquad\qquad$
(23)
$\displaystyle\times\exp\left(-\frac{(\theta-\theta_{0})^{2}}{4D_{r}|t|}\right)\exp\left(-\frac{\sin^{2}\theta_{0}(\varphi-\varphi_{0})^{2}}{4D_{r}|t|}\right).$

Then starting from the formula

$\displaystyle G_{s}^{r}({\bm{\kappa}},t)=\int_{-\infty}^{\infty}\sin\theta
d\theta\exp[-i\kappa r\sin\theta_{0}(\theta-\theta_{0})]$
$\displaystyle\times\int_{-\infty}^{\infty}d\varphi W({\bf n},{\bf n}_{0};t),$
(24)

we arrive at the Fourier-transform of the rotation correlation function:

$G_{s}^{r}({\bm{\kappa}},t)=\exp\left[-\kappa^{2}D_{r}r^{2}\sin^{2}\theta_{0}|t|\right].$
(25)

### 3.3 Magnetization relaxations

Following Ref. Blume  we suppose that there is an electric field gradient
along the magnetic field ${\bf h}_{0}$ at the nucleus 57Fe. The constant field
${\bf h}_{0}$ gives rise to Zeeman splitting of sublevels $1/2,M_{g}$ and
$3/2,M_{e}$ of the ground and excited nuclear states, respectively. Here
$M_{g}$ and $M_{e}$ are the projections of the nuclear spin on the direction
${\bf h}_{0}$. In the fluctuating field ${\bf h}_{0}(t)$ the Mössbauer
spectrum is described by the correlator $G_{N}(t)$ Blume

$\displaystyle
G_{N}(t)=\sum_{j=1}^{6}J_{j}(\beta)\exp[-iQ(3M_{e}^{2}-15/4)t/\hbar]$ (26)
$\displaystyle\times\left(\exp\left[i\alpha_{j}\int_{0}^{t}f(t^{\prime})dt^{\prime}\right]\right)_{av},$

where $(...)_{av}$ implies the stochastic averaging, $f(t)=\pm 1$, the
parameter $Q$ determines a quadrupole shift of the lines, the factors
$J_{j=M_{g}\to M_{e}}(\beta)$ are relative intensities of the Zeeman sextet:

$\displaystyle J_{1=-1/2\to-3/2}(\beta)=J_{6=1/2\to
3/2}(\beta)=\frac{3}{16}(1+\cos^{2}\beta),$ $\displaystyle
J_{2=-1/2\to-1/2}(\beta)=J_{5=1/2\to
1/2}(\beta)=\frac{1}{4}\sin^{2}\beta,\qquad$ (27) $\displaystyle J_{3=-1/2\to
1/2}(\beta)=J_{4=1/2\to-1/2}(\beta)=\frac{1}{16}(1+\cos^{2}\beta),$

depending on the angle $\beta$ between the wave vector of $\gamma$-quanta
${\bm{\kappa}}$ and magnetization $+{\bf M}$. Here the lines of the Zeeman
sextet are numerated in the order of growing energy. In the absence of
external magnetic fields, when the particles are oriented randomly, the
averaged relative intensities of the lines $J_{i}=\langle J_{i}(\beta)\rangle$
are

$J_{1}=J_{6}=\frac{3}{12},\;\;J_{2}=J_{5}=\frac{2}{12},\;\;J_{3}=J_{4}=\frac{1}{12}.$
(28)

Then the stochastic averaging results in Blume

$\displaystyle\left(\exp\left[i\alpha_{j}\int_{0}^{t}f(t^{\prime})dt^{\prime}\right]\right)_{av}=$
(29) $\displaystyle=(\cos x_{j}wt+x_{j}^{-1}\sin x_{j}wt)\exp(-wt),$

with parameters

$\displaystyle x_{j}=[(\alpha_{j}/w)^{2}-1]^{1/2},$
$\displaystyle\alpha_{j}=\alpha_{eg}=(g_{g}M_{g}-g_{e}M_{e})\mu_{N}h_{0}/\hbar,$
(30)

depending on the nuclear magneton $\mu_{N}$ and gyromagnetic ratios $g_{g}$,
$g_{e}$ of the ground and excited states, respectively. From now on, for
brevity, we shall omit the exponential in Eq. (29), associated with the
quadrupole splitting. Once $Q\neq 0$, in all equations below the Doppler shift
$s$ is to be replaced by $s-Q(3M_{e}^{2}-15/4)$.

##  IV Absorption cross section

Substituting (9), (25) and (29) into (6) one finds the absorption cross
section

$\sigma_{a}(s)=\frac{\sigma_{0}\Gamma_{a}}{4}e^{-2W}\mbox{Re}\sum_{j=1}^{6}J_{j}\left[\left(1-\frac{i}{x_{j}}\right)\frac{i}{s+x_{j}\hbar
w+i\Gamma_{\textrm{\scriptsize{eff}}}/2}+\left(1+\frac{i}{x_{j}}\right)\frac{i}{s-x_{j}\hbar
w+i\Gamma_{\textrm{\scriptsize{eff}}}/2}\right],$ (31)

where the effective width

$\Gamma_{\textrm{\scriptsize{eff}}}=\Gamma_{0}+\Delta\Gamma_{B}+\Delta\Gamma_{N},$
(32)

$\Delta\Gamma_{N}=2\hbar w$ is a broadening due to the Néel relaxations, while
the Brownian broadening is given by a sum

$\Delta\Gamma_{B}=\Delta\Gamma_{t}+\Delta\Gamma_{r}({\bf r})$ (33)

with the rotational contribution

$\Delta\Gamma_{r}({\bf r})=2\hbar\kappa^{2}D_{r}r^{2}\sin^{2}\theta_{0},$ (34)

depending on the coordinates of the Mössbauer nucleus 57Fe.

For uniform distribution of these nuclei in nanoparticle the averaged cross
section is defined by

$\langle\sigma_{a}(s)\rangle=\frac{3}{2a_{c}^{3}}\int_{0}^{a_{c}}r^{2}dr\int_{0}^{\pi}\sigma_{a}(s)\sin\theta_{0}d\theta_{0}.$
(35)

Having substituted here the expression (31) we introduce new variables
$\xi=\cos\theta_{0}$ and $\rho=r/a_{c}$ to obtain

$\displaystyle\langle\sigma_{a}(s)\rangle=\frac{3\sigma_{0}\Gamma_{a}}{4}e^{-2W}\mbox{Re}\sum_{j=1}^{6}J_{j}\int_{0}^{1}\rho^{2}d\rho$
(36)
$\displaystyle\times\left[\left(1-\frac{i}{x_{j}}\right)I_{j}^{+}(\rho)+\left(1+\frac{i}{x_{j}}\right)I_{j}^{-}(\rho)\right],$

where $I_{j}^{\pm}(\rho)$ stands for the integral

$I_{j}^{\pm}(\rho)=\int_{0}^{1}\frac{d\xi}{A_{j}^{\pm}(\rho)-B(\rho)\xi^{2}}$
(37)

with

$\displaystyle
B=i\hbar\kappa^{2}D_{r}a_{c}^{2}\rho^{2}=(3i/8)\Delta\Gamma_{t}(a_{c}/a_{h})^{2}\rho^{2}$
(38)

and

$A_{j}^{\pm}=s\pm x_{j}w+i(\Gamma_{0}+\Delta\Gamma_{t}+\Delta\Gamma_{N})/2+B.$
(39)

Trivial integration gives

$I_{j}^{\pm}(\rho)=\frac{1}{\sqrt{A_{j}^{\pm}B}}\operatorname{Arth}\left({\sqrt{B/A_{j}^{\pm}}}\right).$
(40)

In the case of slow relaxations, when $w<<|\alpha_{j}|$ and respectively
$x_{j}>>1$ as well as $x_{j}w\approx\alpha_{j}$, the cross section reduces to
Zeeman’s pattern with broadened lines:

$\displaystyle\langle\sigma_{a}(s)\rangle=\frac{\sigma_{0}\Gamma_{a}}{2}e^{-2W_{a}}\mbox{Re}\sum_{j=1}^{6}J_{j}$
(41)
$\displaystyle\times\int_{0}^{1}\rho^{2}d\rho\frac{1}{\sqrt{A_{j}B}}\operatorname{Arth}\left(\sqrt{B/A_{j}}\right),$

where $B$ is again determined by formula (38), while $A_{j}$ takes the form

$A_{j}=s-\hbar\alpha_{j}+i(\Gamma_{0}+\Delta\Gamma_{t}+\Delta\Gamma_{N})/2+B.$
(42)

In the opposite limit of very rapid relaxations as $w>>|\alpha_{j}|$ the
nucleus only feels an average zero magnetic field. In this case the cross
section collapses to single line 111if $Q\neq 0$ the spectrum collapses to a
quadrupole doublet.

$\displaystyle\langle\sigma_{a}(s)\rangle=\frac{3\sigma_{0}\Gamma_{a}}{2}e^{-2W_{a}}\mbox{Re}\int_{0}^{1}\rho^{2}d\rho$
(43)
$\displaystyle\times\frac{1}{\sqrt{AB}}\operatorname{Arth}\left(\sqrt{B/A}\right),$

where $B$ remains the same, while $A$ becomes

$A=s+i(\Gamma_{0}+\Delta\Gamma_{t})/2+B.$ (44)

Note also that the same expression (43) describes the Mössbauer spectra of
nonmagnetic Brownian particles.

The formulas considerably simplify, if we average only $\Delta\Gamma_{r}({\bf
r})$ instead of the whole cross section (31). Then

$\langle r^{2}\sin^{2}\theta_{0}\rangle=0.4a_{c}^{2}$ (45)

and

$\langle\Delta\Gamma_{r}({\bf
r})\rangle=0.6\hbar\kappa^{2}D_{t}(a_{c}/a_{h})^{2}.$ (46)

In this case the averaged cross section is determined by the same formula (31)
but with the Brownian broadening $\Delta\Gamma_{B}$ replaced by

$\langle\Delta\Gamma_{B}\rangle=2\hbar\kappa^{2}D_{t}[1+0.3(a_{c}/a_{h})^{2}].$
(47)

As to the cross section (43, it reduces to

$\langle\sigma_{a}(s)\rangle=\frac{\sigma_{0}\Gamma_{a}}{4}e^{-2W_{a}}\frac{\Gamma_{0}+\langle\Delta\Gamma_{B}\rangle}{s^{2}+(\Gamma_{0}+\langle\Delta\Gamma_{B}\rangle)^{2}/4}.$
(48)

If a contribution of the rotation is neglected, Eq. (48) coincides with the
result of Singwi and Sjolander Singwi .

The role of rotational diffusion is illustrated by Fig.1, where all the cross
sections are calculated in units $(3\sigma_{0}\Gamma_{a}/2)e^{-2W}$ as a
function of the dimensionless parameter $2s/(\Gamma_{0}+\Delta\Gamma_{t})$ for
the case, when $\Delta\Gamma_{t}=4\Gamma_{0}$ and $a_{c}=a_{h}$. The cross
section (43) is drawn by the solid line. The Singwi-Sjolander’s curve, given
by Eq. (48) with rotational contribution $\langle\Delta\Gamma_{r}({\bf
r})\rangle=0$, by the dashed one. In addition, the approximate curve (48) with
$\langle\Delta\Gamma_{B}\rangle=2.6\hbar\kappa^{2}D_{t}$, is represented by
the dash-dotted line. It is seen that it surprisingly well approximates the
exact result (43).

Figure 1: Dependance of the absorption cross sections on the Doppler shift
$s$, expressed in units $(\Gamma_{0}+\Delta\Gamma_{t})/2$. The cross section
calculated by Eq. (43) is drawn by solid line, the result of Singwi and
Sjolander Singwi  by dashed one, the approximate expression (48) by dash-
dotted.

##  V Approach based on Langevin’s equation

The correlation function (8) is not valid at small times. More correctly self-
diffusion is described by the correlation function Singwi ; Chan

$G_{t}({\bf
R},t)=[2\pi\rho_{t}(t)]^{-3/2}\exp\left[-\frac{R^{2}}{2\rho_{t}(t)}\right],$
(49)

where

$\rho_{t}(t)=\frac{2D_{t}}{\beta_{t}}[\beta_{t}t-1+\exp(-\beta_{t}t)],$ (50)

with

$\beta_{t}=6\pi a\eta/m=\frac{kT}{D_{t}m}.$ (51)

The parameter $\beta_{t}^{-1}$ means the characteristic (relaxation) time for
the Brownian translational motion.

At $t>>\beta_{t}^{-1}$ the correlation functions (8) and (49) coincide. In the
opposite limit of $t<<\beta_{t}^{-1}$ an employment of diffusion approach
leads to paradox, remarked in Ref. Chuev1 . Really, from the diffusion
equation it follows that at $t\geq 0$ the mean-square displacement $\langle
x^{2}\rangle=2D_{t}t$, and therefore the root mean-square velocity along the
axis $x$, given by ${\sqrt{\langle v^{2}_{x}\rangle}}=\sqrt{2D_{t}/t}$, tends
to infinity if $t\to 0$. At the same time, in correspondence with (50), at
$t<<\beta_{t}^{-1}$ the mean-square displacement $\langle
x^{2}\rangle=D_{t}\beta_{t}t^{2}$, hence $\sqrt{\langle
v^{2}_{x}\rangle}=\sqrt{D_{t}\beta_{t}}=(kT/m)^{1/2}$. Thus, the mean kinetic
energy of the Brownian particle $\bar{E}_{\textrm{\scriptsize{kin}}}$ at $t\to
0$ occurs to be determined by the same expression as
$\bar{E}_{\textrm{\scriptsize{kin}}}$ for the molecules of the ideal gas:

$\bar{E}_{\textrm{\scriptsize{kin}}}=3kT/2.$ (52)

The function (49) was derived by Chandrasekhar Chan  from Langevin’s equation

$\frac{d^{2}{\bf R}}{dt^{2}}=-\beta_{t}\frac{d{\bf R}}{dt}+{\bf
F}_{\textrm{\scriptsize{rand}}}(t)/m.$ (53)

Here on the right-hand side the first term is responsible for the dynamical
friction, ${\bf F}_{\textrm{\scriptsize{rand}}}(t)/m$ are random forces acting
on the Brownian particle.

Let us find now analogous correlation function for the rotational Brownian
motion. For this aim we shall first derive the equation similar to Eq. (53),
starting from well-known relationship for the angular momentum ${\bf L}$ of
the rotating rigid body and the total torque ${\bf K}$ acting on it Landau :

$\frac{d{\bf L}}{dt}={\bf K}.$ (54)

We take into account that the angular momentum for a rigid sphere of the
radius $a$ is bound to its angular frequency of rotation ${\bm{\omega}}$ by

${\bf L}=\mathfrak{J}{\bm{\omega}},$ (55)

where the inertia moment of the sphere

$\mathfrak{J}=0.4ma_{h}^{2}.$ (56)

The ${\bf K}$ equals a sum of the friction torque Kirch

${\bf K}_{\textrm{\scriptsize{fr}}}=-8\pi\eta a_{h}^{3}{\bm{\omega}}.$ (57)

and torques ${\bf K}_{\textrm{\scriptsize{rand}}}(t)$ due to random forces.

Inserting (55), (57) into (54) one gets the equation, governing the stochastic
rotational motion:

$\frac{d{\bm{\omega}}}{dt}=-\beta_{r}{\bm{\omega}}+{\bf
K}_{\textrm{\scriptsize{rand}}}(t)/\mathfrak{J},$ (58)

where

$\beta_{r}=\frac{20\pi\eta a_{h}}{m}=\frac{5}{2}\frac{kT}{D_{r}ma_{h}^{2}}.$
(59)

Keeping in mind that for small rotations ${\bm{\omega}}$ is perpendicular to
${\bf n}_{0}$ and equals $\omega=d\vartheta/dt$, we transform (58) to the
equation

$\frac{d^{2}\vartheta}{dt^{2}}=-\beta_{r}\frac{d\vartheta}{dt}+K_{\textrm{\scriptsize{rand}}}(t)/\mathfrak{J},$
(60)

formally equivalent to Langevin’s equation (53) in one-dimensional case. Here
in the same approximation we ignore the boundary conditions for the angle
$\vartheta$ and accept that it ranges from $-\infty$ to $\infty$. Further
repeating derivation, done by Chandrasekhar Chan , one gets the correlation
function

$G_{r}(\vartheta,t)=[2\pi\rho_{r}(t)]^{-1/2}\exp[-\vartheta^{2}/2\rho_{r}(t)],$
(61)

where the function $\rho_{r}(t)$ is again defined by Eq. (50), but with index
$t$ replaced by $r$. As to the Fourier-transform, it is given now by

$G_{r}({\bm{\kappa}},t)=\exp\left[-\kappa^{2}\rho_{r}(t)r^{2}\sin^{2}\theta_{0}/2\right].$
(62)

Combining these equations we get in the slow-relaxation limit the cross
section as a superposition of six lines:

$\sigma_{a}(s)=\sum_{j=1}^{6}\sigma_{a}^{(j)}(s),$ (63)

each of them is given by

$\displaystyle\sigma_{a}^{(j)}(s)=\frac{\sigma_{0}\Gamma_{a}}{4}e^{-2W}e^{b_{t}+b_{r}}J_{j}$
(64)
$\displaystyle\times\sum_{n,k=0}^{\infty}\frac{(-b_{t})^{n}}{n!}\frac{(-b_{r})^{k}}{k!}\frac{\Gamma_{n,k}}{(s-\hbar\alpha_{j})^{2}+\left(\frac{\Gamma_{n,k}}{2}\right)^{2}},$

where we introduced the parameters

$b_{t}=\frac{\kappa^{2}D_{t}}{\beta_{t}},\qquad
b_{r}=\frac{\kappa^{2}D_{r}r^{2}\sin^{2}\theta_{0}}{\beta_{r}}$ (65)

and the widths

$\Gamma_{n,k}=\Gamma_{\textrm{\scriptsize{eff}}}+2n\hbar\beta_{t}+2k\hbar\beta_{r}.$
(66)

For estimations it is sufficient to replace averaging of the cross section
over ${\bf r}$ by the averaging of $r^{2}\sin^{2}\theta_{0}$ i.e., to replace
this product by $\langle r^{2}\sin^{2}\theta_{0}\rangle=0.4a_{c}^{2}$. Then
one has the relation

$\langle b_{r}({\bf r})\rangle=0.045b_{t},$ (67)

which allows to set in (64) $b_{r}=0$.

Figure 2:  Brownian broadening of the absorption lines vs temperature for the
magnetite nanoparticles with diameter $a_{h}=a_{c}=700$ nm in the 60%
glycerol-water mixture, calculated in the diffusion approach (1) and
Langevin’s one (2)

In order to find now the integral width of $j$th Zeeman line we employ
standard formula

$\Gamma_{\textrm{\scriptsize{int}}}=\frac{2}{\pi\sigma^{(j)}_{a}(\hbar\alpha_{j})}\int_{-\infty}^{\infty}\sigma^{(j)}_{a}(s)ds.$
(68)

Simple calculation gives

$\displaystyle\Gamma_{\textrm{\scriptsize{int}}}=e^{-b_{t}}\left[\sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\frac{(b_{t})^{n}}{\langle\Gamma_{n}\rangle}\right]^{-1},$
(69)

where the width
$\langle\Gamma_{n}\rangle=\langle\Gamma_{\textrm{\scriptsize{eff}}}\rangle+2n\hbar\beta_{t}$
and the averaged effective width is determined by

$\langle\Gamma_{\textrm{\scriptsize{eff}}}\rangle=\Gamma_{0}+2\hbar\kappa^{2}D_{t}[1+0.3(a_{c}/a_{h})^{2}]+\Delta\Gamma_{N}.$
(70)

Significance of Langevin’s approach is illustrated by Fig.2, where the
Brownian broadening $\langle\Delta\Gamma_{B}\rangle$ for magnetite MNPs with
diameter 700 nm is calculated by Eq. (47) (curve 1) and by Eq. (69) (curve 2).
In the last case we put $b_{r}=0$, that is describe the translational Brownian
motion by Langevin’s equation and the rotational one by simple diffusion
equation. Note that in previous articles K ; B2  the Brownian rotation has not
been taken into consideration at all.

##  VI Discussion

It can be easily shown that the considered relaxations do not influence on the
overall probability of the absorption without recoil. Really, since

$\int_{-\infty}^{\infty}\exp(ist/\hbar)ds=2\pi\hbar\delta(t)$ (71)

and $G_{B}({\bm{\kappa}},0)G_{N}(0)=1$ one gets

$\int_{-\infty}^{\infty}\sigma_{a}(s)ds=\frac{\pi}{2}\sigma_{0}\Gamma_{a}e^{-2W}.$
(72)

Thus, the square of the Mössbauer spectrum keeps constant value.

Figure 3:  Brownian broadening of the absorption lines vs temperature for the
magnetite nanoparticles with diameters 40 nm and 140 nm dissolved in 96%
glycerol-water mixture. The experiment Chuev1  is presented by circles, our
calculations by solid curves.

We have derived the expression (31) for the absorption cross section
$\sigma_{a}(s)$ of Mössbauer radiation by MNPs suspended in a liquid. It will
coincide with the result of Blume and Tion Blume , if in the effective width
$\Gamma_{\textrm{\scriptsize{eff}}}$ we omit a contribution $\Delta\Gamma_{B}$
of the Brownian motion. But the cross section, averaged over uniform
distribution of the nuclei 57Fe over the particle, takes a combersome form.
The situation significantly simplifies if we replace averaging of
$\sigma_{a}(s)$ by the averaging only of $\Delta\Gamma_{B}$. As it is seen in
Fig.1, this procedure provides very good result.

We compared our calculations with the experimental data Cher  of Cherepanov et
al. (see Fig.3). From the Mössbauer spectra of ferrofluids they subtracted the
spectra of dried samples. It enabled them to extract the contribution only of
Brownian motion into the broadening of the spectral lines. The calculated
dependence of the Brownian broadening $\Delta\Gamma_{B}$ on the temperature is
presented in Fig.3 by solid lines and the experimental data by circles. The
calculations very well agree with the experiment for large particles having
diameter 140 nm, and at the same time they terribly deviate for small ones
with the dimensions 40 nm. This deviation contradicts to the fact that the
viscosity of liquid exponentially decreases with growing temperature, and as a
consequence this should ensure for small MNPs the same behavior of
$\Delta\Gamma_{B}$ as for large MNPs.

The next part of the paper will be addressed to Mössbauer spectra in
ferrofluids subjected to external magnetic fields.

## References

  * (1) R. E. Rosensweig, Ferrohydrodynamics, Cambridge Univ. Press, Cambridge, London (1985); republished by Dover.Publ.Inc., New York (1997).
  * (2) Magnetic Fluids and Applications Handbook , edit. by B. Berkovski and V. Bashtovoy, (Begell House, Wallingford (1996)).
  * (3) C. Scherer and A. M. Figueiredo Neto, Brazil. J. Phys. 35, 718 (2005).
  * (4) D. H. Wang and W. H. Liao, Smart Mater. Struct. 20, 023001 (2011).
  * (5) K. Ray and R. Moskowitz, J. Magn. Magn. Mater. 85, 233 (1990).
  * (6) K. Ray, R. Moskowitz, and R. Casciari, J. Magn. Magn. Mater. 149, 174 (1995).
  * (7) R. Pérez-Castillejos, J. A. Plaza, J. Esteve, P. Losantos, M. C. Acero, C. Cane, and F. Serra-Mestres, Sens. Actuators A 84, 176 (2000).
  * (8) E. H. Kim, H. S. Lee, B. K. Kwak, and B. K. Kim, J. Magn. Magn. Mater. 289, 328 (2005).
  * (9) V. V. Mody, A. Singh, and B. Wesley, Eur. J. Nanomed. 5, 11 (2013).
  * (10) K. M. Krishnan, IEEE Trans. Magn. 46, 2523 (2010).
  * (11) B. Gleich and J. Weizenecker, Nature 435, 1214 (2005).
  * (12) R. M. Ferguson, K. R. Minard, A. P. Khandhar and K. M. Krishnan, Med. Phys. 38, 1619 (2011).
  * (13) R. J. Deissler, Yong Wu, and M. A. Martens, Med. Phys. 41, 012301 (2014).
  * (14) Yufen Xiao and Jianzhong Du, J. Mater. Chem. B 8, 354 (2020).
  * (15) H. T. Kim Duong, A. Abdibastami, L. Gloag, L. Barrera, J. J. Gooding, and R. D. Tilley, Nanoscale 14, 13890 (2022).
  * (16) L. Néel, Ann. Geophys. 5, 99 (1949).
  * (17) J. Landers, S. Salamon, H. Remmer, F.Ludwig, H.Wende, Nano Lett. 16, 1150 (2016).
  * (18) P. P. Craig and N. Sutin, Phys. Rev. Let., 11, 460 (1963).
  * (19) D. St. P. Bunbury, J. A. Elliott, H. E. Hall, and J. M. Williams, Phys. Let. 6, 34 (1963).
  * (20) V. I. Lisichenko, Ukr. J. Phys. 9, 1376 (1964).
  * (21) V. N. Dubinin, et al., Ukr. J. Phys 11, 619 (1966).
  * (22) T. Bonchev, P. Aidemirski, I. Mandzhukov, N. Nedyalkova, B. Skorchev, and A. Strigachev, Sov. Phys. JETP 23, 42 (1966).
  * (23) K. P. Singh and J. G. Mullen, Phys. Rev. A 6, 2354 (1972).
  * (24) V. N. Dubinin, Ukr. J. Phys. 13, 1547 (1968).
  * (25) C. L. Kordyuk, V. I. Lisichenko, 0. L. Orlov, N. N. Polovina, and A. N. Smoilovskii, Sov. Phys. JETP 25, 400 (1967).
  * (26) V. G. Bhide, R. Sundaram, H. C. Bhasin, and T. Bonchev, Phys. Rev. 3, 673 ((1971).
  * (27) K. S. Singwi and A. Sjölander, Phys. Rev. 120, 1093 (1960).
  * (28) A. V. Zatovskiǐ, Sov. Phys. JETP 32, 274 (1971).
  * (29) A. Ya. Dzyublik, Ukr. J. Phys. 18, 1454 (1973).
  * (30) A. Ya. Dzyublik, Sov. Phys. JETP 40, 763 (1975).
  * (31) A. M. Afanas’ev, P. V. Hendriksen, and S. Mørup, Hyperf. Inter. 88, 35 (1994).
  * (32) A. Ya. Dzyublik, Ukr. J. Phys. 23, 881 (1978).
  * (33) Joachim Landers, Soma Salamon, Hilke Remmer, Frank Ludwig, and Heiko Wende, ACS Appl. Mater. Interfaces 11, 3160 (2019).
  * (34) M. A. Chuev, V. M. Cherepanov, M. A. Polikarpov, R. R. Gabbasov, and A. Yu. Yurenya, JETP Lett. 108, 59 (2018).
  * (35) R. Gabbasov, A. Yurenya, A. Nikitin, V. Cherepanov, M. Polikarpov, M. Chuev, A. Majouga, and V. Panchenko, J. Magn. Magn. Mater. 475, 146 (2019).
  * (36) V. M. Cherepanov, et al., Crystal. Rep. 65, 398 (2020).
  * (37) M. Blume, J. A. Tjon, Phys. Rev. 165, 446 (1968).
  * (38) S. Chandrasekhar, Rev. Mod. Phys. 15, 1 (1943).
  * (39) A. Einstein, Ann. Phys. 19, 371 (1906).
  * (40) M. A. Leontovich, Introduction to Thermodynamics. Statistical Physics (High School, Moscow, 1983) (in Russian).
  * (41) L.D.Landau, E.M.Lifshitz, Mechanics 3d Edition (Elsevier Ltd., Oxford, 1976).
  * (42) G. Kirchhoff, Vorlesungen über Mechanik (B.G.Teubner, Leipzig, 1897).


# Novel Symmetry-preserving Neural Network Model for Phylogenetic Inference

Xudong Tang
Wisconsin Institute for Discovery
University of Wisconsin-Madison
Madison, WI 53706
Equal contribution; authors in alphabetical order Leonardo Zepeda-Nuñez∗
Department of Mathematics
University of Wisconsin-Madison
Madison, WI 53706 Now at Google Research Shengwen Yang
Department of Statistics
University of Wisconsin-Madison
Madison, WI 53706
Zelin Zhao
Department of Mathematics
University of Wisconsin-Madison
Madison, WI 53706
Claudia Solís-Lemus
Wisconsin Institute for Discovery
University of Wisconsin-Madison
Madison, WI 53706
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Scientists world-wide are putting together massive efforts to understand how
the biodiversity that we see on Earth evolved from single-cell organisms at
the origin of life and this diversification process is represented through the
Tree of Life. Low sampling rates and high heterogeneity in the rate of
evolution across sites and lineages produce a phenomenon denoted “long branch
attraction" (LBA) in which long non-sister lineages are estimated to be
sisters regardless of their true evolutionary relationship. LBA has been a
pervasive problem in phylogenetic inference affecting different types of
methodologies from distance-based to likelihood-based. Here, we present a
novel neural network model that outperforms standard phylogenetic methods and
other neural network implementations under LBA settings. Furthermore, unlike
existing neural network models, our model naturally accounts for the tree
isomorphisms via permutation invariant functions which ultimately result in
lower memory and allows the seamless extension to larger trees.

_K_ eywords phylogenetic tree $\cdot$ protein sequence evolution $\cdot$ long-
branch attraction $\cdot$ recurrent network $\cdot$ tree isomorphisms

##  1 Introduction

The Tree of Life is the representation of the evolutionary process that
originated the world-wide diversity of life from single-cell organisms.
Existing approaches to reconstruct the Tree of Life involve the collection of
genomic data from living organisms, and then, the use of statistical
techniques to estimate the tree which is the best fit to the observed genomic
data under a specified model of evolution. While the cost of sequencing
collected organisms continues to decrease, many clades in the tree still do
not have a high sampling proportion due to limited or inaccessible species
ranges or poor allocation of resources. Furthermore, in addition to the 1.4
million estimated species, it is estimated that there are from 10 to 100
million undiscovered and undescribed species [1].

Low sampling of species can cause long branches in the tree which could be
incorrectly grouped in a phenomenon denoted “long branch attraction" (LBA) [2,
3, 4, 5] where long branches in the tree are estimated to be sisters
regardless of their true evolutionary relationship. More formally, LBA arises
when the probability that close relatives share character states due to common
ancestry is exceeded by the probability that more distantly related taxa share
states due to convergent evolution [2]. It has been studied that increasing
the sampling frequency can potentially divide those long branches and thus,
overcome the LBA problem [6, 7], yet it is not always possible to increase
sampling of some clades in the Tree of Life.

Aside from low sampling, LBA also arises when there are inequalities in the
rate of evolution among branches [3]. In this scenario, commonly referred to
as “Felsenstein zone" [8, 9], an incorrect tree that groups taxa belonging to
high-rate lineages will be more likely to be estimated. This implies that LBA
would appear in rapidly evolving lineages which produce long branches leading
to some taxa which are then artificially estimated as sister taxa by most
traditional phylogenetic methods. These traditional methods tend to ignore
heterogeneities such as within-site rate variation (heterotachy) [10] or rate
variation across lineages. Site-heterogeneous models [11] that allow different
rates of mutation for different sites have successfully overcome LBA for some
datasets [12]. Yet the complexity to infer site-heterogeneous models has
hampered its use for large datasets with many sites and many taxa.

The first implementations of neural network models in phylogenetic inference
[13, 14] have found that neural networks are able to accurately estimate
phylogenetic trees with long branches and thus, overcome the LBA phenomenon.
These implementations share three similarities: 1) they all focus on neural
network models (mainly CNNs); 2) they can only estimate 4-taxon trees, and 3)
they have been proven to outperform standard phylogenetic inference methods
like maximum likelihood, Bayesian inference and maximum parsimony under a
variety of heterogeneous simulated scenarios, especially those related to the
anomaly zone (LBA) [13, 14], Felsenstein trees, Farris trees [15] or scenarios
with highly gapped alignments [13].

The main shared weakness of these phylogenetic neural network models is their
limitation to 4-taxon datasets. While Zou et al [14] was able to analyze
larger real datasets in combination with Quartet Puzzling [16], extending the
neural network model to $n$ taxa would allow a direct estimation of a
$n$-taxon tree as opposed to breaking up an $n$-taxon dataset into all
combinations of 4 taxa to later merge the resulting quartets. One of the
challenges to extend these neural network models to more taxa is the symmetric
nature of the tree structure. For example, in the quartet ((A,B),(C,D)), the
sequences A and B (and the sequences C and D) are interchangeable, as are the
clades (A,B) and (C,D) with ((C,D),(A,B)) or ((C,D),(B,A)) representing the
same tree structure. Zou et al [14] account for all the symmetries by
artificially creating all 24 permutations of the sequences to guarantee that
their neural network model is invariant to the symmetries. This process,
however, is not scalable for larger trees both in terms of computing time and
memory allocation for all the artificial permutations. More recent
implementations of neural networks in phylogenetics can estimate trees with up
to 15 taxa [17], yet this model relies on comparing simulated data with real
data, rather than classifying the best tree. Other deep learning approaches,
such as the Phyloformer [18] could infer much larger trees. However, this
approach apply deep learning on constructing the pairwise distance matrix,
thus it is not a direct estimation of the tree topologies by the model itself.

Here, we exploit the inherent symmetries of the tree objects with the first
symmetry-preserving neural network model which provides a classification
framework that can be extended to more than four taxa. We show that our NN
model not only outperforms standard phylogenetic inference methods like
neighbor joining, maximum likelihood and Bayesian inference, it also
outperforms existing neural network models [14] all while reducing the memory
needs and training time by not requiring the creation of artificial
permutations of the data. Our work opens the door to neural network models
capable of directly classifying $n$-taxon trees through the invariant
transformations that are built into our model.

##  2 Description of symmetry-preserving neural network model for four taxa

The particular structure of our proposed method leverages two main features of
the problem: the positional information within the sequence and the
equivariance of the tree-topology with respect to the order of the sequences.
These two main features are then embedded within the architecture of the
network. In what follows we will describe in general terms the architecture
and how is able to preserve the two properties mentioned above. For the sake
of clarity, we will present them in a modular fashion, similarly as it is
implemented in the code.

Symmetry-preserving neural network. Even though the number of possible trees
increases super-exponentially with the number of species, there exist several
symmetries that one can leverage. Suppose that we have four different species,
denoted A, B, C, and D, and for each of these species we have data in the form
of sequences $s_{\texttt{A}}$, $s_{\texttt{B}}$, $s_{\texttt{C}}$, and
$s_{\texttt{D}}$. The task at hand is to estimate the most likely tree
topology that explains the data fed as a matrix of the form
$[s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},s_{\texttt{D}}]$. In this case
of four taxa, we only have _three_ possible unrooted tree topologies which are
depicted in Fig. 1 (top). We then seek a function
$p([s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},s_{\texttt{D}}])\rightarrow[p_{I},p_{II},p_{III}]$,
where $p_{I}$ is the probability that the data is explained by the tree of
type I (and similarly for $p_{II}$ and $p_{III}$). Even though we can obtain
$24$ different trees with four leaves by permutation of the leaf labels, many
of these are equivalent in the sense that they share the same labelled
topology, as shown in Fig. 1 (bottom).

Figure 1: Top: Possible unrooted trees for 4 species (quartets). Bottom: All
the possible permutations on the order of the leaves that are isomorphic to a
tree of type I.

This simple example already sheds light on the properties we will exploit to
reduce the search space. The order in which the input is fed is completely
arbitrary, and it is subject to several invariances. For example, if the tree
to be estimated is of type I, then exchanging the position of A and B (or C
and D) leads to the same topology (Fig. 1 top). Thus, if we were to feed our
algorithm the data ordered as
$[s_{\texttt{B}},s_{\texttt{A}},s_{\texttt{C}},s_{\texttt{D}}]$, or
alternatively as
$[s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{D}},s_{\texttt{C}}]$, we would
expect the score $p_{I}$ to be the same. In addition, if we exchange each
clade, i.e., we swap A and B by C and D, then the topology remains the same,
and so should the score $p_{I}$. In addition, if one tree is predicted to be
certain type, and then we permute the order in which the sequences are feed,
it may change the type that was predicted originally. In fact, there is a
correspondence between how the labels should be modified under any permutation
of the sequences.

Creating a network that is invariant to these transformations will reduce the
number of data points necessary for training which in turn would reduce the
computational cost. More formally, we will take advantage of the newly
developed theory of invariant networks [19], which states that any permutation
invariant function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ accepts a
representation of the form

$f(x_{1},x_{2},x_{3},...,x_{n})=\Phi\left(\sum_{i=1}^{n}\phi(x_{i})\right)\,\,\,$

where

$\displaystyle\begin{split}\phi&:\mathbb{R}\rightarrow\mathbb{R}^{m}\,\,\\\
\Phi&:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}\end{split}$

where $m$ is the embedding dimension, which is case dependent. In our study,
we found out that $m=80$ offers the best balance between performance and
scalability for our input dimension. We will exploit this theory of invariant
networks and extend it to the case of phylogenetic trees.

For example, suppose that $\phi:\mathbb{D}^{L}\rightarrow\mathbb{R}^{L\times
m}$, where $\mathbb{D}$ is the dictionary for the sequences (e.g., the set of
nucleotides or amino acids), and $L$ is the sequence length. Suppose that
$\Phi:\mathbb{R}^{L\times m}\rightarrow\mathbb{R}^{L\times m}$, then we have
that $\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{B}}))$ will be invariant if we
permute A and B, and $\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{D}}))$ will be
invariant if we permute C and D.

Then, we define the descriptor

$\begin{split}\mathcal{D}_{I}([s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},s_{\texttt{D}}])&=\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{B}}))\\\
&+\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{D}}))\end{split}$

which is clearly invariant to the transformations which are invariant for the
tree of type I.

Thus, if we have a third function $\Psi:\mathbb{R}^{L\times
m}\rightarrow[0,1]$, and we can write

$\displaystyle\begin{split}p_{I}([s_{\texttt{B}},s_{\texttt{A}},s_{\texttt{C}},s_{\texttt{D}}])&=\Psi(\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{B}}))+\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{D}})))\\\
&=\Psi(\mathcal{D}_{I}([s_{\texttt{B}},s_{\texttt{A}},s_{\texttt{C}},s_{\texttt{D}}]))\end{split}$

which is a score function that satisfies all the symmetries required. We use
the same strategy to define the scores for trees of type II and III, resulting
in

$\displaystyle
p_{II}([s_{\texttt{B}},s_{\texttt{A}},s_{\texttt{C}},s_{\texttt{D}}])$
$\displaystyle=\Psi(\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{C}}))$
$\displaystyle+\Phi(\phi(s_{\texttt{B}})+\phi(s_{\texttt{D}}))),$
$\displaystyle
p_{III}([s_{\texttt{B}},s_{\texttt{A}},s_{\texttt{C}},s_{\texttt{D}}])$
$\displaystyle=\Psi(\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{D}}))$
$\displaystyle+\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{B}}))).$

The scores can be written in a more compact fashion by aggregating them into a
vector:

$p([s_{\texttt{B}},s_{\texttt{A}},s_{\texttt{C}},s_{\texttt{D}}])=\left[\begin{array}[]{c}p_{I}\\\
p_{II}\\\
p_{III}\end{array}\right]=\left[\begin{array}[]{r}\Psi(\mathcal{D}_{I}([s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},s_{\texttt{D}}]))\\\
\Psi(\mathcal{D}_{II}([s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},s_{\texttt{D}}]))\\\
\Psi(\mathcal{D}_{III}([s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},s_{\texttt{D}}]))\end{array}\right]$
(1)

where the descriptors are given by

$\left[\begin{array}[]{c}\mathcal{D}_{I}\\\ \mathcal{D}_{II}\\\
\mathcal{D}_{III}\end{array}\right]=\left[\begin{array}[]{r}\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{B}}))+\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{D}}))\\\
\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{C}}))+\Phi(\phi(s_{\texttt{B}})+\phi(s_{\texttt{D}}))\\\
\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{D}}))+\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{B}}))\end{array}\right].$
(2)

This scoring function satisfies all the symmetries, and it only requires us to
find three functions, which will be parametrized using NNs. This technique can
be generalized to larger numbers of species seamlessly.

In what follows we provide extra details on how the descriptor $\mathcal{D}$
and the scoring function $\Psi$ are implemented.

Figure 2: Diagram of the network. We start with sequences
$[s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},s_{\texttt{D}}]$ in the left,
then we apply the same functions $\phi$ and $\Phi$, but in a different order
to obtain the descriptors $\mathcal{D}_{I}$, $\mathcal{D}_{II}$, and
$\mathcal{D}_{III}$, which are described in (b). Then each descriptor is fed
to the function $\Psi$ (shown in (a)) which outputs the score for each
particular tree topology. (a) Diagram of $\Psi$ which is described by a LSTM
followed by a dense network. (b) Diagram of the computation of a slice of
$\mathcal{D}_{II}$, relying on the functions $\alpha$, which is an embedding
layer, and $\beta$, which is a shallow dense ResNet.

Symmetry-preserving descriptors. We now provide the exact form of the
descriptors for each tree topology. The main motivation of each descriptor is
to respect both the symmetries and the ordering imposed by the sequences. The
descriptor will apply the same transformation to each site but across the four
sequences independently of the rest of the sites. In a nutshell, the $i$-th
slice of the descriptor will only depend on the $i$-th site of the input
sequences.

In particular, the descriptor take an input of the form

$\mathcal{D}_{I}:\mathbb{D}^{L\times 4}\rightarrow\mathbb{R}^{L\times m}$

by applying the same transformation to each of the slices of the input
$[s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},s_{\texttt{D}}]$ independently.

By the sake of simplicity, we will denote by a generic $s\in\mathbb{D}^{L}$
one of the four sequences mentioned before and $s_{i}\in\mathbb{D}$ its $i$-th
component (or site).

As shown in Fig. 2 (b), the first operation is an embedding layer, which will
use the map $\alpha:\mathbb{D}\rightarrow\mathbb{R}^{m}$, and it will be
applied to each site of the sequence, i.e., the function $\phi$ mentioned
above

$\left(\phi(s)\right)_{[i,:]}=\alpha(s_{i})$

where we have used Matlab notation (i.e., if $A\in\mathbb{R}^{m\times m}$ then
$A_{[i,:]}$ corresponds to the $i$-th row of $A$).

Then, we can do the same for $\Phi$, which will be written as the application
of a function $\beta:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ to each slice of
the output of $\phi$, i.e.,

$\begin{split}\left(\Phi(\phi(s)+\phi(s^{\prime}))\right)_{[i,:]}&=\beta\left(\left(\phi(s)_{[i,:]}+\phi(s^{\prime})\right)_{[i,:]}\right)\\\
&=\beta(\phi(s)_{[i,:]}+\phi(s^{\prime})_{[i,:]})\\\
&=\beta(\alpha(s_{i})+\alpha(s^{\prime}_{i})).\end{split}$

Then, we can write a generic descriptor by slice as

$\mathcal{D}_{[i,:]}=\beta(\alpha(s_{i})+\alpha(s^{\prime}_{i}))+\beta(\alpha(\tilde{s}_{i})+\alpha(\hat{s}_{i})),$

where $s,s^{\prime},\tilde{s},$ and $\hat{s}$ are a certain ordering of the
generic sequences $s_{\texttt{A}},s_{\texttt{B}},s_{\texttt{C}},$ and
$s_{\texttt{D}}$, operation depicted in Fig. 2 (b).

We add the correct ordering of the sequences to take place in the expression
above. This produces three different sequences of descriptors, each one of
dimension $L\times m$. In this case, $\alpha$ is implemented by an embedding
layer. This choice allows us to avoid the hot-encoding, which allows in turn
to feed the data to the network more efficiently. On the other hand, $\beta$
implemented by dense ResNet block [20]. One could add a deeper network for the
$\alpha$ and $\beta$ functions. However, empirically, we have found that gains
seem to be marginal.

For simplicity $\phi$ and $\Phi$ are implemented with one-dimensional
convolutions with a kernel size of one, instead a dense network. This produces
the same effect as mentioned above, but allows us to avoid any unnecessary
reshape and efficiently batch the computations.

We point out that the descriptors maintain the order in the sequence, which
will be then leveraged by the recurrent network in what follows.

Recurrent network. Given the ordered nature of the data and the descriptor
network that was built to preserve such order, the scoring function $\Psi$ is
implemented using a standard recurrent LSTM network [21] using the PyTorch
built-in models to process the sequence of each descriptor to another sequence
and then use the last element of the new sequence to compute the score as
depicted in Fig. 2 (a).

As shown in Fig. 2 (a), the LSTM network takes each of the slices of the
descriptor at each site, and the hidden states of the precedent slice, and
outputs a new descriptor that contains the information of the current slice
and the information carried over by the hidden states of the precedent slices,
while preserving some of the information in a hidden state for the next slice.
This is repeated throughout the sequence, resulting in a new sequence of
descriptors in which at each site, unlike the input, has information about all
the sites before it and thus accounting for correlation among sites in one
direction. Then we only consider the last row (or slice) of the output of the
LSTM, which is fed to a feed-forward neural network, and ultimately outputs
one single number: the score of the network to be of certain tree type.

The same process is repeated for the descriptors of each type of tree in
Equation 2 and their outputs are then concatenated, thus obtaining the scores
following Equation 1, which is shown in Fig. 2.

Loss function and optimizer. We used a typical multi class cross entropy loss,
and the optimizer used was Adam [22] with the default parameters in PyTorch
[23]. We used an exponential scheduler for the learning rate.

###  2.1 Extension to five taxa

Let $\texttt{A},\texttt{B},\texttt{C},\texttt{D}$, and E be the five taxa with
sequences $s_{\texttt{A}}$, $s_{\texttt{B}}$, $s_{\texttt{C}}$,
$s_{\texttt{D}}$, and $s_{\texttt{E}}$. While for four taxa there are 3
possible unrooted tree topologies, there are now 15 possibilities for the case
of five taxa (Fig. 3). We can group these trees by the one species in the
middle (in Fig. 3, it is taxon E) having three different topologies.

Figure 3: Three possible 5-taxon unrooted tree topologies out of the 15
possibilities. Note that we can still exploit the symmetry characteristic of
these topologies.

We can still exploit the symmetric characteristic in quartet trees when
working with quintet trees. For example, in Fig. 3, species E is in the
middle. We could simply ignore species E and represent the two binary tree
tips on both sides using $\Phi$ and $\phi$ functions as we did in the quartet
tree cases. Essentially, we would create 15 descriptors on 15 quartet trees
coming from the original 15 quintet trees excluding the middle species as in
Equation 3:

$\left[\begin{array}[]{c}\mathcal{D}_{I}\\\ \mathcal{D}_{II}\\\
\mathcal{D}_{III}\\\ \mathcal{D}_{IV}\\\ \text{...}\\\
\mathcal{D}_{XV}\end{array}\right]=\left[\begin{array}[]{r}\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{B}}))+\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{D}}))\\\
\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{C}}))+\Phi(\phi(s_{\texttt{B}})+\phi(s_{\texttt{D}}))\\\
\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{D}}))+\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{B}}))\\\
\Phi(\phi(s_{\texttt{A}})+\phi(s_{\texttt{B}}))+\Phi(\phi(s_{\texttt{C}})+\phi(s_{\texttt{E}}))\\\
\text{...}\\\
\Phi(\phi(s_{\texttt{B}})+\phi(s_{\texttt{C}}))+\Phi(\phi(s_{\texttt{D}})+\phi(s_{\texttt{E}}))\end{array}\right].$
(3)

The model components and structures are exactly the same as in Fig. 2 except
now we have 15 Descriptors and 15 $p$ values.

###  2.2 Open-source implementation

Our NN model was implemented using Python leveraging Pytorch. We implement two
version of the network, dubbed in what follows, LSTM and OptLSTM. The first
one showcases the proposed network, whereas the second is an optimized version
of the pipeline. All the scripts for the model implementation can be found in
the GitHub repository: https://github.com/crsl4/nn-phylogenetics.

##  3 Simulations

###  3.1 Simulation strategy on four taxa

We simulate protein sequences on known four-taxon trees. Then, we use the
sequences as input data and the tree topology as labels to train the neural
network model described in Section 2.

We follow the simulation procedure in [14]. Zou et al simulate protein
sequences under the CTMC model and quartet gene trees, and then train a neural
network model to classify sequences based on the quartet gene tree they came
from. A tree out of the three possible unrooted quartets was selected at
random and we simulate branch lengths under long-branch attraction cases:
trees with two long branches, each sister to a short branch and separated by a
short internal branch [24]). In this setting, we simulate the short external
branches $b=0.1,0.5,1.0$ and then set the two long branches as $a=\kappa_{1}b$
for $\kappa_{1}=2,10,40$ and the internal branch as $c=\kappa_{2}b$ for
$\kappa_{2}=0.01,0.1,1.0$ (see Fig. 4 top). These settings mimic the long-
branch attraction simulations in [14].

The quartet was rooted in the internal branch always. A sequence of length
$L=1550$ (randomly chosen length in the window $(100,3000)$ used in [14]) is
simulated following a continuous-time Markov model using PAML [25, 26]. We use
the amino acid substitution empirical model of Dayhoff [27] which was randomly
chosen among the empirical models available in PAML. For every site, a
relative evolutionary rate $r$ was factored into the branch lengths. This rate
was sampled from a uniform distribution in $(0.05,1.0)$ as in [14] and used in
all the sites (that is, no site-specific rates were used).

The raw input data consisted in four amino acid sequences of length $L$ which
are fed through an embedding layer. Unlike [14], we do not have to permute the
data to account for the symmetries. Our neural network model incorporates the
symmetries directly (as described in the model description).

For these experiments we fixed the embedding dimension $m=80$, and an hidden
dimension of the LSTM network of $20$ with a drop out probability of $0.2$, a
hidden layer between the LSTM layer and the output layer, and RELU activation
functions. We considered a batch size of $16$.

In [14], the training data consisted of 100,000 samples, yet we used only
30,000 samples for our simulations and achieved comparable prediction accuracy
(see Results). We compare the prediction accuracy of our neural network models
described in Section 2 (LSTM and OptLSTM) with four phylogenetic inference
methods: three standard phylogenetic methods 1) Neighbor-joining (NJ) under a
Dayhoff amino acid model of evolution [28, 29, 30]; 2) Maximum likelihood (ML)
with RAxML under a Dayhoff amino acid model of evolution [31]; 3) Bayesian
inference (BI) with MrBayes under a Dayhoff amino acid model of evolution with
flat priors [32, 33], and 4) a previous neural network model by Zou et al
[14]. The “LSTM" and “OptLSTM" in our results differ only in the network
structures: for “LSTM", we feed each descriptor output into an LSTM network,
resulting in 3 separate NN models, while for “OptLSTM", we run all 3 output of
descriptors in one LSTM network by reshaping the tensors. The idea for both
model is exactly the same except that “OptLSTM" is more efficient. Since the
development of this work, two new neural network models have been published
[13, 15], yet these two models were implemented for DNA sequences, not
proteins, so we were unable to compare the performance of our models with
these two neural network implementations (see Discussion).

In each experiment, we set the first 9,000 samples as the training data and
the last 1,000 samples as the test data. We train the model using all training
data from different dataset, thus including all sets of branch lengths in our
simulation. We then test our model for each testing set separately. For our
implementation, we train the data for 300 epochs and test every 10 epochs to
record the best accuracy and parameters of our model. We use the Adam
optimizer with a learning rate of 0.001 to update our model and we decay the
learning rate of each parameter group by 0.9 every 10 epochs. For the Zou NN
implementation [14], we use the same parameters, except for the loss function:
the reduction applied to the output is summed instead of averaged. For the
three standard phylogenetic methods (NJ, ML, BI), we test on the last 500
samples.

We highlight that for the standard phylogenetic methods, we use the correct
model specification. That is, we use the Dayhoff model to estimate the
evolutionary distances in NJ with phangorn package in R [34], and we specify
the substitution model to be Dayhoff for RAxML and MrBayes. This means that
the performance of these methods reflects ideal conditions that are likely
unattainable with real data for which the real model is unknown.

###  3.2 Results on four taxa

Fig. 4 (bottom) shows the performance of the six methods to estimate a 4-taxon
quartet (top) under a variety of normal and LBA scenarios. We observe that for
cases with low LBA ($a=2b,c=1.0$), all six methods perform relatively well
with testing accuracy of 100% in most cases. However, when the branch length
$a$ grows ($a=10b,40b$), the standard phylogenetic methods (NJ, ML, BI) are no
longer able to estimate the correct quartet tree with prediction accuracy of
0.0 in many of these cases. The neural network implementations, on the other
hand, are able to correctly predict 60% trees even in extreme cases such as
$a=40b,c=0.1$, vastly exceeding the naive prediction of 33%.

Among the neural network implementations, our model slightly outperforms Zou’s
[14], but while the improvement on prediction accuracy is small, our NN
implementation is far more efficient in terms of memory given that we do not
create all possible permutations of the input sequences. Using the A100-SXM4
GPU on the platform of Center for High Throughput Computing at University of
Wisconsin-Madison, our NN model used on average 11700 megabytes of memory and
3772 megabytes of GPU memory during the training and testing process, which is
18% and 50% of the memory and GPU memory used respectively by the Zou’s NN
model (62670 megabytes of memory and 7415 megabytes of GPU memory on average)
with identical training and testing samples and same number of epoches. The
efficiency also translated into much faster training with an average of 3.3
hours to train our NN implementation compared to 40.2 hours for the case of
Zou’s, only 8% of the time used by Zou’s model.

Figure 4: Comparison of testing accuracy for six methods: Neighbor-joining
(NJ), Maximum Likelihood (ML), Bayesian Inference (BI), Neural Network model
in [14] (Zou), our Neural Network models (LSTM and OptLSTM). We observe that
neural network models (Zou, LSTM, OptLSTM) outperform standard phylogenetic
inference methods (NJ, ML, BI) on LBA cases ($a=10b,a=40b$). Among those
cases, our NN outperforms Zou’s with a slightly higher testing accuracy.

Fig. 5 shows the learning dynamics of our OptLSTM model on three chosen
datasets: $c=1,b=1$ with $a=2b,10b,40b$ that represent the three levels of
LBA. In these plots, we present as dashed lines the final testing accuracies
of the standard phylogenetic methods (NJ, ML, BI) and Zou’s NN model. We also
present as a gray dotted line the expected accuracy for a naive model that
selects quartets at random (33% for 3 possible quartets).

For the case of low LBA ($a=2b$), our model reaches a testing accuracy of 100%
rapidly and all methods in this case have a testing accuracy of 100% as well.
As LBA worsens (left to right in Fig. 5), the testing accuracy decreases for
all methods, but more dramatically for the standard phylogenetic methods (NJ,
ML, BI) that have a testing accuracy of 0.0% (overlapped dashed lines at
zero). Only the two NN models (Zou’s and our OptLSTM) are able to reach a
testing accuracy greater than the naive predictor with our model slightly
outperforming Zou’s in both cases ($a=10b,40b$).

Figure 5: Testing accuracy of our NN model implementation (OptLSTM) by epoch
for three datasets of increased severity of LBA (top to bottom). Dashed lines
represent final prediction accuracy for the other phylogenetic inference
methods. For the low LBA case (left: $a=2,b=c=1$), all methods reach 100%
accuracy quite fast while for the other two cases: medium LBA (center:
$a=10,b=c=1$) and high LBA (right: $a=40,b=c=1$), all standard phylogenetic
methods (NJ, ML, BI) have a 0% accuracy and the two NN models (Zou in dashed
pink line and our OptLSTM model in solid yellow line) have a higher accuracy
than the naive random prediction (dotted gray line at 33%). We note that our
method slightly outperforms Zou’s NN.

###  3.3 Simulation strategy on five taxa

Since LBA is less of a concern on more than four taxa [35], we use the same
simulation strategy as with four taxa, but with a different setup for the
branch lengths. We simulate branch lengths by uniformly sampling from $(0,1)$.
All other simulation parameters remain the same.

The hyperparameters of the model, however, are modified to achieve better
performance. We increase the embedding dimension $m$ from 80 to 160 and the
hidden dimension of the LSTM network from 20 to 40, with 6 layers of LSTM
instead of 3. We decrease the batch size from 16 to 8. These changes would
increase the training time but could offer better model quality overall [36].

We use 10,000 samples for simulations to match the training set of quartet
trees. Since the model in [14] could only work in a quartet tree setting, we
could not conduct a comparison of performance like we did in the quartet tree
case. While novel deep learning implementations are suitable for more than
four taxa [17, 18], these approaches are not one-shot learning. That is, these
models have non-machine learning steps in the pipeline and are thus not very
comparable with our approach.

In this experiment, we take 1,000 samples from the training set to train the
model for each epoch. While in theory having more training samples would
result in better model quality, in practice, we find no significant difference
in model performance between training with the whole set and training with a
subset. Training with a subset vastly improves the training time. We test the
model performance with 50 samples randomly selected every 10 epochs. The
setting for the optimizer and learning rate is the same as the quartet tree
cases. We train the model for a total of 1300 epochs.

###  3.4 Results on five taxa

Fig. 6 shows a maximum achieved accuracy of 0.72 with the best model, as shown
by the red dashed line. We also show the result from standard phylogenetic
methods (ML and NJ) with the black and green dashed lines. It is clear that in
this case, the performance of our model is inferior to the standard
phylogenetic methods, mostly due to the weakness of the classification task.
While the performance is not as good as non-LBA cases for quartet trees, it is
worth noting that the classification for a quintet tree is a much harder task
than a quartet tree. For a quartet tree, the tree space is 3, thus the naive
prediction of the model is 33%. For a quintet tree, the tree space increases
to 15, lowering the naive prediction to 6%, marked by the blue dashed line in
Fig. 6. Thus, the model has reasonable performance given the number of classes
it is trying to classify from.

Figure 6: Testing accuracy of all 1300 epochs. The black dashed line indicates
the result of the Maximum Likelihood method; the green dashed line indicates
the result of the Neighbor-Joining method; the red dashed line indicates the
maximum accuracy achieved among all epochs; the blue dashed line indicates the
random naive prediction.

##  4 Analysis of Zika virus data

Five viral aminoacid sequences for Zika virus were downloaded using NCBI [37]
(query link in GitHub repository https://github.com/crsl4/nn-phylogenetics).
See Table 1 for more details on the sequences. We estimate the best 5-taxon
tree with our OptLSTM NN model. For this case, we train our NN model on a
heterogeneous dataset that includes all different scenarios of $a,b,c$
described in the simulations (total sample size 10000). The heterogeneous
dataset is constructed by 1000 samples from all different scenarios of $a,b,c$
dataset, and has 10000 samples in total. We compare our tree with the
estimated tree from IQ-Tree [38]. We use the default model selection option
within IQ-Tree which selected the FLU model [39] as the best fit for the data.
This model, however, is different from the model we used to train the NN model
(Dayhoff). We find that the estimated trees are quite similar regardless of
the use of a different model (Figure 7). We also find that the IQ-Tree
estimated tree has some very short branches which could represent artificial
polytomies due to lack of phylogenetic signal to resolve evolutionary
relationships.

Table 1: Zika viral aminoacid sequences downloaded from NCBI [37]. Accession | Length | Host | Country | Year
---|---|---|---|---
QIH53581 | 3423 | Homo sapiens | Brazil | 2017
ANG09399 | 3423 | Homo sapiens | Honduras | 2016
AXF50052 | 3423 | Mus Musculus | Colombia | 2016
AWW21402 | 3423 | Simiiformes | Cambodia | 2016
AYI50274 | 3423 | Macaca mulatta |  | 2015
Figure 7: Estimated phylogenetic Zika virus tree using IQ-Tree (left) and our
OptLSTM model (right). Discrepancies between the trees could be explained by
the low phylogenetic signal to resolve the artificial polytomies visible in
the IQ-Tree tree. Our tree does not have meaningful branch lenghts as we only
infer the tree topology.

##  5 Discussion

We present here the first neural network model that fully exploits the
symmetries of trees to avoid permutation operations as other neural models do
in phylogenetics [14]. We showed that our method outperforms standard
phylogenetic inference methods like distance-based (NJ), likelihood (ML) and
Bayesian (BI) in cases related to LBA, yet it lags behind in the case of five
taxa.

While it is difficult to justify the existence of a model that underperforms
compared to other phylogenetic inference methods, we argue that the novelty in
the mathematical definition of our model is still worthy of recognition. Even
when the accuracy of our model to infer 5-taxon trees is below standard
phylogenetic methods such as maximum likelihood or neighbor joining, our model
is still accurate despite the complexity of the classification task.
Furthermore, our model is more efficient in terms of memory and computing time
compared to other deep learning implementations. Last, our model is the only
existing mathematical representation of a neural network model that
intrinsically exploits the tree symmetries, property that can open the door to
more complex – yet efficient – models in the future.

For example, the main limitation of existing neural network classifiers in
phylogenetics [15, 13, 14] is that inference is possible only for 4-taxon
datasets. Unlike other neural network models, however, our model manages to
extend the inference task to 5-taxon trees given the built-in invariance to
tree permutations. Zou et al [14] account for these tree isomorphisms by
explicitly adding more samples with all the possible permutations of the
sequences ($4!=24$ permutations). This data extension imposes demanding memory
requirements which become unsustainable for more taxa in the data. Our model,
on the contrary, accounts for the tree isomorphisms internally via invariant
descriptors which is memory efficient as the number of taxa increases.

Here, we also shed light into the infeasibility of neural network classifiers
for larger trees. The limitation of all existing NN classifiers is that it is
very hard to go beyond 5-taxon trees, as our model is an optimized complete
search of the tree space. For 6-taxon trees, the tree space is 105, while for
bigger trees such as 10-taxon trees, the tree space is 2,027,025. It quickly
becomes infeasible in terms of computation time and resources for deep
Learning models to search the whole tree space. Most current deep learning
models use quartet amalgamation methods such as in Zou et al’s model [14] that
combines quartet trees to estimate larger trees. Zaharias et al [40] showed
that such two stage approach performs worse than direct estimation methods
such as ML and NJ. The limitation of deep Learning methods showed in Zaharias
et al’s work calls for models that avoid doing a complete search on the tree
space. While other approaches such as Phyloformer [18] could estimate larger
trees, it is not a direct estimation of trees from sequence data using only
the deep learning model. Recently Smith et al [17] proposed a General
Adversarial Network (GAN) that does a heuristic search approach of tree space
without visiting all possible trees. This approach is very similar to standard
phylogenetic inference methods and has achieved good results for up to 10
taxon trees. Future developments of deep learning models for phylogenetic
inference should follow this direction by doing heuristic search on the tree
space or conditional generation from a learned tree distribution under an
unsupervised setting. Another direction would be deep learning models that are
designed to work with non-Euclidean space, namely the Graph Neural Network
(GNN) [41, 42, 43, 44]. Such models would more effectively learn the tree
topologies of phylogenetic datasets and could potentially offer better
performance compared to existing classifiers, which mostly transform
phylogenetic trees into Euclidean space datasets and feed into a conventional
Deep Learning model such as CNN and RNN.

In addition to extending the model to more taxa, future work will also involve
allowing for nucleotide sequences as input in addition to proteins. Currently,
our model is only suited to handle protein sequences so that we could compare
our performance to [14]. The other two neural network models applied to
phylogenetics [15, 13] use nucleotide sequences as input, so we were unable to
compare them to our model. Thus, one future extension of our model could be
the modification of the model to allow for the set of nucleotides to be used
as the dictionary for the sequences. Furthermore, Suvorov et al [13] allow for
gaps in the alignment and proved that a neural network model is better suited
to handle highly gapped sequences (which are quite common on bacterial or
viral sequences due to difficulty of alignment) compared to standard
phylogenetic methods. Therefore, another extension of our model will be the
incorporation of gaps in the sequences to test its performance to a variety of
gap scenarios.

## Data Availability

All the scripts for the model implementation, simulations and real data
analysis can be found in the GitHub repository: https://github.com/crsl4/nn-
phylogenetics.

## Acknowledgements

This work was partially supported by the Department of Energy [DE-SC0021016 to
C. S.-L.], National Science Foundation [DEB-2144367 to CSL] and the National
Science Foundation [DMS-2012292 to L. Z.-N.].

## References

  * [1] B. Groombridge.  Global biodiversity: status of the Earth’s living resources.  World Conservation Monitoring Centre, 1992.
  * [2] M. D. Hendy and D. Penny.  A framework for the quantitative study of evolutionary trees.  Systematic zoology, 38(4):297–309, 1989.
  * [3] J. Felsenstein.  Cases in which parsimony or compatibility methods will be positively misleading.  Systematic zoology, 27(4):401–410, 1978.
  * [4] J. Bergsten.  A review of long-branch attraction.  Cladistics, 21(2):163–193, 2005.
  * [5] F. E. Anderson and D. L. Swofford.  Should we be worried about long-branch attraction in real data sets? Investigations using metazoan 18S rDNA.  Molecular phylogenetics and evolution, 33(2):440–451, 2004.
  * [6] J. Gauthier, A. G. Kluge, and T. Rowe.  Amniote phylogeny and the importance of fossils.  Cladistics, 4(2):105–209, 1988.
  * [7] J. J. Wiens.  Can incomplete taxa rescue phylogenetic analyses from long-branch attraction?  Systematic Biology, 54(5):731–742, 2005.
  * [8] J. Felsenstein and J. Felenstein.  Inferring phylogenies, volume 2.  Sinauer associates Sunderland, MA, 2004.
  * [9] J. P. Huelsenbeck and D. M. Hillis.  Success of phylogenetic methods in the four-taxon case.  Systematic Biology, 42(3):247–264, 1993.
  * [10] H. Philippe, Y. Zhou, H. Brinkmann, N. Rodrigue, and F. Delsuc.  Heterotachy and long-branch attraction in phylogenetics.  BMC evolutionary biology, 5(1):1–8, 2005.
  * [11] N. Lartillot and H. Philippe.  A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process.  Molecular biology and evolution, 21(6):1095–1109, 2004.
  * [12] N. Lartillot, H. Brinkmann, and H. Philippe.  Suppression of long-branch attraction artefacts in the animal phylogeny using a site-heterogeneous model.  BMC evolutionary biology, 7(1):1–14, 2007.
  * [13] A. Suvorov, J. Hochuli, and D. R. Schrider.  Accurate Inference of Tree Topologies from Multiple Sequence Alignments Using Deep Learning.  Syst. Biol., 69(2):221–233, March 2020.
  * [14] Z. Zou, H. Zhang, Y. Guan, and J. Zhang.  Deep Residual Neural Networks Resolve Quartet Molecular Phylogenies.  Mol. Biol. Evol., 37(5):1495–1507, May 2020.
  * [15] A. F. Leuchtenberger, S. M. Crotty, T. Drucks, H. A. Schmidt, S. Burgstaller-Muehlbacher, and A. von Haeseler.  Distinguishing Felsenstein Zone from Farris Zone Using Neural Networks.  Molecular Biology and Evolution, 37(12):3632–3641, December 2020.
  * [16] K. Strimmer and A. Von Haeseler.  Quartet puzzling: a quartet maximum-likelihood method for reconstructing tree topologies.  Molecular biology and evolution, 13(7):964–969, 1996.
  * [17] M. L. Smith and M. W. Hahn.  Phylogenetic inference using Generative Adversarial Networks.  bioRxiv, 2023.
  * [18] L. Nesterenko, B. Boussau, and L. Jacob.  Phyloformer: towards fast and accurate phylogeny estimation with self-attention networks.  bioRxiv, 2022.
  * [19] M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Poczos, R. R. Salakhutdinov, and A. J. Smola.  Deep sets.  pages 3391–3401, 2017.
  * [20] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [21] S. Hochreiter and J. Schmidhuber.  Long Short-Term Memory.  Neural Computation, 9(8):1735–1780, 11 1997.
  * [22] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In The 3rd International Conference for Learning Representations (ICLR), 2015.  arXiv:1412.6980v8.
  * [23] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [24] J. Felsenstein.  Cases in which Parsimony or Compatibility Methods will be Positively Misleading.  Syst. Biol., 27(4):401–410, December 1978.
  * [25] Z. Yang, R. Nielsen, and M. Hasegawa.  Models of amino acid substitution and applications to mitochondrial protein evolution.  Mol. Biol. Evol., 15(12):1600–1611, December 1998.
  * [26] Z. Yang.  PAML 4: phylogenetic analysis by maximum likelihood.  Mol. Biol. Evol., 24(8):1586–1591, August 2007.
  * [27] M. O. Dayhoff and R. M. Schwartz.  Chapter 22: A model of evolutionary change in proteins.  Atlas of Protein Sequence and Structure, 1978.
  * [28] E. Paradis and K. Schliep.  ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R.  Bioinformatics, 35:526–528, 2019.
  * [29] N. Saitou and M. Nei.  The neighbor-joining method: a new method for reconstructing phylogenetic trees.  Molecular biology and evolution, 4(4):406–425, 1987.
  * [30] J. A. Studier and K. J. Keppler.  A note on the neighbor-joining algorithm of Saitou and Nei.  Molecular biology and evolution, 5(6):729–731, 1988.
  * [31] A. Stamatakis.  RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies.  Bioinformatics, 30(9):1312–1313, 2014.
  * [32] J. P. Huelsenbeck and F. Ronquist.  MRBAYES: Bayesian inference of phylogenetic trees.  Bioinformatics, 17(8):754–755, 2001.
  * [33] F. Ronquist and J. P. Huelsenbeck.  MrBayes 3: Bayesian phylogenetic inference under mixed models.  Bioinformatics, 19(12):1572–1574, 2003.
  * [34] K. Schliep.  phangorn: phylogenetic analysis in R.  Bioinformatics, 27(4):592–593, 2011.
  * [35] A. Graybeal.  Is it better to add taxa or characters to a difficult phylogenetic problem?  Systematic biology, 47(1):9–17, 1998.
  * [36] N. S. Keskar, D. Mudigere, J. Nocedal, M. Smelyanskiy, and P. T. P. Tang.  On Large-Batch Training for Deep Learning: Generalization Gap and Sharp Minima.  CoRR, abs/1609.04836, 2016.
  * [37] National Center for Biotechnology Information.  https://www.ncbi.nlm.nih.gov/.  Accessed: 2021-11-20.
  * [38] L.-T. Nguyen, H. A. Schmidt, A. Von Haeseler, and B. Q. Minh.  IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies.  Molecular biology and evolution, 32(1):268–274, 2015.
  * [39] C. C. Dang, Q. S. Le, O. Gascuel, and V. S. Le.  FLU, an amino acid substitution model for influenza proteins.  BMC evolutionary biology, 10(1):1–11, 2010.
  * [40] P. Zaharias, M. Grosshauser, and T. Warnow.  Re-evaluating Deep Neural Networks for Phylogeny Estimation: The Issue of Taxon Sampling.  J Comput Biol., 29(1):74–89, January 2022.
  * [41] C. Zhang.  Learnable Topological Features for Phylogenetic Inference via Graph Neural Networks, 2023.
  * [42] W. Jin, R. Barzilay, and T. Jaakkola.  Junction Tree Variational Autoencoder for Molecular Graph Generation.  In J. Dy and A. Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2323–2332. PMLR, 10–15 Jul 2018.
  * [43] D. Jiang, Z. Wu, C.-Y. Hsieh, G. Chen, B. Liao, Z. Wang, C. Shen, D. Cao, J. Wu, and T. Hou.  Could graph neural networks learn better molecular representation for drug discovery? A comparison study of descriptor-based and graph-based models.  Journal of Cheminformatics, 13(1):12, Feb 2021.
  * [44] Y. Kwon, J. Yoo, Y.-S. Choi, W.-J. Son, D. Lee, and S. Kang.  Efficient learning of non-autoregressive graph variational autoencoders for molecular graph generation.  Journal of Cheminformatics, 11(1):70, Nov 2019.