###### Abstract.

In this article we present a remarkable concyclicity of four centroids related
to the orthocenters of the triangles $ABC,$ $BPC,$ $CPA,$ and $PAB$ of a
quadrangle $ABCP$. Furthermore, we establish a result about orthologic
triangles inscribed in a rectangular hyperbola. We present an application of
the concyclic centroids theorem when $P$ corresponds to the incenter $I$ of
$ABC$. This specialization leads to the discovery of a novel triangle
intimately related to the Fuhrmann triangle.

Concyclic Centroids and the Anti-Fuhrmann Triangle

CONCYCLIC CENTROIDS AND THE ANTI-FUHRMANN TRIANGLE

SUDHARSHAN K V, SHANTANU NENE

Abstract. In this article we present a remarkable concyclicity of four
centroids related to the orthocenters of the triangles $ABC,$ $BPC,$ $CPA,$
and $PAB$ of a quadrangle $ABCP$. Furthermore, we establish a result about
orthologic triangles inscribed in a rectangular hyperbola. We present an
application of the concyclic centroids theorem when $P$ corresponds to the
incenter $I$ of $ABC$. This specialization leads to the discovery of a novel
triangle intimately connected to the Fuhrmann triangle.

##  1\. Introduction

Let $ABC$ be a triangle, and let $P$ be a point, and let $H,H_{A},H_{B},$ and
$H_{C}$ be the orthocenters of triangles $ABC,BPC,CPA,$ and $PAB$
respectively. Let $\mathcal{H}_{P}$ be a circumrectangular hyperbola of $ABC$
passing through $P$. Let $Z$ denote the Poncelet point of $ABCP$, which is
also the center of $\mathcal{H}_{P}$. The following two lemmas are well-known.

######  Lemma 1.1.

A circumconic $\mathcal{H}$ of $ABC$ is a rectangular hyperbola if and only if
$\mathcal{H}$ contains $H$.

Proof. Suppose $\mathcal{H}$ is a rectangular hyperbola. Then, $\mathcal{H}$
passes through two points at infinity along perpendicular directions, which we
denote by $\infty_{1},\infty_{2}$. We have an equality of cross ratios

$(BA,BH;B\infty_{1},B\infty_{2})=(CH,CA;C\infty_{2},C\infty_{1}),$

since each line in the second pencil is perpendicular to the corresponding
line in the first pencil. It immediately follows that
$(BA,BH;B\infty_{1},B\infty_{2})$ $=(CA,CH;C\infty_{1},C\infty_{2})$. The
point $H$ thus lies in the conic $\mathcal{H}$.

Conversely, suppose that $\mathcal{H}$ contains $H$. Let $\infty_{1}$ be a
point at infinity on $\mathcal{H}$, and let $\infty_{2}$ be a point at
infinity in the direction perpendicular to $\infty_{1}$. Once again, from an
equality of cross ratios $(BA,BH;B\infty_{1},B\infty_{2})=$
$(CH,CA;C\infty_{2},C\infty_{1})$, we obtain an equality

$(BA,BH;B\infty_{1},B\infty_{2})=(CA,CH;C\infty_{1},C\infty_{2}),$

which implies that $\infty_{2}$ lies on $\mathcal{H}$. $\square$

######  Corollary 1.1.

Let $A,B,C$ be three points on a rectangular hyperbola $\mathcal{H}$. If the
perpendicular from $A$ to $BC$ intersects $\mathcal{H}$ again at $D$, then $D$
is the orthocenter of $ABC$.

————————————–

Keywords and phrases: Poncelet Points, Fuhrmann Triangle, Nagel Point

(2020)Mathematics Subject Classification: 51M04, 51P99

######  Lemma 1.2.

Let $\mathcal{H}_{P}$ intersect the circumcircle of $ABC$ at $H^{\prime}$.
Then, $H,H^{\prime}$ are antipodal points in $\mathcal{H}_{P}$.

Proof. Because of Lemma 1.1, $\mathcal{H}_{P}$ is the unique rectangular
hyperbola containing $H^{\prime}$. Hence, $Z$ is also the Poncelet point of
$ABCH^{\prime}$. The Poncelet point of a cyclic quadrilateral is its
anticenter (See [4], Theorem 3.4). But the anticenter of $ABCH^{\prime}$ is
the midpoint of $HH^{\prime}$. $\square$

##  2\. Preliminaries

Let $\triangle_{A},\triangle_{B}$, and $\triangle_{C}$ denote the triangles
$AH_{B}H_{C},BH_{C}H_{A}$, and $CH_{A}H_{B}$ respectively.

######  Proposition 2.1.

The triangles $\triangle_{A},\triangle_{B},\triangle_{C}$ have a common
orthocenter.

Proof. See that $AH_{C}$ and $CH_{A}$ are parallel since they both are
perpendicular to $BP$. Let $X$ be the orthocenter of $\triangle_{A}$. Then,
$XH_{B}\perp AH_{C}$, so $XH_{B}\perp CH_{A}$. From Lemma 1.1, $X$ lies on
$\mathcal{H}_{P}$. Applying Corollary 1.1, we get that $X$ is the orthocenter
of $\triangle_{C}$. A similar argument shows that $X$ is also the orthocenter
of $\triangle_{B}$. $\square$

######  Remark 2.1.

The usage of Cartesian coordinates provides straightforward proofs to the
above and many of the results in this section.

We retain the notation from the proof, and denote the common orthocenter by
$X$.

######  Corollary 2.1.

The circumcircles of $\triangle_{A},\triangle_{B},$ and $\triangle_{C}$ are
concurrent at a point on $\mathcal{H}_{P}$.

Figure 1. Concurrent Circumcircles

We can further characterize the common point of these circumcircles.

######  Proposition 2.2 (Navneel Singhal).

Let $Q$ be the isogonal conjugate of $P$, and let $Q^{*}$ denote the inverse
of $Q$ in the circumcircle of $ABC$. Then the circumcircles of
$\triangle_{A},\triangle_{B},$ and $\triangle_{C}$ pass through the midpoint
of $QQ^{*}$.

Let $P^{\prime}$ denote the antipode of $P$ in $\mathcal{H}_{P}$. We require a
Lemma by AoPS user enhanced for the proof of the proposition.

######  Lemma 2.1.

Let $X^{\prime}$ be the antipode of $X$ in $\mathcal{H}_{P}$. Then
$HX^{\prime}\parallel PP^{\prime}$.

Proof. Let $\tau$ denote the tangent line to $\mathcal{H}_{P}$ at $P$.
Pascal’s Theorem applied to $PXAHPH_{A}$ shows that the points $HX\cap\tau,$
$AP\cap XH_{A},$ and $AH\cap PH_{A}$ are collinear. Now, $AP$ and $XH_{A}$ are
both perpendicular to $BH_{C}$, so their intersection is at infinity.
Similarly, $AH$ and $PH_{A}$ are perpendicular to $BC$ and so their
intersection is also at infinity. Hence, $HX$ is parallel to $\tau$.
The line joining the midpoints of parallel chords of a conic passes through
the center. So, $PZ$ bisects the segment $HX$. But, $HXH^{\prime}X^{\prime}$
is a parallelogram with center $Z$. So, $PZ$ $\parallel HX^{\prime}$.
$\square$

Proof. (Of Proposition 2.2) The circumconic $\mathcal{H}_{P}$ is the isogonal
conjugate of the line $OQ$. Choose points $Q_{1},\dots,Q_{4}$ on the line
$OQ$, and let $P_{1},\dots,P_{4}$ denote their isogonal conjugates. Then we
have an equaliy of cross ratios $(AQ_{1},AQ_{2};AQ_{3},AQ_{4})$
$=(AP_{1},AP_{2};AP_{3},AP_{4})$ from which we obtain an equality

(1)
$(Q_{1},Q_{2};Q_{3},Q_{4})_{OQ}=(P_{1},P_{2};P_{3},P_{4})_{\mathcal{H}_{P}},$

where the first and second cross ratios are taken with respect to line $OQ$
and conic $\mathcal{H}_{P}$ respectively.
Let $Q_{\infty},M$ denote the point at infinity along $OQ$ and the midpoint of
$QQ^{*}$ respectively. The pairs $(P^{\prime},Q^{*})$ and
$(H^{\prime},Q_{\infty})$ are isogonal conjugates (see [4], Theorem 2.2). Let
$M^{*}$ denote the isogonal conjugate of $M$. Then from (1), we obtain
$(P,P^{\prime};M^{*},H^{\prime})_{\mathcal{H}_{P}}=(Q,Q^{*};M,Q_{\infty})_{OQ}=-1$.

If $P_{\infty}$ is the point at infinity along $PP^{\prime}$, the Lemma gives

$(P,P^{\prime};X^{\prime},H^{\prime})_{\mathcal{H}_{P}}\stackrel{{\scriptstyle
H}}{{=}}(P,P^{\prime},P_{\infty},Z)=-1,$

where the first equality is obtained by projecting from $H$. So, $X^{\prime}$
is the isogonal conjugate of $M$. But $X^{\prime}$ is common to the
circumcircles of $\triangle_{A},\triangle_{B},$ and $\triangle_{C}$ from Lemma
1.2. $\square$

###  2.1. On Orthologic Triangles inscribed in a Rectangular Hyperbola

We present a result about orthologic triangles which follows immediately from
our discussions above.

######  Proposition 2.3.

Let $A,B,C,D,E,F$ be six points on a rectangular hyperbola $\mathcal{H}$.
Suppose that triangles $ABC$ and $DEF$ are orthologic, and suppose further
that one of the orthology centers, $P$, lies on $\mathcal{H}$. Then the other
orthology center lies on $\mathcal{H}$ as well.

Proof. Without losing generality, assume that $DP\perp BC$, $EP\perp AC$, and
$FP\perp AB$. Then, Corollary 1.1 says that $D,E,F$ are the orthocenters of
triangles $BPC,CPA,APB$ respectively. Let $Q$ be the common orthocenter of
triangles $AEF,BFD,CDE$ (see Proposition 2.1). Then, $Q$ lies on $\mathcal{H}$
and $AQ\perp EF,BQ\perp FD,CQ\perp DE$, so $Q$ is the other orthology center
of $ABC$ and $DEF$. $\square$

##  3\. Concyclic Centroids

Let $G_{A},G_{B},G_{C},$ and $G$ denote the centroids of
$\triangle_{A},\triangle_{B},\triangle_{C},$ and $ABC$ respectively. We have a
remarkable concyclicity of the points $G,G_{A},G_{B},$ and $G_{C}$, and we
state it as a Theorem.

Figure 2. Concyclic Centroids

######  Theorem 3.1 (Concyclic Centroids).

$G_{A},G_{B},G_{C},$ and $G$ lie on a circle.

Proof. First, we claim that $GG_{A}\perp AP$. Indeed, the quadrilateral
$BH_{C}CH_{B}$ is a trapezoid, and so the line joining the midpoints $M_{1}$
and $M_{2}$ of the segments $BC$ and $H_{C}H_{B}$ is parallel to $BH_{C}$ and
$CH_{B}$, which are in turn perpendicular to $AP$. But the segment $GG_{A}$ is
obtained by scaling $M_{1}M_{2}$ about point $A$ by a factor of $\frac{2}{3}$.
So, $GG_{A}\perp AP$. Similarly, one shows $GG_{B}$ and $GG_{C}$ are
perpendicular to $BP$ and $CP$. It is now clear that $\measuredangle
G_{C}GG_{B}=\measuredangle CPB$.

Let $D_{1}$ and $E$ denote the midpoints of segments $AH_{B}$ and $BH_{A}$.
Then, $G_{A}G_{B}$ is obtained by scaling $D_{1}E$ about $H_{C}$ by a factor
of $\frac{2}{3}$. Hence, $G_{A}G_{B}\parallel D_{1}E$. Similarly, if we let
$D_{2}$ and $F$ denote the midpoints of $AH_{C}$ and $CH_{A}$, then
$G_{A}G_{C}\parallel D_{2}F$. However, $D_{1}E$ and $D_{2}F$ are lines joining
the midpoints of parallel chords of a conic ($\mathcal{H}_{P}$). Thus, the
point $Z$ is common to lines $D_{1}E$ and $D_{2}F$. So, we have an equality of
directed angles (modulo $\pi$), $\measuredangle
G_{C}G_{A}G_{B}=\measuredangle(D_{2}F,D_{1}E)=\measuredangle FZE$.

Now, let $D_{0}$ be the midpoint of $PH_{A}$. But $Z$ lies on the nine point
circle of $BPC$, which is the circumcircle of $EFD_{0}$. Hence, we have an
equality of angles

$\begin{split}\measuredangle G_{C}G_{A}G_{B}&=\measuredangle FZE\\\
&=\measuredangle FD_{0}E\\\ &=\measuredangle CPB\,\text{(scaling at $H_{A}$ by
a factor of $2$)}\\\ &=\measuredangle G_{C}GG_{B},\end{split}$

after which the theorem is immediate. $\square$

##  4\. A Special Case

Let $I$ be the incenter of $ABC$. The orthocenters of triangles $BIC,$ $CIA,$
and $AIB$ are well-studied. We shall study the orthocenters in detail in the
coming sections. We make use of the notation from the previous section, and we
replace $P,$ $P^{\prime},$ and $\mathcal{H}_{P}$ by $I,$ $I^{\prime},$ and
$\mathcal{H}_{I}$.

The Nagel point $N_{a}$ is very relevant in this setting (see [2]). In fact,
we have the following:

######  Proposition 4.1.

The common orthocenter of $\triangle_{A},\triangle_{B},$ and $\triangle_{C}$
is $N_{a}$.

Figure 3. Nagel Point is the common Orthocenter

In view of the proof of Proposition 2.3, it suffices to show that $AN_{a}$ is
perpendicular to $H_{B}H_{C}$. We require a well-known Lemma to show this.

######  Lemma 4.1.

The polar of the midpoint of $BC$ with respect to the incircle of $ABC$ is
$H_{B}H_{C}$.

Proof. Let $A^{\prime},B^{\prime},$ and $C^{\prime}$ denote the midpoints of
$BC,CA,$ and $AB$. Let $D,E,$ and $F$ denote the points where the incircle
touches $BC,CA,$ and $AB$. It suffices to show that $B^{\prime}C^{\prime}$ is
the polar of $H_{A}$, because of La Hire’s Theorem.

Let $BI$ meet $B^{\prime}C^{\prime}$ at $T$ and $A^{\prime}B^{\prime}$ at $S$.
It is well-known that $T$ lies on $DE$ (see [5], Lemma 8).
$B^{\prime}T\parallel BC$ and $CE=CD$ imply that $B^{\prime}E=B^{\prime}T$.
Similarly $B^{\prime}E=B^{\prime}S$ and hence $B^{\prime}$ is the circumcenter
of $ETS$. Since $IE\perp B^{\prime}E$, the incircle and the circumcircle of
$ETS$ are orthogonal. Thus, $T$ and $S$ are inverses of each other in the
incircle. It is also well-known that $S,I,E,D,C$ are concyclic (see [5], Lemma
8), and hence $IS\perp CS$.

This means that $CS$ is the polar of $T$ in the incircle, and $H_{A}$ lies on
$CS$. By La Hire’s theorem, polar of $H_{A}$ passes through $T=BI\cap
B^{\prime}C^{\prime}$. A similar argument shows that the polar also passes
through $CI\cap B^{\prime}C^{\prime}$. Therefore the polar is
$B^{\prime}C^{\prime}$. $\square$

Proof. (Of Proposition 4.1) Let the other tangent from $A^{\prime}$ intersect
the incircle at $D^{\prime}$. By Lemma 1.2, the points $D,D^{\prime}$ lie on
$H_{B}H_{C}$. Let $AN_{a}$ intersect the incircle first at $D_{1}$ and let
$D_{2}$ be the intersection of $AN_{a}$ with $BC$. It is well-known that
$D_{1}$ is the antipode of $D$, and that $D_{2}$ and $D$ are isotomic in $BC$
(see [5], Lemma 2). Hence, by midpoint theorem, $IA^{\prime}\parallel
D_{1}D_{2}$. But, $IA^{\prime}\perp DD^{\prime}$ since $DD^{\prime}$ is the
polar of $A^{\prime}$. Hence we obtain the perpendicularity

$AN_{a}\perp DD^{\prime}=H_{B}H_{C}.$

$\square$

###  4.1. The Fuhrmann Triangle

We present a few results about the Fuhrmann triangle (see [3]). Let
$\triangle_{a}$, $\triangle_{b}$, and $\triangle_{c}$ denote the triangles
$H_{A}BC$,$H_{B}CA$ and $H_{C}BA$ respectively.

######  Definition 4.1.

The triangle formed by circumcenters $O_{a},O_{b},$ and $O_{c}$ of
$\triangle_{a}$, $\triangle_{b}$, and $\triangle_{c}$ is called the Fuhrmann
Triangle $\triangle_{f}$ of $ABC$. The circumcircle $\Omega_{f}$ of
$\triangle_{f}$ is called the Fuhrmann Circle of $ABC$.

Let $\Omega$ be the circumcircle of $ABC$, and let $AI,BI,CI$ meet $\Omega$
again at $A_{0},B_{0},C_{0}$ respectively. It is easy to see that $A_{0}$ is
the circumcenter of $BIC$. Hence, $O_{a}$ is the reflection of $A_{0}$ in
$BC$. Similarly, $O_{b}$ and $O_{c}$ are reflections of $B_{0}$ and $C_{0}$ in
$AC$ and $AB$ respectively. Hence, $\Omega_{f}$ is the Hagge circle ([1]) of
$I$ with respect to $ABC$. All of the following four propositions follow from
the properties of the Hagge circle proved in [1]. However, we shall provide
alternative proofs for some of them (due to AoPS user enhanced).

######  Proposition 4.2.

The orthocenter $H$ lies on $\Omega_{f}$.

Proof. We shall not retain the names of points introduced in this proof. Let
$D,E,$ and $F$ denote the feet of the altitudes from $A,B,$ and $C$ in $ABC$.
Let $\Phi$ denote the negative inversion centered at $H$ swapping the pairs
$(A,D),(B,E),$ and $(C,F)$.

The circumcircle of $BHC$ is the reflection of $\Omega$ about $BC$. Since the
midpoint of arc $BC$ reflects about to $O_{a}$, it follows that $HO_{a}$ is
bisects $\angle BHC$. Under $\Phi$, the circumcircle of $BHC$ maps to line
$EF$. Consequently, $\Phi(O_{a})$ is the intersection of $HO_{a}$ with $EF$.

From the angle bisector theorem, we have

$\frac{\overline{E\Phi(O_{a})}}{\overline{F\Phi(O_{a})}}=-\frac{EH}{FH}.$

Similarly, one gets

$\frac{\overline{F\Phi(O_{b})}}{\overline{D\Phi(O_{b})}}=-\frac{FH}{DH},$

and

$\frac{\overline{D\Phi(O_{c})}}{\overline{E\Phi(O_{c})}}=-\frac{DH}{EF}.$

From Menelaus’ theorem, we conclude that the images of $O_{a},O_{b},O_{c}$
under $\Phi$ are collinear. Therefore, $H$ lies on the circumcircle
$\Omega_{f}$ of $O_{a}O_{b}O_{c}$. $\square$

######  Proposition 4.3.

The center $O_{f}$ of $\Omega_{f}$ is the reflection of $I$ in the nine point
center $N$.

Proof. Once again, we shall discard the object names introduced in the proof.
Let $A^{*}$ denote the antipode of $A$ in $\Omega$. The segments
$HA^{*},O_{a}A_{0},$ and $BC$ have the same midpoint, so in particular,
$HO_{a}A^{*}A_{0}$ is a parallelogram. It follows from $AA_{0}\perp
A_{0}A^{*}$ that $AI\perp HO_{a}$.

The reflections of $H,O_{a}$ about $BC$ lie on $\Omega$. Hence, the
perpendicular bisector $\ell_{1}$ of $HO_{a}$ passes through the reflection
$O_{1}$ of $O$ about $BC$. But, $AHO_{1}O$ is a parallelogram, so $N$ is the
midpoint of $AO_{1}$. Both $AI$ and $\ell_{1}$ are perpendicular to $HO_{a}$,
so the reflection of $I$ about $N$ lies on $\ell_{1}$. Similarly, the
reflection lies on the perpendicular bisectors of $HO_{b}$ and $HO_{c}$, from
which it is immediate that this reflection must be $O_{f}$. $\square$

######  Proposition 4.4.

$N_{a}$ is the antipode of $H$ in $\Omega_{f}$.

Proof. From the well-known property of the Nagel line, we have

$\frac{\overline{IG}}{\overline{GN_{a}}}=\frac{\overline{OG}}{\overline{GH}}=\frac{1}{2}.$

From Proposition 4.3, we have the relation

$\frac{\overline{IN}}{\overline{NO_{f}}}=\frac{\overline{ON}}{\overline{NH}}=1.$

Consequently, the lines $HN_{a}$ and $HO_{f}$ are parallel to $OI$, which
implies that $H,N_{a},$ and $O_{f}$ are collinear.

Applying Menelaus’ theorem to $IN_{a}O_{f}$ with the transversal $GNH$, we
find that

$\frac{\overline{O_{f}H}}{\overline{HN_{a}}}\cdot\frac{\overline{IG}}{\overline{GN_{a}}}\cdot\frac{\overline{IN}}{\overline{NO_{f}}}=-1.$

It follows that $\overline{HN_{a}}=-2\overline{O_{f}H}$, so that $N_{a}$ and
$H$ are antipodes in $\Omega_{f}$. $\square$

######  Theorem 4.1 (Stevanovik).

The orthocenter of $\triangle_{f}$ is $I$.

Proof. This follows from Property 3 of [1], after noticing that $I$ is the
orthocenter of $A_{0}B_{0}C_{0}$. $\square$

######  Proposition 4.5.

The Euler reflection point of $\triangle_{f}$ is $H$.

Proof. This is a direct consequence of Corollary 4.1 of [1]. $\square$

##  5\. The Anti-Fuhrmann Triangle

In this section, we introduce a new triangle, which bears very similar
structure to the Fuhrmann triangle.

######  Definition 5.1.

We call the triangle formed by circumcenters of $\triangle_{A}$,
$\triangle_{B}$, and $\triangle_{C}$ the Anti-Fuhrmann Triangle
$\triangle_{r}$ of $ABC$. The circumcircle $\Omega_{r}$ of $\triangle_{r}$ is
called the Anti-Fuhrmann Circle of $ABC$.

Figure 4. Fuhrmann and Anti-Fuhrmann Triangle

######  Theorem 5.1.

The Anti-Fuhrmann triangle $\triangle_{r}$ is the reflection of the Fuhrmann
triangle $\triangle_{f}$ about the midpoint of $HI$.

The proof uses the following simple lemma:

######  Lemma 5.1.

Let $G_{A},G_{B},G_{C},G,G_{a},G_{b},G_{c},G^{\prime}$ be the centroids of
$\triangle_{A}$, $\triangle_{B}$, $\triangle_{C}$, $ABC$, $\triangle_{a}$,
$\triangle_{b}$, $\triangle_{c}$, $H_{A}H_{B}H_{C}$ respectively. Then
$GG_{A}G_{B}G_{C}$ and $G^{\prime}G_{a}G_{b}G_{c}$ are reflections about a
point.

Proof. Let $T$ denote the centroid of the six points
$A,B,C,H_{A},H_{B},H_{C}$. Then $T$ is the midpoint of segments
$GG^{\prime},G_{A}G_{a},G_{B}G_{b},G_{C}G_{c}$. Therefore $GG_{A}G_{B}G_{C}$
and $G^{\prime}G_{a}G_{b}G_{c}$ are reflections of each other in $T$.
$\square$

Proof. (Of Theorem 5.1) Note that $I$ is the common orthocenter of
$\triangle_{a},\triangle_{b}$ and $\triangle_{c}$ . It follows from taking a
homothety at $I$ with ratio $\frac{3}{2}$ that
$\overrightarrow{O_{a}O_{b}}=\frac{3}{2}\cdot\overrightarrow{G_{a}G_{b}}$.
Similarly, since $\triangle_{A},\triangle_{B},\triangle_{C}$ share an
orthocenter (Proposition 4.1), it follows that
$\overrightarrow{O_{A}O_{B}}=\frac{3}{2}\cdot\overrightarrow{G_{A}G_{B}}$.
From Lemma 5.1, $\overrightarrow{G_{A}G_{B}}=-\overrightarrow{G_{a}G_{b}}$.
Thus, $\overrightarrow{O_{A}O_{B}}=-\overrightarrow{O_{a}O_{b}}$. Analogous
relations imply $\triangle_{r}$ is a reflection of $\triangle_{f}$ in some
point.

We proceed to show that this point is actually the midpoint of $HI$. Let
$H^{*}$ be a point such that $IH^{*}$, $O_{A}O_{a},O_{B}O_{b},O_{C}O_{c}$
share a midpoint. Then it follows that $IO_{A}\parallel H^{*}O_{a}$. From
Theorem 3.1, $G$ lies on the circumcircle of $G_{A}G_{B}G_{C}$. It is a well-
known property of the Nagel line that
$\frac{\overline{IG}}{\overline{GN_{a}}}=\frac{1}{2}$. Taking a homothety at
$N_{a}$ with ratio $\frac{3}{2}$ shows that $I$ lies on $\Omega_{r}$. Hence,
$H^{*}$ lies on $\Omega_{f}$. The claim in the proof of Theorem 3.1 implies
that $GG_{A}\perp AI$. The homothety at $N_{a}$ with ratio $\frac{3}{2}$ also
implies that $O_{A}I\perp AI$. Hence $H^{*}O_{a}\perp AI$. Therefore $H^{*}$
is a point on $\Omega_{f}$, such that $H^{*}O_{a}\perp AI$. In view of the
proof of Proposition 4.3, $HO_{a}\perp AI$. From Proposition 4.2, $H$ lies on
$\Omega_{f}$. Thus $H^{*}=H$. It immediately follows that $\triangle_{r}$ and
$\triangle_{f}$ are reflections of each other in the midpoint of $HI$.
$\square$

######  Corollary 5.1.

The points $H$ and $I$ are the orthocenter and Euler reflection point of
$\triangle_{r}$ respectively.

###  5.1. Further Properties

######  Proposition 5.1.

The center of $\Omega_{r}$ is the reflection of $O$ about $I$.

Proof. Let $O_{r}$ denote the center of $\Omega_{r}$ and let $M$ denote the
midpoint of $HI$. It follows from Theorem 5.1 that $M$ is also the midpoint of
$O_{r}O_{f}$. Then, $O_{f}IO_{r}H$ is a parallelogram. Also, as a consequence
of Proposition 4.3, $OIHO_{f}$ is a parallelogram. It follows that $O_{r}$ is
the reflection of $O$ about $I$. $\square$

Recall that the circumcircles of $\triangle_{A},\triangle_{B},$ and
$\triangle_{C}$ are concurrent at a point because of Corollary 2.1.

######  Proposition 5.2.

The common point $X$ of the circumcircles of $\triangle_{A},\triangle_{B},$
and $\triangle_{C}$ lies on the line $O_{r}H$.

Proof. Since $N$ is the midpoint of $OH$ and $I$ is the midpoint of $OO_{r}$,
we have

$IN\parallel O_{r}H.$

The point $O_{f}$ lies on $IN$ from Proposition 4.3, and so does $F_{e}$
because of Feuerbach’s theorem. Using Lemma 1.2 and Proposition 4.1, we find
that $F_{e}$ is the midpoint of $O_{r}X$. Moreover, $O_{f}$ is the midpoint of
$N_{a}H$ (Proposition 4.4), so we have

$HX\parallel O_{f}F_{e}=IN.$

It follows that $X$ lies on $O_{r}H$. $\square$

Finally, we go a little bit further than Proposition 5.1.

######  Proposition 5.3.

The reflections of the Euler line of $\triangle_{r}$ about its sidelines are
precisely the lines $H_{A}I,H_{B}I,$ and $H_{C}I$.

Proof. Using Proposition 5.2 and Corollary 5.1, we find that $X$ lies on the
Euler line of $\triangle_{r}$. Since $I$ is the Euler reflection point of
$\triangle_{r}$ (again Corollary 5.1), we find that the reflection $XH$ about
$O_{B}O_{C}$ passes through $I$. Since $X$ and $H_{A}$ are common to the
circumcircles of $\triangle_{A}$ and $\triangle_{B}$, they are also
reflections about $O_{B}O_{C}$. Consequently, the reflection of $XH$ about
$O_{B}O_{C}$ is $IH_{A}$, after which the proof is immediate. $\square$

## References

  * [1] T. Cohl, Properties of hagge circle, https://artofproblemsolving.com/community/c284651h1272116, (2016).
  * [2] J. Vonk, On the nagel line and a prolific polar triangle, http://forumgeom.fau.edu/FG2008volume8/FG200827.pdf, (2008).
  * [3] E. W. Weisstein, Fuhrmann triangle. from mathworld–a wolfram web resource, https://mathworld.wolfram.com/FuhrmannTriangle.html.
  * [4] I. Zelich, The poncelet point and its applications, https://www.academia.edu/6251353/Poncelet_point_and_its_applications_3_, (2014).
  * [5] Y. Zhao, Lemmas in euclidean geometry, http://yufeizhao.com/olympiad/geolemmas.pdf, (2007).

DEPARTMENT OF MATHEMATICS

INDIAN INSTITUTE OF SCIENCE

BANGALORE, INDIA

E-mail address<EMAIL_ADDRESS>

DEPARTMENT OF MATHEMATICS

INDIAN INSTITUTE OF TECHNOLOGY

BOMBAY, INDIA

E-mail address<EMAIL_ADDRESS>


#  _FMM_ : An R Package for Modeling Rhythmic Patterns in Oscillatory Systems

Itziar Fernández
Department of Statistics and Operations Research
Universidad of Valladolid, Valladolid, Spain
<EMAIL_ADDRESS>
Alejandro Rodríguez-Collado
Department of Statistics and Operations Research
Universidad of Valladolid, Valladolid, Spain
<EMAIL_ADDRESS>
Yolanda Larriba
Department of Statistics and Operations Research
Universidad of Valladolid, Valladolid, Spain
<EMAIL_ADDRESS>
Adrián Lamela
Department of Statistics and Operations Research
Universidad of Valladolid, Valladolid, Spain
<EMAIL_ADDRESS>
Christian Canedo
Department of Statistics and Operations Research
Universidad of Valladolid, Valladolid, Spain
<EMAIL_ADDRESS>
Cristina Rueda
Department of Statistics and Operations Research
Universidad of Valladolid, Valladolid, Spain
<EMAIL_ADDRESS>

###### Abstract

This paper is dedicated to the R package FMM which implements a novel approach
to describe rhythmic patterns in oscillatory signals. The frequency modulated
Möbius (FMM) model is defined as a parametric signal plus a gaussian noise,
where the signal can be described as a single or a sum of waves. The FMM
approach is flexible enough to describe a great variety of rhythmic patterns.
The FMM package includes all required functions to fit and explore single and
multi-wave FMM models, as well as a restricted version that allows equality
constraints between parameters representing a priori knowledge about the shape
to be included. Moreover, the FMM package can generate synthetic data and
visualize the results of the fitting process. The potential of this
methodology is illustrated with examples of such biological oscillations as
the circadian rhythm in gene expression, the electrical activity of the
heartbeat and neuronal activity.

_Keywords_ oscillation $\cdot$ rhythmic signal $\cdot$ frequency modulated
Möbius model $\cdot$ FMM $\cdot$ restricted FMM model $\cdot$ R

##  1 Introduction

Oscillations naturally occur in a multitude of physical, chemical, biological,
and even economic and social processes. Periodic signals appear, for example,
during the cell-cycle, in biological time-keeping processes, in human
heartbeats, in neuronal signals, in light emissions from certain types of
stars, or in business cycles in economics, among many others. Three features
typically describe the periodic nature of the oscillatory motion: period,
amplitude and phase. The period is the time required for one complete
oscillation. Within a period, a sum of monocomponent models characterized by
the phase and amplitude parameters can be used to describe the rhythmic
pattern of a signal (Boashash, 2016). By varying the number of monocomponents
and considering phase and amplitude parameters as fixed or variable, a large
number of rhythmic signal representations can be found.

One of the most popular representations of oscillating signals is the Fourier
decomposition (FD): a multicomponent representation with a fixed amplitude
parameter. Its monocomponent version, the cosinor model (COS) (Cornelissen,
2014), is widely used, in particular in chronobiology, with acceptable results
when a sinusoidal shape response within a period is expected. Due to its
widespread use, many software utilities are available. Particularly in R (R
Core Team, 2020), the estimation of a COS model can be performed using cosinor
(Sachs, 2014) and cosinor2 packages (Mutak, 2018). In addition, other packages
from widely differing areas of knowledge have specific functions for fitting
COS models. Such is the case of, for example, the function CATCosinor in the
CATkit package (Gierke _et al._ , 2018), which implements tools for
periodicity analysis; the function cosinor in the psych package (Revelle,
2018), dedicated to personality and psychological research; or the function
cosinor contained in a very recent package, card (Shah, 2020), which is
dedicated to the assessment of the regulation of cardiovascular physiology.
Recently, it has also been implemented in other languages such as CosinorPy, a
cosinor python package (Moskon, 2020). The COS model is easy to use and
interpret with symmetrical patterns. However, asymmetric shapes are not
captured properly by COS. When the waveform is non-sinusoidal, the use of
multiple components analysis to fit a model consisting of a sum of several
periodical functions is recommended. However, the multicomponent FD models,
developed to provide flexibility from COS, often require the use of a large
number of components resulting in serious over-fitting issues.

In recent years, alternative methods, mostly nonparametric statistical
methods, have been developed and used for analyzing rhythmicity, especially in
biological data sets. Some very popular ones, such as the JTK_Cycle (Hughes
_et al._ , 2010), wrongly assume that any underlying rhythms have symmetric
waveforms. Others, such as RAIN (Thaben and Westermark, 2014), designed to
detect more diverse wave shapes including asymmetric patterns, are not focused
on modeling but on detecting rhythmic behavior in sets of data. Thus, they are
not useful to describe the underlying oscillatory phenomena. The proliferation
of methodology in this field has been accompanied by software developments.
This is the case, for example, of the DiscoRhythm R package (Carlucci _et al._
, 2020), very recently available on Bioconductor
(https://bioconductor.org/packages/DiscoRhythm) with a web interface based on
the R Shiny platform (Chang _et al._ , 2020). This tool allows four popular
approaches to be used, including the COS model and JTK_Cycle, to discover
biological rhythmicity. Another recent example is the CircaCompare (Parsons
_et al._ , 2019), an R package implemented for modeling cosinusoidal curves by
nonlinear regression and available from GitHub
(https://github.com/RWParsons/circacompare). Also hosted on GitHub, we can
find the LimoRhyde R package (https://github.com/hugheylab/limorhyde) for the
differential analysis of rhythmic transcriptome data, based on fitting linear
models (Singer and Hughey, 2019).

Motivated by the need for a flexible, interpretable and parametric methodology
to fit rhythmic patterns, our research group recently proposed the frequency
modulated Möbius (FMM) model (Rueda _et al._ , 2019). The FMM is an additive
nonlinear parametric regression model capable of adapting to non-sinusoidal
shapes and whose parameters are easily interpretable. The single component
model has been shown to successfully fit data as diverse as circadian clock
signals, hormonal levels data or light data from distant stars. In addition,
for more complex oscillatory signals, a multicomponent model of order $m$,
denoted as FMMm, which includes $m$ single FMM components, can be used. This
is, for example, the case for describing electrocardiography (ECG) signals.
The FMMecg signal, presented in Rueda _et al._ (2021b), is defined as the
combination of five single FMM components. Another interesting area where the
FMM approach has already shown its usefulness is in electrophysiological
neuroscience. Specifically, we have proposed FMM methodology for modeling
neuronal action potential (AP) curves, oscillating signals that measure the
difference between the electrical potential inside and outside the cell (see
Rueda _et al._ , 2021c; Rodríguez-Collado and Rueda, 2021b). An FMM2 model
provides an accurate fitting for a single AP curve; whereas series of AP
curves with similar repetitive spikes can be efficiently fitted by the FMMST
model, a restricted version of the multicomponent FMM model.

In this work we introduce the FMM package, programmed in R and available from
the Comprehensive R Archive Network (CRAN) at
http://CRAN.R-project.org/package=FMM. The package implements all required
functions to fit and explore single and multicomponent FMM models, as well as
a restricted multicomponent version. In addition, the FMM package provides
functions to generate synthetic data and visualize the results of the fitted
model. Furthermore, its use is illustrated in the aforementioned applications.
The remainder of this paper is organized as follows: Section 2 provides a
brief overview of both mono and multicomponent FMM models (2.1 and 2.2,
respectively). The FMMm model with equality constraints is presented in
Section 2.3. Section 3 is dedicated to the implementation details of the FMM
package. Through simulated examples, Section 4 introduces the basic usage of
the package, including the data generation and the fitting, as well as the
visualization of the results. In Section 5, the FMM package performance is
shown through three application areas governed by oscillatory systems:
chronobiology, ECG and neuroscience. Finally, a summary is provided in Section
6.

##  2 Frequency modulated Möbius (FMM) model

All the methodological details that justify the mathematical formulation of
the FMM models are given in Rueda _et al._ 2019.

At the time point $t$, a single FMM wave is defined as
$W(t;\upsilon)=A\cos(\phi(t;\alpha,\beta,\omega))$ where
$\upsilon=(A,\alpha,\beta,\omega)^{\prime}$, $A\in\Re^{+}$ represents the wave
amplitude and,

$\phi(t;\alpha,\beta,\omega)=\beta+2\arctan(\omega\tan(\frac{t-\alpha}{2}))$
(1)

the wave phase. The phase angle $\phi$ of an FMM wave is defined using the
Möbius link (see Downs and Mardia, 2002; Kato _et al._ , 2008) rather than the
linear link function as in the COS model. The Möbius link provides much more
flexibility to describe non-sinusoidal patterns. Without loss of generality,
we assume that the time point $t\in[0,2\pi]$. Otherwise it can be transformed
into $t^{\prime}\in[t_{0},T+t_{0}]$ by $t=\frac{(t^{\prime}-t_{0})2\pi}{T}$.

Each of the four parameters of an FMM wave characterizes some aspect of a
rhythmic pattern. $A$ describes the amplitude of the signal, while $\alpha$,
$\beta$ and $\omega$ describe the wave phase. $\alpha\in[0,2\pi]$ is a
translation parameter and a wave location parameter in the real space, whereas
$\beta\in[0,2\pi]$ and $\omega\in[0,1]$ describe the wave shape. To be
precise, assuming $\alpha=0$, the unimodal symmetric waves are characterized
by values of $\beta$ close to $0$, $\pi$ or $2\pi$. When $\beta=\frac{\pi}{2}$
or $\beta=\frac{3\pi}{2}$, extreme asymmetric patterns are described.
Moreover, a value of $\omega$ close to zero describes an extreme spiked wave
and, as $\omega$ value increases, the pattern is increasingly smoother. When
$\omega=1$, a sinusoidal wave is described and the FMM model matches the COS
model where $\varphi=\beta-\alpha$ is the acrophase parameter.

Two important features of a wave are the peak and trough, defined as the
highest and lowest points above and below the rest position, respectively. In
many applications, the peak and trough times could be very useful tools to
extract practical information of a wave, since they capture important aspects
of the dynamics. These two interesting parameters can be directly derived from
the main parameters of an FMM wave as,

$\displaystyle t^{U}$ $\displaystyle=$
$\displaystyle\alpha+2\arctan(\frac{1}{\omega}\tan(-\frac{\beta}{2}))$ (2)
$\displaystyle t^{L}$ $\displaystyle=$
$\displaystyle\alpha+2\arctan(\frac{1}{\omega}\tan(\frac{\pi-\beta}{2}))$

where $t^{U}$ and $t^{L}$ denote the peak and trough times, respectively.

###  2.1 Monocomponent FMM model

Let $X(t_{i})$, $t_{1}<t_{2}<\dots<t_{n}$ be the vector of observations. The
monocomponent FMM model is defined as follows:

$X(t_{i})=M+W(t_{i};\upsilon)+e(t_{i}),\quad i=1,\dots,n$ (3)

where $M\in\Re$ is an intercept parameter describing the baseline level of the
signal, $W(t_{i};\upsilon)$ is an FMM wave, and it is assumed that the errors
$e(t_{i})$ are independent and normally distributed with zero mean and a
common variance $\sigma^{2}$.

####  2.1.1 Estimation algorithm

A two-step algorithm to estimate monocomponent FMM model parameters is
proposed. We now describe the substantial details of each stage of the
algorithm.

Step 1: Initial parameter estimation. A two-way grid search over the choice of
$(\alpha,\omega)$ parameters is performed. For each pair of $(\alpha,\omega)$
fixed values, the estimates for $M$, $A$ and $\beta$ are obtained by solving a
least square problem as detailed below.

The model for a single FMM component can be written as:

$X(t_{i})=M+A\cos(t_{i}^{*}+\varphi)+e(t_{i})$ (4)

where $t_{i}^{*}=\alpha+2\arctan(\omega\tan(\frac{t_{i}-\alpha}{2}))$,
$\varphi=\beta-\alpha$, and $e(t_{i})\sim N(0,\sigma^{2})$ for $i=1,\dots,n$.

Using trigonometric angle sum identity, the model can be rewritten as:

$X(t_{i})=M+\delta z_{i}+\gamma w_{i}+e(t_{i})$ (5)

where $\delta=A\cos(\varphi)$, $\gamma=-A\sin(\varphi)$,
$z_{i}=\cos(t_{i}^{*})$ and $w_{i}=\sin(t_{i}^{*})$.

Since $\alpha$ and $\omega$ are fixed, the estimates for $M$, $\delta$ and
$\gamma$ are obtained by minimizing the residual sum of the squares (RSS),

$RSS=\sum_{i=1}^{n}(X(t_{i})-(\hat{M}+\hat{\delta}z_{i}+\hat{\gamma}w_{i}))^{2}$
(6)

And the estimates for $M$, $A$ and $\beta$ are straightforward to derive as
follows,

$\displaystyle\hat{M}$ $\displaystyle=$
$\displaystyle\bar{X}-\hat{\delta}\sum_{i=1}^{n}z_{i}-\hat{\gamma}\sum_{i=1}^{n}w_{i}$
(7) $\displaystyle\hat{A}$ $\displaystyle=$
$\displaystyle\sqrt{\hat{\delta}^{2}+\hat{\gamma}^{2}}$ (8)
$\displaystyle\hat{\beta}$ $\displaystyle=$ $\displaystyle\alpha+\varphi$ (9)

The best combination of $(\alpha,\omega)$ values, with the lowest RSS, is
retained and the corresponding estimates are the initial parameter estimation
values.

Step 2: Optimization. In the second step, the Nelder-Mead optimization method
(Nelder and Mead, 1965) is used to obtain the final FMM parameter estimates.
Nelder-Mead is a heuristic search method commonly used to minimize or maximize
an objective function in a multidimensional space.

###  2.2 Multicomponent FMM model

A multicomponent FMM model of order $m$, denoted by FMMm, is defined as

$\displaystyle X(t_{i})=M+\sum_{J=1}^{m}W(t_{i};\upsilon_{J})+e(t_{i})$ (10)
$\displaystyle t_{1}<t_{2}<\dots<t_{n};i=1,\dots,n$

where $W(t_{i};\upsilon_{J})$, hereinafter denoted by $W_{J}(t_{i})$, is the
Jth FMM wave and,

  * •

$M\in\Re$

  * •

$\upsilon_{J}=(A_{J},\alpha_{J},\beta_{J},\omega_{J})^{\prime}\in\Re^{+}\times[0,2\pi]\times[0,2\pi]\times[0,1]$;
$J=1,\dots,m$

  * •

$\alpha_{1}\leq\alpha_{2}\leq\dots\leq\alpha_{m}\leq\alpha_{1}$

  * •

$(e(t_{1}),\dots,e(t_{n}))^{\prime}\sim N_{n}(0,\sigma^{2}I_{n})$

####  2.2.1 Model adequacy

The goodness of fit of an FMM model is measured with the $R^{2}$ statistic
that represents the proportion of the variance explained by the model out of
the total variance, that is:

$R^{2}=1-\frac{\sum_{i=1}^{n}(X(t_{i})-\hat{X}(t_{i}))^{2}}{\sum_{i=1}^{n}(X(t_{i})-\bar{X})^{2}}$
(11)

where $\hat{X}(t_{i})$ represents the fitted value at $t_{i},i=1,\dots,n$.

####  2.2.2 Estimation algorithm

An iterative backfitting algorithm is proposed to derive estimates for the FMM
parameters. Let $\hat{W}_{J}^{(k)}(t_{i})$ denote the fitted values from the
Jth FMM wave at $t_{i},i=1,...,n$ in the kth iteration. The algorithm is
structured as follows:

  1. 1

Initialize. Set $\hat{W}_{1}^{(0)}(t_{i})=\dots=\hat{W}_{m}^{(0)}(t_{i})=0$.

  2. 2

Backfitting step. For $J=1,\dots,m$, calculate

$r_{J}^{(k)}(t_{i})=X(t_{i})-\sum_{I=1}^{J-1}\hat{W}_{I}^{(k)}(t_{i})-\sum_{I=J+1}^{m}\hat{W}_{I}^{(k-1)}(t_{i})$
(12)

and fit a monocomponent FMM model to $r_{J}^{(k)}(t_{i})$ obtaining
$\hat{\alpha}_{J}^{(k)}$, $\hat{\beta}_{J}^{(k)}$, $\hat{\omega}_{J}^{(k)}$
and $\hat{W}_{J}^{(k)}(t_{i})$.

  3. 3

Repeat the backfitting step until the stopping criterion is reached. The
stopping criterion is defined as the difference between the explained
variability in two consecutive iterations: $R_{k}^{2}-R_{k-1}^{2}\leq C$,
where $R_{k}^{2}$ (defined in Equation 11) is the proportion of variance
explained by the model in the kth iteration and $C$ a constant.

  4. 4

$\hat{M}$ and $\hat{A_{J}}$ are derived by solving

$\min_{M\in\Re;A_{J}\in\Re^{+}}\sum_{i=1}^{n}(X(t_{i})-M-\sum_{J=1}^{m}A_{J}\cos(\hat{\phi}_{J}(t_{i})))^{2}$
(13)

where
$\hat{\phi}_{J}(t_{i})=\phi(t_{i};\hat{\alpha}_{J},\hat{\beta}_{J},\hat{\omega}_{J})$
defined in Equation 1.

###  2.3 Restricted multicomponent FMM model

Modeling signals with repetitive shape-similar waves can be very useful in
some applications (see Rodríguez-Collado and Rueda, 2021b). In order to obtain
more efficient estimators, equality constraints are imposed on the $\beta$ and
$\omega$ parameters of an FMMm model. In particular, we add $d$ blocks of
restrictions on $\beta$ parameters:

$\displaystyle\beta_{1}=\dots=\beta_{m_{1}}$ (14)
$\displaystyle\beta_{m_{1}+1}=\dots=\beta_{m_{2}}$ $\displaystyle\dots$
$\displaystyle\beta_{m_{d-1}+1}=\dots=\beta_{m_{d}}$

and/or $d^{\prime}$ blocks of restrictions on $\omega$ parameters:

$\displaystyle\omega_{1}=\dots=\omega_{m_{1}^{\prime}}$
$\displaystyle\omega_{m_{1}^{\prime}+1}=\dots=\omega_{m^{\prime}_{2}}$
$\displaystyle\dots$
$\displaystyle\omega_{m_{d^{\prime}-1}^{\prime}+1}=\dots=\omega_{m^{\prime}_{d^{\prime}}}$

The parameter estimation problem is solved by an adaptation of the standard
procedure.

####  2.3.1 FMMm estimation algorithm with restrictions on the $\beta$
parameters

Given the unrestricted estimates obtained in step $3$, the estimates for
$\beta_{1},\beta_{m_{1}+1},\dots,\beta_{m_{d-1}+1}$ under equality
restrictions (Equation 14) are computed as follows:

$\displaystyle\hat{\beta}_{J}^{*}=\mathrm{angularMean}(\hat{\beta}_{1},\dots,\hat{\beta}_{m_{1}})$
$\displaystyle J=1,\dots,m_{1}$
$\displaystyle\hat{\beta}_{J}^{*}=\mathrm{angularMean}(\hat{\beta}_{m_{1}+1},\dots,\hat{\beta}_{m_{2}})$
$\displaystyle J=m_{1}+1,\dots,m_{2}$ $\displaystyle\dots$
$\displaystyle\hat{\beta}_{J}^{*}=\mathrm{angularMean}(\hat{\beta}_{m_{d-1}+1},\dots,\hat{\beta}_{m_{d}})$
$\displaystyle J=m_{d-1}+1,\dots,m_{d}$

Then, the algorithm continues to the next step.

####  2.3.2 FMMm estimation algorithm with restrictions on the $\omega$
parameters

A block is defined as the set of FMM waves with common $\omega$ parameter.
When constraints for the $\omega$ parameters are incorporated, the grid search
for the different $\omega$ values is outside the backfitting loops. When the
number of blocks is large, the estimation procedure can be computationally
unaffordable. In order to reduce the execution time, a two-nested backfitting
algorithm is proposed. In the outer backfitting loop, a block is fitted. In
the inner loop, the FMM waves belonging to the same block are estimated. This
procedure generates a close to optimal solution and is a less computationally
expensive alternative.

##  3 FMM package: Implementation details

The FMM code makes use of the doParallel package (Corporation and Weston,
2020) to embed parallelization for the restricted version of the fitting
process. Several utilities from the ggplot2 (Wickham, 2016) and RColorBrewer
(Neuwirth, 2014) packages are occasionally necessary for the visualization of
the fitted models.

The implementation of FMM is divided into four main functionalities described
from Sections 3.1 to 3.4: the fitting of the FMM models, the new $S4$ object
of class ‘FMM’, the graphical visualization of the fittings and the simulation
of synthetic data.

Some general details about the functions contained in the FMM package are
shown in Table 1.

Function |  Description
---|---
Fitting function |
$\mathrm{fitFMM(vData,timePoints,nback,...)}$ |  Estimates an FMMnback model to $\mathrm{vData}$ observed at $\mathrm{timePoints}$.
Utility functions |
$\mathrm{plotFMM(objFMM,...)}$ |  Graphically displays an object of class ‘FMM’.
$\mathrm{generateFMM(M,A,alpha,beta,omega,...)}$ |  Simulates values from an FMM model.
$\mathrm{getFMMPeaks(objFMM,...)}$ |  Estimates peak and trough times, together with signal values at those times, for each FMM wave.
$\mathrm{extractWaves(objFMM)}$ |  Extracts individual contribution to the fitted values of each FMM wave.
Standard methods for objects of class ‘FMM’ |
$\mathrm{summary()}$, $\mathrm{show()}$, $\mathrm{coef()}$,
$\mathrm{fitted()}$
Table 1:  Summary of the fitting, utility functions and standard methods
implemented in FMM package.

###  3.1 Fitting an FMM model

An FMM model can be fitted using the main function fitFMM(). The description
and default values of its inputs arguments are shown in Table 2.

The fitting function fitFMM() requires the $\mathrm{vData}$ input argument,
which contains the data to be fitted. Two other arguments can be used to
control a basic fitting: $\mathrm{timePoints}$, which contains the specific
time points of the single period; and $\mathrm{nback}$, with the number of FMM
components to be fitted. For some applications, such as the study of circadian
rhythms, data are collected over multiple periods. This information is
received by the fitFMM() function through the input argument
$\mathrm{nPeriods}$. When $\mathrm{nPeriods}>1$, the FMM fitting is carried
out by averaging the data collected at each time point across all considered
periods.

Argument |  Default value |  Description
---|---|---
$\mathrm{vData}$ |  no default value |  Data to be fitted by an FMM model.
$\mathrm{nPeriods}$ |  $1$ |  Number of periods at which $\mathrm{vData}$ is observed.
$\mathrm{timePoints}$ |  $\mathrm{NULL}$ |  Vector of time points per period at which data is observed. When $\mathrm{timePoints=NULL}$ an equally spaced sequence from $0$ to $2\pi$ will be assigned.
$\mathrm{nback}$ |  $1$ |  Number of FMM components to be fitted.
$\mathrm{betaRestrictions}$ |  $1:\mathrm{nback}$ |  $\beta$ constraint vector.
$\mathrm{omegaRestrictions}$ |  $1:\mathrm{nback}$ |  $\omega$ constraint vector.
$\mathrm{maxiter}$ |  $\mathrm{nback}$ |  Maximum number of iterations for the backfitting algorithm.
$\mathrm{stopFunction}$ |  $\mathrm{alwaysFalse}$ |  Function to check the criterion convergence for the backfitting algorithm.
$\mathrm{lengthAlphaGrid}$ |  $48$ |  Grid resolution of the parameter $\alpha$.
$\mathrm{lengthOmegaGrid}$ |  $24$ |  Grid resolution of the parameter $\omega$.
$\mathrm{numReps}$ |  $3$ |  Number of times $(\alpha,\omega)$ parameters are refined.
$\mathrm{showProgress}$ |  $\mathrm{TRUE}$ |  $\mathrm{TRUE}$ to display a progress indicator on the console.
$\mathrm{showTime}$ |  $\mathrm{TRUE}$ |  $\mathrm{TRUE}$ to display execution time on the console.
$\mathrm{parallelize}$ |  $\mathrm{FALSE}$ |  $\mathrm{TRUE}$ to use parallelized procedure to fit restricted FMMm model.
Table 2:  Description of the input arguments of the fitFMM() function and
their default values.

There are three key issues in the fitting process: the grid search of the pair
$(\alpha,\omega)$ to solve the estimation problem of a single FMM wave, the
backfitting algorithm used for the estimation of the multicomponent models,
and the incorporation of restrictions on $\beta$ and $\omega$ parameters. Each
of these issues is controlled by several arguments described below.

  * •

Grid search of the pair $(\alpha,\omega)$. The $\mathrm{lengthAlphaGrid}$ and
$\mathrm{lengthOmegaGrid}$ arguments are used to set the grid resolution by
specifying the number of equally spaced $\alpha$ and $\omega$ values. Thus,
the objetive function will be evaluated a total number of
$(\mathrm{lengthAlphaGrid})\times(\mathrm{lengthOmegaGrid})$ times, so when
both arguments are large, the computational demand can be high. By reducing
the size of the sequences of the $\alpha$ and $\omega$ parameters, the
algorithm will be computationally more efficient. However, it may fail to
obtain an accurate estimation if the grid resolution is too sparse. An
implemented option to fine-tune the estimation of the parameters is to repeat
the fitting process a $\mathrm{numReps}$ of times, in such a way that, at each
repetition, a new two-dimensional grid of $(\alpha,\omega)$ points is created
around the previous estimates.

  * •

Backfitting algorithm. The argument $\mathrm{maxiter}$ sets the maximum number
of backfitting iterations. Through the argument $\mathrm{stopFunction}$, it is
possible to set a stopping criterion. Two criteria have been implemented as
stop functions in this package. When $\mathrm{stopFunction=alwaysFalse}$,
$\mathrm{maxiter}$ iterations will be forced. If $\mathrm{stopFunction=R2()}$,
the algorithm will be stopped when the difference between the explained
variability in two consecutive iterations is less than a value pre-specified
in the $\mathrm{difMax}$ argument.

  * •

Restrictions. The arguments $\mathrm{betaRestrictions}$ and
$\mathrm{omegaRestrictions}$ set the equality constraints for the $\beta$ and
$\omega$ parameters, respectively. For the unrestricted case,
$\mathrm{betaRestrictions}=1:\mathrm{nback}$ and
$\mathrm{omegaRestrictions}=1:\mathrm{nback}$. To add restrictions, ‘integer’
vectors of length $m$ can be passed to these arguments, so that positions with
the same numeric value correspond to FMM waves whose parameters, $\beta$
and/or $\omega$, are forced to be equal. Since restricted fitting can be
computationally intensive, especially for constraints on $\omega$ parameters,
a parallel processing implementation can be used when the argument
$\mathrm{parallelize=TRUE}$.

###  3.2 Object of class ‘FMM’

The fitFMM() function outputs an $S4$ object of class ‘FMM’ which contains the
slots presented in Table 3.

Slot |  Description
---|---
$\mathrm{timePoints}$ |  A ‘numeric’ vector containing the time points for each data point if one single period is observed.
$\mathrm{data}$ |  A ‘numeric’ vector containing the data to be fitted to an FMM model. Data could be collected over multiple periods.
$\mathrm{summarizedData}$ |  A ‘numeric’ vector containing the summarized data at each time point across all considered periods.
$\mathrm{nPeriods}$ |  A ‘numeric’ value containing the number of periods in $\mathrm{data}$.
$\mathrm{fittedValues}$ |  A ‘numeric’ vector of the fitted values by the FMM model.
$\mathrm{M}$ |  A ‘numeric’ value of the estimated intercept parameter $M$.
$\mathrm{A}$ |  A $m$-element ‘numeric’ vector of the estimated FMM wave amplitude parameter(s) $A$.
$\mathrm{alpha}$ |  A $m$-element ‘numeric’ vector of the estimated FMM wave phase translation parameter(s) $\alpha$.
$\mathrm{beta}$ |  A $m$-element ‘numeric’ vector of the estimated FMM wave skewness parameter(s) $\beta$.
$\mathrm{omega}$ |  A $m$-element ‘numeric’ vector of the estimated FMM wave kurtosis parameter(s) $\omega$.
$\mathrm{SSE}$ |  A ‘numeric’ value of the residual sum of squares values.
$\mathrm{R2}$ |  A $m$-element ‘numeric’ vector specifying the explained variance by each of the fitted FMM components.
$\mathrm{nIter}$ |  A ‘numeric’ value containing the number of iterations of the backfitting algorithm.
Table 3:  Summary of the slots of the $S4$ object of class ‘FMM’ resulting
from fitting an FMM model with $m$ components.

The standard methods implemented for the class ‘FMM’ include the functions
summary(), show(), coef() and fitted(). These methods display relevant
information of the FMM fitting, and provide the estimated parameters and
fitted values. In addition, two more specific functions have been implemented.
Through the extractWaves() function, the individual contribution of each FMM
wave to the fitted values can be extracted. Finally, the location of the peak
and trough of each FMM wave, as well as the value of the signal at these time
points, can be estimated using the getFMMPeaks() function. The required
argument of all these methods and functions is an object of the class ‘FMM’.
Particularly, getFMMPeaks() has an optional argument:
$\mathrm{timePointsIn2pi}$, that forces the peak and trough locations to be
returned into the interval from $0$ to $2\pi$ when it is $\mathrm{TRUE}$.

###  3.3 Plotting FMM models

The FMM package includes the function plotFMM() to visualize the results of an
FMM fit. The arguments of this function are summarized in Table 4.

Argument |  Default value |  Description
---|---|---
$\mathrm{objFMM}$ |  no default value |  The object of class ‘FMM’ to be plotted.
$\mathrm{components}$ |  $\mathrm{FALSE}$ |  $\mathrm{TRUE}$ to display a plot of components.
$\mathrm{plotAlongPeriods}$ |  $\mathrm{FALSE}$ |  $\mathrm{TRUE}$ to plot more than 1 period.
$\mathrm{use\\_ggplot2}$ |  $\mathrm{FALSE}$ |  $\mathrm{TRUE}$ to plot with ggplot2 package.
$\mathrm{legendInComponentsPlot}$ |  $\mathrm{TRUE}$ |  $\mathrm{TRUE}$ to indicate if a legend should be plotted in the component plot.
$\mathrm{textExtra}$ |  empty string |  Extra text to be added to the title of the plot.
Table 4:  Description of the input arguments of the plotFMM() function and
their default values.

An object of class ‘FMM’ can be plotted in two ways (see Figure 1). The
default graphical representation will be a plot on which original data (as
points) and the fitted model (as a line) are plotted together (left panel in
Figure 1). The other possible representation is a component plot for
displaying each centered FMM wave separately (right panel in Figure 1). Set
the boolean argument $\mathrm{components=TRUE}$ to show a component plot. When
$\mathrm{legendInComponentsPlot=TRUE}$, a legend appears at the bottom of the
component plot to indicate the represented waves. The argument
$\mathrm{textExtra}$ allows an extra text to be added to the title of both
graphical representations.

As mentioned above, in some cases, data are collected from different periods.
All periods can be displayed simultaneously on the default plot using
$\mathrm{plotAlongPeriods=TRUE}$. For the component plot, this argument is
ignored.

The argument $\mathrm{use\\_ggplot2}$ provides a choice between building the
plot using base R graphics or ggplot2 packages. By default, the graphics
package is used. When $\mathrm{use\\_ggplot2=TRUE}$, a more aesthetic and
customizable plot is created using the ggplot2 package.

###  3.4 Simulating data from an FMM model

Data from an FMM model can be easily simulated using the function
generateFMM() of the package FMM. All input arguments of this function are
shown in Table 5, along with a short description and their default values.

Argument |  Default value |  Description
---|---|---
$\mathrm{M}$ |  no default value |  Value of the intercept parameter $M$.
$\mathrm{A}$ |  no default value |  Vector of the FMM wave amplitude parameter $A$.
$\mathrm{alpha}$ |  no default value |  Vector of the FMM wave phase translation parameter $\alpha$.
$\mathrm{beta}$ |  no default value |  Vector of the FMM wave skewness parameter $\beta$.
$\mathrm{omega}$ |  no default value |  Vector of the FMM wave kurtosis parameter $\omega$.
$\mathrm{from}$ |  $0$ |  Initial time point of the simulated data.
$\mathrm{to}$ |  $2\pi$ |  Final time point of the simulated data.
$\mathrm{length.out}$ |  $100$ |  Desired length of the simulation.
$\mathrm{timePoints}$ |  $\mathrm{seq(from,to,length=length.out)}$ |  Time points at which the data will be simulated.
$\mathrm{plot}$ |  $\mathrm{TRUE}$ |  $\mathrm{TRUE}$ when the simulated data should be drawn on a plot.
$\mathrm{outvalues}$ |  $\mathrm{TRUE}$ |  $\mathrm{TRUE}$ when the numerical simulation should be returned.
$\mathrm{sigmaNoise}$ |  $0$ |  Standard deviation of the gaussian noise to be added.
Table 5:  Description of the input arguments of the generateFMM() function and
their default values.

The main arguments of this function are $\mathrm{M}$, $\mathrm{A}$,
$\mathrm{alpha}$, $\mathrm{beta}$ and $\mathrm{omega}$, whereby the values of
the FMM model parameters are passed to the function. All these arguments are
‘numeric’ vectors of length $m$, except $\mathrm{M}$, which has length $1$.
Longer and smaller vectors will be truncated or replicated as appropriate.

By default, the data will be simulated at a sequence of $100$ equally spaced
time points from $0$ to $2\pi$. The arguments $\mathrm{from}$, $\mathrm{to}$
and $\mathrm{length.out}$ control such sequences. The sequence can also be
manually set using the argument $\mathrm{timePoints}$, in which case
$\mathrm{from}$, $\mathrm{to}$ and $\mathrm{length.out}$ will be ignored.

The user can add a gaussian noise by argument $\mathrm{sigmaNoise}$. A
positive ‘numeric’ value sets the corresponding standard deviation of the
gaussian noise to be added. To create the normally distributed noise, the
rnorm() function is used.

The arguments $\mathrm{plot}$ and $\mathrm{outvalues}$, both boolean values,
determine the output of the generateFMM() function. When
$\mathrm{outvalues=TRUE}$, a ‘list’ with input parameters, time points and
simulated data is returned. These elements are named $\mathrm{input}$,
$\mathrm{t}$ and $\mathrm{y}$, respectively. In addition, a scatter plot of
$\mathrm{y}$ against $\mathrm{t}$ can be drawn by setting
$\mathrm{plot=TRUE}$.

##  4 Basic usage of the FMM package

The two examples bellow, based on FMM synthetic data, illustrate the basic
uses and capabilities of the functions implemented in the FMM package. The
first example is based on a simulated data set from a multicomponent FMM model
and illustrates the fitting process and the presentation of its results. The
second is devoted to the restricted version of FMMm models and serves to
exemplify the control of the grid search related arguments to reduce the
execution time.

###  4.1 Example 1: A multicomponent FMM data set

The first simulated data set consists of $100$ observations from an FMM2 model
with intercept parameter $M=0$ and two FMM waves defined by the parameters
$\upsilon_{1}=(A_{1}=2,\alpha_{1}=1.5,\beta_{1}=0.2,\omega_{1}=0.1)^{\prime}$
and
$\upsilon_{2}=(A_{2}=2,\alpha_{2}=3.4,\beta_{2}=2.3,\omega_{2}=0.2)^{\prime}$,
respectively. The standard deviation of the error term is set at $\sigma=0.3$.
We use the function generateFMM() to simulate this data set. A set.seed()
statement is used to guarantee the reproducibility of the results.

⬇

R> library("FMM")

R> set.seed(15)

R> fmm2.data <- generateFMM(M = 0, A = rep(2,1),

alpha = c(1.5, 3.4),

beta = c(0.2, 2.3),

omega = c(0.1, 0.2),

plot = FALSE, outvalues = TRUE, sigmaNoise = 0.3)

R> str(fmm2.data)

List of 3

$ input:List of 5

..$ M : num 0

..$ A : num [1:2] 2 2

..$ alpha: num [1:2] 1.5 3.4

..$ beta : num [1:2] 0.2 2.3

..$ omega: num [1:2] 0.1 0.2

$ t : num [1:100] 0 0.0635 0.1269 0.1904 0.2539 …

$ y : num [1:100] 1.26 2.19 2.08 3.03 3.44 …

The estimation of an FMM2 can be performed by setting $\mathrm{nback=2}$ in
the fitting function fitFMM().

⬇

R> fit.fmm2 <- fitFMM(vData = fmm2.data$y, timePoints = fmm2.data$t,

nback = 2 , showTime = FALSE)

|————————————————–|

|==================================================|

Stopped by reaching maximum iterations ( 2 iterations )

When the argument $\mathrm{showProgress=TRUE}$ (its default value), a progress
indicator and information about the stopping criterion of the backfitting
algorithm will be displayed on the console. By default, $\mathrm{maxiter}$
iterations will be forced, but we can use the argument $\mathrm{stopFunction}$
to modify the criteria. For example, to continue the search until there is an
improvement, in terms of explained variability, of less than $1\%$, the R2()
function can be used specifying the $\mathrm{difMax=0.01}$ argument.

⬇

R> fit.fmm2 <- fitFMM(vData = fmm2.data$y, timePoints = fmm2.data$t, nback =
2,

showTime = FALSE, maxiter = 5, stopFunction = R2(difMax = 0.01))

|————————————————–|

|==================================================|

Stopped by the stopFunction ( 4 iterations )

The results are displayed by the function summary():

⬇

R> summary(fit.fmm2)

Title:

FMM model with 2 components

Coefficients:

M (Intercept): -0.2188

A alpha beta omega

FMM wave 1: 1.9615 3.3743 2.2275 0.2053

FMM wave 2: 2.2006 1.4930 0.2474 0.0976

Peak and trough times and signals:

t.Upper Z.Upper t.Lower Z.Lower

FMM wave 1: 0.4339 3.9237 5.7250 0.0146

FMM wave 2: 5.9656 0.1348 4.6103 -4.0110

Residuals:

Min. 1st Qu. Median Mean 3rd Qu. Max.

-0.80445 -0.21776 -0.01272 0.00000 0.17982 0.66971

R-squared:

Wave 1 Wave 2 Total

0.6907 0.2778 0.9685

The FMM wave parameter estimates, as well as the peak and trough times,
together with the signal values at those times, are presented in tabular form,
where each row corresponds to a component and each column to an FMM wave
parameter. As part of the summary, a brief description of the residuals, the
proportion of variance explained by each FMM component and by the global model
are also shown. From the above summary results, we can see that the variance
explained by the fitted model is $96.85\%$ and that the FMM waves are labelled
in decreasing order according to the $R^{2}$ value: the explained variance is
$69.07\%$ and $27.78\%$ by “FMM wave 1” and “FMM wave 2”, respectively. The
$\mathrm{summary}$ output can be assigned to an object to get a ‘list’ of all
the displayed results.

Other options to return the results are the functions coef(), getFMMPeaks()
and resid(). The first two return a ‘list’ similar to those obtained with
summary(). The resid() method can be used to obtain the complete residuals
vector. In addition, the fitted values can be extracted by the function
fitted(), which returns a ‘data.frame’ with two columns: time points and
fitted values.

⬇

R> fitted.fmm2 <- fitted(fit.fmm2)

R> head(fitted.fmm2)

timePoints fitted

1 0.00000000 1.251937

2 0.06346652 1.726570

3 0.12693304 2.269928

4 0.19039955 2.826444

5 0.25386607 3.319683

6 0.31733259 3.681024

The code given below creates Figure 1.

⬇

R> library("RColorBrewer")

R> library("ggplot2")

R> library("gridExtra")

R> titleText <- "Example 1: two simulated FMM waves"

R> defaultFMM2 <- plotFMM(fit.fmm2, use_ggplot2 = TRUE ,textExtra = titleText)
\+

theme(plot.margin=unit(c(1,0.25,1.3,1), "cm")) \+

ylim(-5, 5)

R> compFMM2 <- plotFMM(fit.fmm2, components = TRUE, use_ggplot2 = TRUE ,

textExtra = titleText) \+

theme(plot.margin=unit(c(1,0.25,0,1), "cm")) \+

ylim(-5, 5) \+

scale_color_manual(values = brewer.pal("Set1",n = 8)[3:4])

R> grid.arrange(defaultFMM2, compFMM2, nrow = 1)

The FMM plots can be written in the R graphics or ggplot2 packages. In this
example, we use $\mathrm{use\\_ggplot2=TRUE}$ to build Figure 1 based on
ggplot2. The use of ggplot2 makes it easier to customize our plots and modify
features, such as scales, margins, axes, etc. In Figure 1, the two possible
FMM plots are arranged via the $\mathrm{grid.arrange}$ function of the
gridExtra package (Auguie, 2017).

Figure 1: Graphical representations of the FMM2 fitted model from the
simulated data set of Example 1. The default plot is displayed on the left and
the component plot is shown on the right.

###  4.2 Example 2: A restricted multicomponent FMM data set

A set of $100$ observations is generated from an FMM4 model with intercept
parameter $M=3$, amplitude parameters: $A_{1}=4.5$, $A_{2}=3$, $A_{3}=1$ and
$A_{4}=1.5$, and phase translation parameters: $\alpha_{1}=1.5$,
$\alpha_{2}=4.2$, $\alpha_{3}=2$ and $\alpha_{4}=4.7$. The symmetry parameters
are equal in all four FMM waves and, with regard to the kurtosis, pairs of
waves are equal. Specifically, the shape parameters satisfy:

$\displaystyle\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=3$
$\displaystyle\omega_{1}=\omega_{2}=0.01$
$\displaystyle\omega_{3}=\omega_{4}=0.15$

Using the function generateFMM():

⬇

R> set.seed(1115)

R> rfmm.data <- generateFMM(M = 3, A = c(4.5,3,1,1.5), alpha =
c(1.5,4.2,2,4.7),

beta = rep(3,4), omega = c(rep(0.01,2),rep(0.15,2)),

plot = FALSE, outvalues = TRUE,

sigmaNoise = 0.3)

The $\mathrm{betaRestrictions}$ and $\mathrm{omegaRestrictions}$ parameters
allow a wide variety of shape restrictions to be incorporated into the fitting
procedure. In this example, to impose the shape restrictiction on the fitting
process, we use $\mathrm{betaRestrictions=c(1,1,1,1)}$ and
$\mathrm{omegaRestrictions=c(1,1,2,2)}$.

⬇

R> fit.rfmm <- fitFMM(vData = rfmm.data$y, timePoints = rfmm.data$t, nback =
4,

betaRestrictions = c(1, 1, 1, 1),

omegaRestrictions = c(1, 1, 2, 2),

showTime = FALSE)

|————————————————–|

|==================================================|

Stopped by reaching maximum iterations ( 4 iterations )

R> summary(fit.rfmm)

Title:

FMM model with 4 components

Coefficients:

M (Intercept): 1.4781

A alpha beta omega

FMM wave 1: 3.9422 1.4806 2.9996 0.0126

FMM wave 2: 1.5654 4.6953 2.9996 0.1366

FMM wave 3: 1.0413 1.9800 2.9996 0.1366

FMM wave 4: 1.8896 4.2109 2.9996 0.0126

Peak and trough times and signals:

t.Upper Z.Upper t.Lower Z.Lower

FMM wave 1: 4.6240 1.3097 4.2721 -6.8088

FMM wave 2: 1.5732 -3.7161 5.6553 -6.3382

FMM wave 3: 5.1411 -4.7470 2.9400 -6.7395

FMM wave 4: 1.0711 -2.5415 0.7192 -6.6951

Residuals:

Min. 1st Qu. Median Mean 3rd Qu. Max.

-0.735225 -0.217362 0.004187 0.000000 0.186612 0.928811

R-squared:

Wave 1 Wave 2 Wave 3 Wave 4 Total

0.3664 0.2765 0.1656 0.1081 0.9166

Figure 2 displays the two graphical representations of the fitted model: the
fitted FMM model and the component FMM plot on the left and right,
respectively.

Figure 2: Graphical representation of the fitted restricted FMM4 model with
$\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}$, $\omega_{1}=\omega_{3}$ and
$\omega_{2}=\omega_{4}$ constraints. A scatter plot of the simulated data
along with the fitted model is displayed on the left (default plot). The
component plot is shown on the right.

The fitting of more complex signals requires more execution time. By combining
the values of the arguments $\mathrm{lengthAlphaGrid}$,
$\mathrm{lengthOmegaGrid}$ and $\mathrm{numReps}$ of fitFMM(), a balance
between the computational cost and the accuracy of the estimates can be
achieved. In this example, we managed to reduce the execution time by setting
half-length sequences of the default ones, for both the
$\mathrm{lengthAlphaGrid}$ and $\mathrm{lengthOmegaGrid}$ arguments, while
increasing the number of repetitions from $3$ to $5$. Note that this
configuration even allows us to slightly improve the accuracy of the solution
from $R^{2}=91.66\%$ to $R^{2}=92.03\%$.

⬇

R> fit.rfmm.eff <- fitFMM(vData = rfmm.data$y, timePoints = rfmm.data$t, nback
= 4,

betaRestrictions = c(1, 1, 1, 1),

omegaRestrictions = c(1, 1, 2, 2),

showTime = FALSE,

lengthAlphaGrid = 24, lengthOmegaGrid = 12, numReps = 5)

|————————————————–|

|==================================================|

Stopped by reaching maximum iterations ( 4 iterations )

R> summary(fit.rfmm.eff)

Title:

FMM model with 4 components

Coefficients:

M (Intercept): 0.9556

A alpha beta omega

FMM wave 1: 3.5811 1.4855 3.0122 0.0136

FMM wave 2: 1.4826 4.6988 3.0122 0.1761

FMM wave 3: 1.0190 2.0158 3.0122 0.1761

FMM wave 4: 1.9404 4.2063 3.0122 0.0136

Peak and trough times and signals:

t.Upper Z.Upper t.Lower Z.Lower

FMM wave 1: 4.6289 0.6142 4.2143 -6.8739

FMM wave 2: 1.5800 -3.9989 5.4040 -5.5034

FMM wave 3: 5.1802 -4.9250 2.7210 -6.7180

FMM wave 4: 1.0664 -2.3001 0.6519 -6.7150

Residuals:

Min. 1st Qu. Median Mean 3rd Qu. Max.

-0.65126 -0.20547 -0.02922 0.00000 0.20787 0.91315

R-squared:

Wave 1 Wave 2 Wave 3 Wave 4 Total

0.3656 0.2605 0.1822 0.1120 0.9203

##  5 Real data analysis using the FMM package

This section illustrates the use of the FMM package on the analysis of real
signals from chronobiology, electrocardiography and neuroscience. To do this,
the package includes four real-world data sets in $\mathrm{RData}$ format
which are described in the following sections.

###  5.1 Example 1: Chronobiology

Chronobiology studies ubiquitous daily variations found in nature and in many
aspects of the physiology of human beings, such as blood pressure or hormone
levels (Mermet _et al._ , 2017). These phenomena commonly display signals with
oscillatory patterns that repeat every 24 hours, usually known as circadian
rhythms. In particular, circadian gene expression data have been deeply
analyzed in the literature as they regulate the vast majority of molecular
rhythms involved in diverse biochemical and cellular functions, see among
others Zhang _et al._ (2014), Cornelissen (2014) and Larriba _et al._ (2020).

The FMM package includes a data set called $\mathrm{mouseGeneExp}$ that
contains expression data of the Iqgap2 gene from mouse liver. The liver
circadian database is widely extended in chronobiology since the liver is a
highly rhythmic organ with moderate levels of noise (Anafi _et al._ , 2017;
Larriba _et al._ , 2018, 2020). The complete database is freely available at
NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/), with GEO accession number
GSE11923. Gene expression values are given along 48 hours with a sampling
frequency of 1hour/2days. Hence, data are collected along two periods, and an
FMM1 model is fitted to the Iqgap2 average expressed values as follows:

⬇

R> data("mouseGeneExp", package = "FMM")

R> fitGene <- fitFMM(mouseGeneExp, nPeriods = 2, showProgress = FALSE)

R> summary(fitGene)

Title:

FMM model with 1 components

Coefficients:

M (Intercept): 10.1508

A alpha beta omega

FMM wave 1: 0.4683 3.0839 1.5329 0.0816

Peak and trough times and signals:

t.Upper Z.Upper t.Lower Z.Lower

FMM wave 1: 0.1115 10.6191 6.0686 9.6825

Residuals:

Min. 1st Qu. Median Mean 3rd Qu. Max.

-9.751e-02 -3.490e-02 2.269e-03 -1.530e-06 2.670e-02 1.890e-01

R-squared:

[1] 0.8752

The behavior of the FMM versus COS model to describe this asymmetric pattern
has been compared in terms of $R^{2}$. The FMM model clearly outperforms the
COS one with an $R^{2}$ of $0.8752$ and $0.2835$, respectively. In addition, a
difference of $4.73$ hours in peak time estimation between both models is
observed, the FMM peak estimate being much more reliable, as is shown in
Figure 3.

Figure 3: Iqgap2 gene expression data along two periods (grey dots); FMM (red
line) and COS (blue line) fitted models.

###  5.2 Example 2: Electrocardiography

ECG records the periodic electrical activity of the heart. This activity
represents the contraction and relaxation of the atria and ventricle,
processes related to the crests and troughs of the ECG waveform. Heartbeats
are decomposed into five fundamental waves, labelled as $P$, $Q$, $R$, $S$ and
$T$, corresponding to the different phases of the heart’s electric activity.
The main features used in medical practice for cardiovascular pathology
diagnosis are related to the location and amplitudes of these waves, and, of
them, those labeled as $P$, $R$ and $T$ are of particular interest (Bayes de
Luna, 2007). Standard ECG signals are registered using twelve leads,
calculated from different electrode locations, Lead II being the reference
one.

The FMM package includes the analysis of a typical ECG beat from the QT
database (Laguna _et al._ , 1997). This recording, from the subject sel100,
belongs to the Normal category, regarding Physionet’s pathology classification
(Goldberger _et al._ , 2000). The data illustrate the voltage of the heart’s
electric activity, measured in $mV$, along the beat with a sampling frequency
of $250Hz$. Specifically, the ECG signal from Lead II in the fifth of the
thirty annotated beats is analysed. Recordings are publicly available on
(http://www.physionet.org). Data are saved as $\mathrm{ecgData}$ in the
package. For an ECG beat, an FMMecg, a fifth order multicomponent FMM model
can be fitted with the instruction:

⬇

R> data("ecgData", package = "FMM")

R> fitEcg <- fitFMM(ecgData, nback = 5, showProgress = FALSE)

R> summary(fitEcg)

Title:

FMM model with 5 components

Coefficients:

M (Intercept): 5.2908

A alpha beta omega

FMM wave 1: 0.6387 5.5134 3.2786 0.0322

FMM wave 2: 0.0995 4.4191 3.7638 0.1365

FMM wave 3: 0.2412 5.3504 0.6575 0.0320

FMM wave 4: 0.3310 5.5906 4.7894 0.0122

FMM wave 5: 0.0668 1.7998 2.1309 0.1641

Peak and trough times and signals:

t.Upper Z.Upper t.Lower Z.Lower

FMM wave 1: 2.3674 6.2328 3.2488 4.7246

FMM wave 2: 1.1898 4.9486 2.0805 4.6806

FMM wave 3: 2.3959 6.0793 2.1870 4.5543

FMM wave 4: 2.4227 5.7845 2.4715 4.7177

FMM wave 5: 5.1224 4.8646 4.3649 4.7226

Residuals:

Min. 1st Qu. Median Mean 3rd Qu. Max.

-0.0694805 -0.0096146 -0.0001619 0.0000000 0.0098442 0.0625571

R-squared:

Wave 1 Wave 2 Wave 3 Wave 4 Wave 5 Total

0.7650 0.0927 0.0559 0.0402 0.0380 0.9919

It is worth noting that the FMM package not only provides ECG signal-fitting
(the left hand panel in Figure 4), but it also does wave decomposition and
fiducial mark annotations on the desired waves (the right hand panel in Figure
4). It is clearly visible how the specific shapes of the five main waves
contribute to drawing and explaining the Lead II ECG waveform from the Normal
morphology. See Rueda _et al._ (2021b) for a complete review of FMMecg.

Figure 4: FMMecg performance on a single beat from patient sel100 from the QT
database. Left: Data (grey dots) and FMM fitting (red line). Black dots locate
the $P$, $R$ and $T$ fiducial marks. Right: ECG decomposition on $P$(orange),
$Q$ (purple), $R$ (green), $S$ (yellow) and $T$ (blue) waves. Dash lines
indicate $P$, $R$ and $T$ peak times.

###  5.3 Example 3: Neuroscience

####  5.3.1 Single AP curve

The study of the electrophysiological activity of neurons is one of the main
research branches in neuroscience. The AP curves are oscillatory signals that
serve as basic information units between neurons. They measure the electrical
potential difference between inside and outside the cell due to an external
stimulus. Gerstner _et al._ (2014) can serve as a basic reference for
electrophysiological neuroscience. Recently, the shape and other features of
the AP have been used in problems such as spike sorting (Rácz _et al._ , 2020;
Souza _et al._ , 2019; Caro-Martín _et al._ , 2018) or neuronal cell type
classification (Teeter _et al._ , 2018; Gouwens _et al._ , 2019; Mosher _et
al._ , 2020; Rodríguez-Collado and Rueda, 2021a).

The package includes an example of a neuronal AP. The data were simulated with
the renowned Hodgkin-Huxley model, first presented in Hodgkin and Huxley
(1952), which is defined as a system of ordinary differential equations and
has been used in a wide array of applications, as it successfully describes
the neuronal activity in various organisms. The simulation has been done using
a modified version of the python package NeuroDynex available at Gerstner _et
al._ (2014). More concretely, a short square stimulus of $12\mu A$ has been
applied to the neuron. The data can be accurately fitted by an FMM2 model as
follows:

⬇

R> data("neuronalSpike", package = "FMM")

R> fitSingleAP <- fitFMM(neuronalSpike, nback = 2, showProgress = FALSE)

R> summary(fitSingleAP)

Title:

FMM model with 2 components

Coefficients:

M (Intercept): 44.9448

A alpha beta omega

FMM wave 1: 52.9023 4.4160 3.0606 0.0413

FMM wave 2: 18.5044 4.6564 4.9623 0.0322

Peak and trough times and signals:

t.Upper Z.Upper t.Lower Z.Lower

FMM wave 1: 1.2777 110.8354 5.9670 -2.5003

FMM wave 2: 1.4320 36.9015 1.5649 -16.2606

Residuals:

Min. 1st Qu. Median Mean 3rd Qu. Max.

-14.3012 -1.0033 0.7473 0.0000 1.3230 24.8616

R-squared:

Wave 1 Wave 2 Total

0.9064 0.0604 0.9669

Figure 5: Neuronal AP simulated with the Hodgkin-Huxley model (parameters:
$C=1,g_{Na}=260,g_{K}=30,g_{L}=0.31,V_{K}=-12,V_{Na}=115,V_{L}=10.6,\widetilde{a}_{n}=1.15,\widetilde{b}_{n}=0.85,\widetilde{a}_{m}=0.9,\widetilde{b}_{m}=1.3,\widetilde{a}_{h}=1,\widetilde{b}_{h}=1$
and applying a current of $12\mu A$ for $1$ millisecond) and the estimated
FMM2 signal in red. An FD model of the same number of degree of freedom has
been fitted and plotted in blue.

The goodness of fit of the FMM2 model can be ascertained in Figure 5. For
comparison purposes, an FD model has been fitted with the same number of
degrees of freedom. While the FD attains an $R^{2}=0.3926$, the FMM model
achieves a better fit with $R^{2}=0.9669$.

####  5.3.2 AP train

Multiple AP curves, denominated spike or AP train, are usually observed as the
response to a stimulus. Various models, such as the widely used leaky-and-fire
models (Lynch and Houghton, 2015), cut the signal into segments, each one
containing an AP curve. Some authors suggest cutting the signal into even
segments (Gerstner _et al._ , 2014). However, the length of the segments turns
out to be significantly different between different types of neurons, as
explained in Teeter _et al._ (2018), and unequal data segments can lessen the
utility of some approaches. An important aspect to take into account is that
the shape of the APs in the spike train is considered to be similar and,
consequently, a restricted FMM model can accurately fit the entire signal.

The FMM package includes the data of a spike train composed of three AP
curves. The proposed model for use with these data is an FMMST model, as
defined in Rodríguez-Collado and Rueda (2021b). Each AP is modeled by two
components. The $\beta$ and $\omega$ parameters are constrained between AP
curves. The code below fits the model.

⬇

R> data("neuronalAPTrain", package = "FMM")

R> nAPs <- 3; restriction <- c(rep(1,nAPs),rep(2,nAPs))

R> fitAPTrain<-fitFMM(neuronalAPTrain, nback = nAPs*2,

betaRestrictions = restriction,

omegaRestrictions = restriction,

showProgress = FALSE)

R> summary(fitAPTrain)

Title:

FMM model with 6 components

Coefficients:

M (Intercept): 135.4915

A alpha beta omega

FMM wave 1: 51.1192 4.2327 2.8190 0.0384

FMM wave 2: 51.7219 6.1358 2.8190 0.0384

FMM wave 3: 52.1105 1.7536 2.8190 0.0384

FMM wave 4: 19.6488 0.0975 4.8644 0.0551

FMM wave 5: 19.2020 1.9978 4.8644 0.0551

FMM wave 6: 20.3642 4.4783 4.8644 0.0551

Peak and trough times and signals:

t.Upper Z.Upper t.Lower Z.Lower

FMM wave 1: 1.1036 111.2685 0.6274 -0.1787

FMM wave 2: 3.0067 111.4205 2.5304 -1.2404

FMM wave 3: 4.9077 111.7226 4.4314 -1.4963

FMM wave 4: 3.1109 58.2307 3.3336 -13.7196

FMM wave 5: 5.0112 58.2938 5.2339 -13.2943

FMM wave 6: 1.2086 58.5306 1.4313 -14.2673

Residuals:

Min. 1st Qu. Median Mean 3rd Qu. Max.

-15.1339 -1.5071 0.4973 0.0000 1.6303 19.9627

R-squared:

Wave 1 Wave 2 Wave 3 Wave 4 Wave 5 Wave 6 Total

0.3502 0.2922 0.2485 0.0411 0.0274 0.0244 0.9838

In Figure 6, the fit of the FMMST model can be visualized. The goodness of fit
of the model is excellent, achieving an $R^{2}=0.9838$.

Figure 6: Neuronal AP simulated with the Hodgkin-Huxley model (parameters:
$C=1,g_{Na}=232,g_{K}=45,g_{L}=0.215,V_{K}=-12,V_{Na}=115,V_{L}=10.6,\widetilde{a}_{n}=0.95,\widetilde{b}_{n}=1.3,\widetilde{a}_{m}=1,\widetilde{b}_{m}=1.15,\widetilde{a}_{h}=1,\widetilde{b}_{h}=1$
and applying a short square current of 4.5 $\mu A$ for $1$ millisecond) and
the estimated FMMST signal in red. The components plot of the model can be
seen on the right hand side of the figure.

##  6 Summary

A general overview on the R package FMM, which implements the estimation of
FMM models, is provided in this paper. The flexibility offered by these models
to fit oscillatory signals of many different shapes makes them a very useful
tool to model complex rhythmic patterns. The FMM methodology and its
application to very diverse biological data has been described in previous
papers (Rueda _et al._ , 2019, 2021b, 2021c) and recently revised in Rueda _et
al._ (2021a).

The package allows both single and multicomponent FMM models to be estimated.
In order to provide greater flexibility, equality constraints for shape
parameters have also been implemented. In addition, graphical representations
of the fitted models and the possibility of generating synthetic data are
available. The functionality of the package has been illustrated by simulated
data and also by real examples from different areas of application related to
present-day biological problems.

Possible future extensions of the FMM package include the implementation of
additional restrictions to suit the model to other real signals; the
possibility to include weights that determine how much each observation
influences the parameter estimates; and the improvement of the computational
efficiency. Currently, one of the most important limitations of the stable
version of the package (version $0.1.2$) is that the fitting process,
especially of the restricted model, can be a time-consuming process.
Nevertheless, the execution times of the development version available on the
GitHub repository (https://github.com/alexARC26/FMM) have already been
significantly reduced, in many cases by an order of $10$ times or more.

## Acknowledgments

This research was partially supported by the Spanish Ministerio de Economía y
Competitividad, grant PID2019-106363RB-I00.

## References

  * Anafi _et al._ [2017] Anafi RC, Francey LJ, Hogenesch JB, Kim J (2017).  “CYCLOPS Reveals Human Transcriptional Rhythms in Health and Disease.”  _Proceedings of the National Academy of Sciences_ , 114(20), 5312–5317.  doi:10.1073/pnas.1619320114.
  * Auguie [2017] Auguie B (2017).  _gridExtra : Miscellaneous Functions for Grid Graphics_.  R package version 2.3, URL https://CRAN.R-project.org/package=gridExtra.
  * Bayes de Luna [2007] Bayes de Luna A (2007).  _Basic Electrocardiography: Normal and Abnormal ECG Patterns_.  John Wiley & Sons, Ltd.  doi:10.1002/9780470692622.
  * Boashash [2016] Boashash B (2016).  _Time-Frequency Signal Analysis and Processing: A Comprehensive Reference_.  2nd edition. Academic Press, San Francisco, CA.
  * Carlucci _et al._ [2020] Carlucci M, Kriščiūnas A, Li H, Gibas P, Koncevičius K, Petronis A, Oh G (2020).  “DiscoRhythm: An Easy-To-Use Web Application And R Package for Discovering Rhythmicity.”  _Bioinformatics_ , 36(6), 1952––1954.  doi:10.1093/bioinformatics/btz834.
  * Caro-Martín _et al._ [2018] Caro-Martín, Delgado-García, Gruart, Sánchez-Campusano (2018).  “Spike Sorting Based on Shape, Phase, and Distribution Features, and K-TOPS Clustering with Validity and Error Indices.”  _Scientific Reports_ , 8(1), 1–28.  doi:10.1038/s41598-018-35491-4.
  * Chang _et al._ [2020] Chang W, Cheng J, Allaire J, Xie Y, McPherson J (2020).  _shiny : Web Application Framework for R_.  R package version 1.5.0, URL https://CRAN.R-project.org/package=shiny.
  * Cornelissen [2014] Cornelissen G (2014).  “Cosinor-Based Rhythmometry.”  _Theoretical Biology and Medical Modelling_ , 11, 16.  doi:10.1186/1742-4682-11-16.
  * Corporation and Weston [2020] Corporation M, Weston S (2020).  _doParallel : Foreach Parallel Adaptor for the ’parallel’ Package_.  R package version 1.0.16, URL https://CRAN.R-project.org/package=doParallel.
  * Downs and Mardia [2002] Downs TD, Mardia KV (2002).  “Circular Regression.”  _Biometrika_ , 89(3), 683–698.  doi:10.1093/biomet/89.3.683.
  * Gerstner _et al._ [2014] Gerstner W, Kistler WM, Naud R, Paninski L (2014).  _Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition_.  Cambridge University Press.  doi:10.1017/CBO9781107447615.
  * Gierke _et al._ [2018] Gierke CL, Helget R, Cornelissen-Guillaume G (2018).  _CATkit : Chronomics Analysis Toolkit (CAT): Periodicity Analysis_.  R package version 3.3.3, URL https://CRAN.R-project.org/package=CATkit.
  * Goldberger _et al._ [2000] Goldberger AL, Amaral LA, Glass L, Hausdorff JM, Ivanov PC, Mark RG, Mietus JE, Moody GB, Peng CK, Stanley HE (2000).  “PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals.”  _Circulation_ , 101(23), E215–20.  doi:10.1161/01.cir.101.23.e215.
  * Gouwens _et al._ [2019] Gouwens NW, Sorense SA, Berg J, Lee C, _et al._ (2019).  “Classification of Electrophysiological and Morphological Neuron Types in the Mouse Visual Cortex.”  _Nature Neuroscience_ , 22, 1182––1195.  doi:10.1038/s41593-019-0417-0.
  * Hodgkin and Huxley [1952] Hodgkin AL, Huxley AF (1952).  “A Quantitative Description of Membrane Current and Its Application to Conduction and Excitation in Nerve.”  _The Journal of Physiology_ , 117(4), 500–544.  doi:10.1113/jphysiol.1952.sp004764.
  * Hughes _et al._ [2010] Hughes ME, Hogenesch JB, Kornacker K (2010).  “JTK_CYCLE: An Efficient Nonparametric Algorithm for Detecting Rhythmic Components in Genome Scale Data Sets.”  _Journal of Biological Rhythms_ , 25(5), 372–380.  doi:10.1177/0748730410379711.
  * Kato _et al._ [2008] Kato S, Shimizu K, Shieh GS (2008).  “A Circular - Circular Regression Model.”  _Statistica Sinica_ , 18(2), 633–645.
  * Laguna _et al._ [1997] Laguna P, Mark RG, Goldberg A, Moody GB (1997).  “A Database for Evaluation of Algorithms for Measurement of QT and Other Waveform Intervals in the ECG.”  In _Computers in Cardiology 1997_ , pp. 673–676. IEEE.
  * Larriba _et al._ [2018] Larriba Y, Rueda C, Fernández MA, Peddada SD (2018).  “A Bootstrap Based Measure Robust to the Choice of Normalization Methods for Detecting Rhythmic Features in High Dimensional Data.”  _Frontiers in Genetics_ , 9, 24.  doi:10.3389/fgene.2018.00024.
  * Larriba _et al._ [2020] Larriba Y, Rueda C, Fernández MA, Peddada SD (2020).  “Order Restricted Inference in Chronobiology.”  _Statistics in Medicine_ , 39(3), 265–278.  doi:10.1002/sim.8397.
  * Lynch and Houghton [2015] Lynch EP, Houghton CJ (2015).  “Parameter Estimation of Neuron Models Using In-Vitro and In-Vivo Electrophysiological Data.”  _Frontiers in Neuroinformatics_ , 9, 10.  doi:10.3389/fninf.2015.00010.
  * Mermet _et al._ [2017] Mermet J, Yeung J, Naef F (2017).  “Systems Chronobiology: Global Analysis of Gene Regulation in a 24-Hour Periodic World.”  _Cold Spring Harbor Perspectives in Biology_ , 9(3).  doi:10.1101/cshperspect.a028720.
  * Mosher _et al._ [2020] Mosher CP, Wei Y, Kamiński J, Nandi A, Mamelak AN, Anastassiou CA, Rutishauser U (2020).  “Cellular Classes in the Human Brain Revealed In Vivo by Heartbeat-Related Modulation of the Extracellular Action Potential Waveform.”  _Cell Reports_ , 30(10), 3536–3551.e6.  doi:10.1016/j.celrep.2020.02.027.
  * Moskon [2020] Moskon M (2020).  “CosinorPy: A python Package for Cosinor-Based Rhythmometry.”  _BMC Bioinformatics_ , 21(1), 485.  doi:10.1186/s12859-020-03830-w.
  * Mutak [2018] Mutak A (2018).  _cosinor2 : Extended Tools for Cosinor Analysis of Rhythms_.  R package version 0.2.1, URL https://CRAN.R-project.org/package=cosinor2.
  * Nelder and Mead [1965] Nelder JA, Mead R (1965).  “A Simplex Method for Function Minimization.”  _The Computer Journal_ , 7(4), 308–313.  doi:10.1093/comjnl/7.4.308.
  * Neuwirth [2014] Neuwirth E (2014).  _RColorBrewer : ColorBrewer Palettes_.  R package version 1.1-2, URL https://CRAN.R-project.org/package=RColorBrewer.
  * Parsons _et al._ [2019] Parsons R, Parsons R, Garner N, Oster H, Rawashdeh O (2019).  “CircaCompare: A Method to Estimate and Statistically Support Differences in Mesor, Amplitude and Phase, Between Circadian Rhythms.”  _Bioinformatics_ , 36(4), 1208–1212.  doi:10.1093/bioinformatics/btz730.
  * Rácz _et al._ [2020] Rácz M, Liber C, Németh E, Fiáth R, Rokai J, Harmati I, Ulbert I, Márton G (2020).  “Spike Detection and Sorting with Deep Learning.”  _Journal of Neural Engineering_ , 17(1), 016038.  doi:10.1088/1741-2552/ab4896.
  * R Core Team [2020] R Core Team (2020).  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria.  URL https://www.R-project.org/.
  * Revelle [2018] Revelle W (2018).  _psych : Procedures for Psychological, Psychometric, and Personality Research_.  R package version 1.8.12, URL https://CRAN.R-project.org/package=psych.
  * Rodríguez-Collado and Rueda [2021a] Rodríguez-Collado A, Rueda C (2021a).  “Electrophysiological and Transcriptomic Features Reveal a Circular Taxonomy of Cortical Neurons.”  doi:10.1101/2021.03.24.436849.  Preprint on webpage at https://www.biorxiv.org/content/10.1101/2021.03.24.436849v3.
  * Rodríguez-Collado and Rueda [2021b] Rodríguez-Collado A, Rueda C (2021b).  “A Simple Parametric Representation of the Hodgkin-Huxley Model.”  doi:10.1101/2021.01.11.426189.  Preprint on webpage at https://www.biorxiv.org/content/10.1101/2021.01.11.426189v1.
  * Rueda _et al._ [2021a] Rueda C, Fernandez I, Larriba Y, Rodriguez-Collado A (2021a).  “The FMM Approach to Analyze Biomedical Signals: Theory, Software, Applications and Future.”  Preprint.
  * Rueda _et al._ [2021b] Rueda C, Larriba Y, Lamela A (2021b).  “The Hidden Waves in the ECG Uncovered Revealing a Sound Automated Interpretation Method.”  _Scientific Reports_ , 11, 3724.  doi:10.1038/s41598-021-82520-w.
  * Rueda _et al._ [2019] Rueda C, Larriba Y, Peddada SD (2019).  “Frequency Modulated Möbius Model Accurately Predicts Rhythmic Signals in Biological and Physical Sciences.”  _Scientific Reports_ , 9(1), 18701.  doi:10.1038/s41598-019-54569-1.
  * Rueda _et al._ [2021c] Rueda C, Rodríguez-Collado A, Larriba Y (2021c).  “A Novel Wave Decomposition for Oscillatory Signals.”  _IEEE Transactionns on Signal Processing_ , 69, 960––972.  doi:10.1109/TSP.2021.3051428.
  * Sachs [2014] Sachs M (2014).  _cosinor : Tools for Estimating and Predicting the Cosinor Model_.  R package version 1.1, URL https://CRAN.R-project.org/package=cosinor.
  * Shah [2020] Shah AS (2020).  _card : Cardiovascular and Autonomic Research Design_.  R package version 0.1.0, URL https://CRAN.R-project.org/package=card.
  * Singer and Hughey [2019] Singer JM, Hughey JJ (2019).  “LimoRhyde: A Flexible Approach for Differential Analysis of Rhythmic Transcriptome Data.”  _Journal of Biological Rhythms_ , 34(1), 5–18.  doi:10.1177/0748730418813785.
  * Souza _et al._ [2019] Souza BC, dos Santos VL, Bacelo J, Tort AB (2019).  “Spike Sorting with Gaussian Mixture Models.”  _Scientific Reports_ , 9(1), 1–14.  doi:10.1038/s41598-019-39986-6.
  * Teeter _et al._ [2018] Teeter C, Iyer R, Menon V, Gouwens N, Feng D, Berg J, Szafer A, Cain N, Zeng H, Hawrylycz M, Koch C, Mihalas S (2018).  “Generalized Leaky Integrate-and-Fire Models Classify Multiple Neuron Types.”  _Nature Communications_ , 9(1), 1–15.  doi:10.1038/s41467-017-02717-4.
  * Thaben and Westermark [2014] Thaben PF, Westermark PO (2014).  “Detecting Rhythms in Time Series with RAIN.”  _Journal of Biological Rhythms_ , 29(6), 391–400.  doi:10.1177/0748730414553029.
  * Wickham [2016] Wickham H (2016).  _ggplot2 : Elegant Graphics for Data Analysis_.  Springer-Verlag New York.  ISBN 978-3-319-24277-4.  URL https://ggplot2.tidyverse.org.
  * Zhang _et al._ [2014] Zhang R, Lahens NF, Ballance HI, Hughes ME, Hogenesch JB (2014).  “A Circadian Gene Expression Atlas in Mammals: Implications for Biology and Medicine.”  _Proceedings of the National Academy of Sciences_ , 111(45), 16219–16224.  doi:10.1073/pnas.1408886111.


11institutetext: Universidade Federal da Bahia, Brazil, 11email:
<EMAIL_ADDRESS>22institutetext: Imperial College London, UK, 22email:
<EMAIL_ADDRESS>
# A scaling limit of the parabolic Anderson model with exclusion interaction

Dirk Erhard1 and Martin Hairer2

###### Abstract

We consider the (discrete) parabolic Anderson model $\partial u(t,x)/\partial
t=\Delta u(t,x)+\xi_{t}(x)u(t,x)$, $t\geq 0$, $x\in\mathbb{Z}^{d}$, where the
$\xi$-field is $\mathbb{R}$-valued and plays the role of a dynamic random
environment, and $\Delta$ is the discrete Laplacian. We focus on the case in
which $\xi$ is given by a properly rescaled symmetric simple exclusion process
under which it converges to an Ornstein–Uhlenbeck process. Scaling the
Laplacian diffusively and restricting ourselves to a torus, we show that in
dimension $d=3$ upon considering a suitably renormalised version of the above
equation, the sequence of solutions converges in law.

As a by-product of our main result we obtain precise asymptotics for the
survival probability of a simple random walk that is killed at a scale
dependent rate when meeting an exclusion particle. Our proof relies on the
discrete theory of regularity structures of [ErhardHairerRegularity] and on
novel sharp estimates of joint cumulants of arbitrary large order for the
exclusion process. We think that the latter is of independent interest and may
find applications elsewhere. Keywords: Simple exclusion process, parabolic
Anderson, stochastic PDE MSC2020 classification: 60K35, 60L30, 60L90, 60H15

###### Contents

  1. 1 Introduction

##  1 Introduction

The (discrete) parabolic Anderson model is the partial differential equation

$\partial_{t}u(x,t)=(\Delta u)(x,t)-\xi_{t}(x)u(x,t)\;.{\hbox{\rm(1)}\kern
5.69046pt\kern-5.69046pt}q:PAM{=}e0{1}$