# Some mathematical models for flagellar activation mechanisms

François Alouges, Irene Anello, Antonio Desimone, Aline Lefebvre-Lepot, Jessie Levillain

► To cite this version:

François Alouges, Irene Anello, Antonio Desimone, Aline Lefebvre-Lepot, Jessie Levillain. Some mathematical models for flagellar activation mechanisms. Mathematical Models and Methods in Applied Sciences, 2025, 35 (11), pp.2395-2424. 10.1142/S0218202525500423 . hal-04689234v2

HAL Id: hal-04689234

https://hal.science/hal-04689234v2

Submitted on 16 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Some mathematical models for flagellar
activation mechanisms

François Alouges<fnref n="1" />, Irene Anello², Antonio DeSimone², Aline Lefebvre-Lepot³, and Jessie Levillain⁴

¹Centre Borelli, ENS Paris-Saclay, CNRS, Université Paris-Saclay, 4 avenue des sciences, 91190 Gif-sur-Yvette, France, Institut Universitaire de France (IUF)

²Scuola Internazionale Superiore di Studi Avanzati, via Bonomea 265, I-34136 Trieste, Italy

³CNRS, Fédération de Mathématiques de CentraleSupélec, 9, rue Joliot Curie, 91190 Gif-sur-Yvette

⁴CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, route de Saclay, 91120 Palaiseau, France

April 16, 2025

Abstract

This paper focuses on studying a model for dyneins, cytoskeletal motor proteins responsible for axonemal activity. The model is a coupled system of partial differential equations inspired by [16, 18] and incorporating two rows of molecular motors between microtubules filaments. Existence and uniqueness of a solution is proved, together with the presence of a super-critical Hopf bifurcation. Additionally, numerical simulations are provided to illustrate the theoretical results. A brief study on the generalization to *N*-rows is also included.

Introduction

Molecular motors are proteins that move along tubulin filaments, converting chemical energy (ATP) into mechanical work. Our focus is on dyneins, a specific type of molecular motor found in the axoneme, the cytoskeletal structure of cilia. The axoneme consists of nine pairs of tubulin filaments<fnref n="1" />, known as microtubules, and functions as an active structure due to the dyneins positioned between

<footnote n="1">¹On some biological systems there is also a central pair of microtubules [32]. We don't consider it in the present study.</footnote>
adjacent microtubule pairs. These molecular motors attach to one microtubule pair and walk along the neighboring pair, generating local sliding between them. When the microtubules are considered elastic, as in cilia, this sliding induces bending. Along the microtubules, the molecular structure is periodic with periodicity $\ell \ll L$ where L is the total length of the filament.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** A microscopic slice of the axoneme (on the right), zoomed in from the whole filament (on the left). The periodic structure is highlighted by an alternation of green and blue sections, modeling preferred binding sites for molecular motors.</figcaption></figure>

We analyze a microscopic slice of thickness $\ell$ of the microtubules (See Figure 1), focusing exclusively on the interaction between molecular motors and the filaments. At this scale, the axoneme is considered *isolated*, i.e. there is no fluid to be taken into account, and the filaments are treated as rigid. The ultimate objective of this article is to model the entire structure of such a slice of the isolated axoneme, composed of nine microtubule pairs.

As a starting point, we adopt the *one-row model* described in [16, 15, 18, 17, 12], which characterizes the collective motion of a single row of molecular motors attached to a moving filament against a fixed counterpart. A key finding of this model is the emergence of mesoscopic order from the independent activity of individual molecular motors. When an elastic force is present, an analysis of the linearized system coming from the model shows that the system might possess a supercritical Hopf bifurcation instability [12]. Numerical simulations [12, 16, 15, 17, 18] confirm this kind of behavior.

Building on the one-row model, and in order to model the circular structure of the slice of the axoneme, we introduce a second row of molecular motors, as shown in Figure 2(b), forming a *two-row model* that symmetrizes the previous framework. This simplified structure was already suggested e.g. in [24] Fig. 1. A & B, or in [33] Fig. 2., A & B, though not studied as an isolated system. We further generalize this to an *N-row model*, where N represents the number of motor rows – typically 8 – and, consequently, the number of microtubule pairs. Both models are analytically described in the next section.

The novelty of our work lies in four key contributions: (i) We generalize the one-row model by first introducing an alternative description for the isolated
<figure><img src="image_3.png" /><figcaption>**Figure 2:** (a) Cross section of the axoneme. (b) Cross section of the axoneme where we group the filaments 1-4 (blue) and 5-9 (red) together. (c) Unfolding of the axoneme in (b) with two opposite microtubule pairs.</figcaption></figure>

axoneme: the two-row model. We then use this model as a first step towards a more realistic *N*-row model. (ii) Then we prove the existence and uniqueness of solutions for the two-row model which is composed of coupled non-linear PDEs. This result fills a gap in the existing theoretical literature and also extends to the one-row or *N*-row models. (iii) We establish rigorously the presence of a supercritical Hopf bifurcation in the system for general parameters. The non-linearity of the equations is treated through the Center Manifold technique. We further compute the amplitude of the periodic solution after the instability has occurred. (iv) Eventually, we introduce a specific upwind scheme to perform numerical simulations on the *N*-row model for *N* ≥ 1.

The paper is organized as follows: the first section presents the derivation of the one-row, two-row and *N*-row models and the motivation underlying the increasing biology complexity. In the second section, we state and prove the existence and uniqueness of the solutions (Theorem 2.1) and the existence of the Hopf bifurcation (Theorem 2.2). In the final section, we describe the numerical scheme employed in the paper, present the numerical simulations, and explore the questions arising from the *N*-row structure.

# 1 Motivation and Modeling

## 1.1 The one-row model

In all the presented models, the molecular motors are described by a generalized version of the stochastic two-state model proposed in [16]. The key difference between the one-, two-, and *N*-row models lies in the number of motor rows, which determines the number of filaments representing the slice of the isolated axoneme.

In the two-state model, each motor lies between two filaments. It is attached to one of them, and has two identical heads (head *A* and head *B*), each of which can be either bound or unbound to the opposite filament. The motor can exist
in two distinct chemical states: state A, where head A is bound and head B is unbound, and state B, where the opposite happens. Importantly, the model assumes that both heads cannot be simultaneously unbound or bound. In the one-row model, motors are attached to the moving top filament and exert a force on the stationary bottom filament. Between the pair of filaments, passive elastic and viscous elements resist the motion. They are modeled by the positive constants k (elastic coefficient) and η (viscous coefficient), respectively.

Along its length the filaments present a microscopic periodic structure of period $ℓ$ that drives the chemical behavior of the heads. Therefore, we consider a slice of thickness $ℓ$ of the isolated axoneme. Then, each state $i = A, B$ experiences a ($ℓ$-periodic) potential energy $W_i(\xi)$ at position $\xi \in [0, ℓ]$ on the fixed bottom filament. The $ℓ$-periodic transition rates $\omega_A(\xi)$ and $\omega_B(\xi)$ represent the probability per unit time for a motor to switch from state A to state B and from state B to state A respectively. The transition rates also depend on the ATP concentration $\Omega$, i.e. $\omega_i = \omega_i(\xi; \Omega)$, as ATP provides the energy for state transitions. A detailed explanation of this dependence is given in [15]. Note that more complex approaches exist in the literature, for which transition rates are also force-dependent. We refer the reader to [34, 35, 33] for detailed explanations on such models.

In the limit of an infinite number of motors on $[0, ℓ]$, we introduce the densities $\tilde{P}_i(\xi) \ge 0$ (for $i = A, B$), which represent the probability of finding a motor at position $\xi$ in state $i$. These densities satisfy $\tilde{P}_A(\xi) + \tilde{P}_B(\xi) = 1/ℓ$, reflecting the fact that a motor is always in either state A or state B. The position of the top filament is denoted by $x(t)$, with velocity $v(t) = d/dt x(t)$.

The one-row model, as described in [16], expresses the motion of the system through a transport equation for $P = \tilde{P}_A$:

$$ \left\{ \begin{array}{l} \partial_t P(\xi, t) + v(t) \partial_\xi P(\xi, t) = -(\omega_A + \omega_B) P(\xi, t) + \omega_B/\ell \\ v(t) = \frac{1}{\eta} \left( \int_0^\ell P(\xi, t) \partial_\xi \Delta W(\xi) d\xi - kx(t) \right), \end{array} \right. \quad (1) $$

where $\Delta W(\xi) = W_B(\xi) - W_A(\xi)$. The second equation in system (1) represents the force balance on the moving filament:

$$ f_{\text{ext}} - \eta v - kx + f_{\text{mot}} = 0, $$

with zero external force ($f_{\text{ext}} = 0$) since we don't consider fluid-structure interaction here and where the force due to the motors is given by

$$ f_{\text{mot}}(t) = \int_{0}^{t} P(\xi, t) \partial_{\xi} \Delta W(\xi) d\xi. $$

The system (1) admits a stationary, non-moving solution:

$$ P_{\text{eq}}(\xi; \Omega) = \frac{\omega_B}{\ell(\omega_A + \omega_B)}, \quad v_{\text{eq}} = 0, \quad x_{\text{eq}} = \frac{1}{k} \int_{0}^{l} P_{\text{eq}}(\xi; \Omega) \partial_{\xi} \Delta W(\xi) d\xi. $$
Depending on the specific forms of the transition rates $\omega_i$ and potentials $W_i$, a threshold ATP concentration $\Omega = \Omega_0$ may exist, above which this stationary solution becomes unstable, leading to oscillatory behavior in the system (see [12, 16, 15, 17, 18]).

## 1.2 The two-row model

We now derive in detail the two-row model, which is depicted in Figure 2(b). The two-row model is a simplified version of the circular structure of the axoneme ([24, 7]), forming a cylinder with only two pairs of microtubules linked by two rows of motors, as shown in Figure 2(a). This model is a first step towards the N-row model presented in section 1.3, where N corresponds to the number of motor rows.

Moreover, the two-row model symmetrizes the one-row model. In the one-row model [16], even without any external force, the moving filaments experiences a non-zero equilibrium displacement $x(t) = x_{eq}$, which disrupts axonemal symmetry and causes bending at rest. Jülicher and Camalet mentioned this issue in [6], solving it by choosing symmetrical potentials, as in our case, and symmetrical transition rates, namely $|\omega_i(\xi + l/2)| = |\omega_i(-\xi + l/2)|$, with $\xi \in [0, \ell]$ and $i = A, B$. Instead, in this work, we consider general non-symmetrical transition rates and introduce a second row of molecular motors, producing a zero equilibrium $x_{eq} = 0$ at rest.

Unfolding the two-row structure, we may consider three pairs of filaments, numbered from 0 to 2 where the pairs 0 and 2 are identified and fixed, while the central pair 1, may move (See Figure 2(c)). We call $x(t)$ its displacement and $v(t) = \frac{d}{dt}x(t)$ its velocity. We define $P_1(\xi, t) = \tilde{P}_{1,A}(\xi, t)$ and $P_2(\xi, t) = \tilde{P}_{2,A}(\xi, t)$ as the probabilities for the motors to be in state A at position $\xi$ and time $t$, when attached to the pair 1 (bottom row) or to the pair 2 (top row), respectively. We assume that the sum of the transition rates is uniform, (see [6, 12]),

$$ \omega_A(\xi; \Omega) + \omega_B(\xi; \Omega) = a_0(\Omega), \quad (2) $$

where $a_0(\Omega) > 0$, and obtain two transport equations for $\xi \in [0, \ell]$ and $t > 0$:

$$ \left\{ \begin{array}{l} \frac{\partial P_1}{\partial t}(\xi, t) + v(t) \frac{\partial P_1}{\partial \xi}(\xi, t) = -a_0(\Omega) P_1(\xi, t) + \frac{\omega_B(\xi; \Omega)}{\ell}, \\ \\ \frac{\partial P_2}{\partial t}(\xi, t) = -a_0(\Omega) P_2(\xi, t) + \frac{\omega_B(\xi - x(t); \Omega)}{\ell}. \end{array} \right. \quad (3) $$

Additionally, the force balance equation on the central filament reads:

$$ -2\eta v - 2kx + f_{\text{mot}} = 0, \quad (4) $$

where as before, no external force is considered, $kx$ and $\eta v$ represent the elastic and viscous resistances, respectively and where $f_{\text{mot}}$ is the active force exerted by the motors, which is now given by

$$ f_{\text{mot}}(t) = \int_{0}^{\ell} (P_{1}(\xi, t) \partial_{\xi} \Delta W(\xi) - P_{2}(\xi, t) \partial_{\xi} \Delta W(\xi - x(t))) d\xi. \quad (5) $$
Let $P(\xi, t) = P_1(\xi, t)$ and $Q(\xi, t) = P_2(\xi + x(t), t)$. Then, (3) becomes

$$ \left\{ \begin{aligned} \frac{\partial P}{\partial t}(\xi, t) + v(t) \frac{\partial P}{\partial \xi}(\xi, t) &= -a_0(\Omega) P(\xi, t) + \frac{\omega_B(\xi)}{\ell}, \\ \frac{\partial Q}{\partial t}(\xi, t) - v(t) \frac{\partial Q}{\partial \xi}(\xi, t) &= -a_0(\Omega) Q(\xi, t) + \frac{\omega_B(\xi)}{\ell}. \end{aligned} \right. \qquad (6) $$

The periodicity of $P_1$ and $P_2$, enables us to rewrite the motor force as

$$ f_{mot}(t) = \int_{0}^{\ell} (P(\xi, t) - Q(\xi, t)) \partial_{\xi} \Delta W(\xi) d\xi, $$

and obtain, from (4),

$$ v(t) = \dot{x}(t) = \frac{1}{2\eta} \left( \int_{0}^{\ell} (P(\xi, t) - Q(\xi, t)) \partial_{\xi} \Delta W(\xi) d\xi - 2kx(t) \right). \qquad (7) $$

Finally, we obtain the system of equations

$$ \left\{ \begin{aligned} \frac{\partial P}{\partial t}(\xi, t) + v(t) \frac{\partial P}{\partial \xi}(\xi, t) &= -(\omega_A(\xi) + \omega_B(\xi)) P(\xi, t) + \frac{\omega_B(\xi)}{\ell}, \\ \frac{\partial Q}{\partial t}(\xi, t) - v(t) \frac{\partial Q}{\partial \xi}(\xi, t) &= -(\omega_A(\xi) + \omega_B(\xi)) Q(\xi, t) + \frac{\omega_B(\xi)}{\ell}, \\ v(t) &= \frac{1}{2\eta} \left( \int_{0}^{\ell} (P(\xi, t) - Q(\xi, t)) \partial_{\xi} \Delta W(\xi) d\xi - 2kx(t) \right), \\ P(0, t) = P(\ell, t), \quad Q(0, t) = Q(\ell, t). \end{aligned} \right. \qquad (8) $$

The equilibrium state for this system is,

$$ \left\{ \begin{aligned} v_{eq} &= 0, \\ x_{eq} &= 0, \\ P_{eq} &= Q_{eq} = \frac{\omega_B}{\ell(\omega_A + \omega_B)}. \end{aligned} \right. $$

Notice that while the probability equilibrium is the same as in the one-row model, the equilibrium position is now zero. This system will then be studied both theoretically and numerically in the next sections, assuming suitable initial conditions for $P$ and $Q$.

## 1.3 The N-row model

Starting from two rows of molecular motors, we extend the model from section 1.2 to $N \ge 2$ rows of molecular motors.

As illustrated in figure 3 The system is then composed of $N + 1$ microtubule
<figure><img src="image_4.png" /><figcaption>Figure 3: Simplified axoneme with 9 microtubule doublets and no central pair [31].</figcaption></figure>

doublets on the outside, arranged in a circle, and no central pair. All filaments have the same polarity, meaning that between two microtubule pairs, motors move towards the base. As before, filament pairs 0 and N are considered to be the same pair in order to take into account the cylindrical periodic structure. They have no shift between each other, as shown in the split view Figure 4. Moreover, all microtubule pairs are assumed inextensible, rigid and at a constant distance from each other. They are held together by elastic and viscous elements that do not permit infinite sliding between them. We look at molecular motors in a periodicity cell of an N-axoneme, of length $\ell$.

<figure><img src="image_5.png" /><figcaption>Figure 4: Motors and tubule shift in the N-row model (here the horizontal axis x shows absolute tubule displacement instead of motor position $\xi$ in a periodicity cell).</figcaption></figure>

For $i \in \{0, \dots, N\}$, we denote by $x_i$ the horizontal shift of the i-th tubule, as shown in Figure 4. We measure displacement with respect to the 0-th filament: $x_0 = x_N = 0$. For $1 \le i \le N$, we introduce $\Delta_i = x_i - x_{i-1}$, the relative displacement between filaments i and i - 1. It follows that $\sum_{i=1}^N \Delta_i(t) = 0$. The shifting speed of the i-th filament is defined by $v_i = \dot{x}_i - \dot{x}_{i-1} = \dot{\Delta}_i$, and we have

$$ \sum_{i=1}^{N} v_i = \sum_{i=1}^{N} \dot{\Delta}_i(t) = 0. \qquad (9) $$
Once again, the variable $\xi \in [0, \ell]$ represents the local variable along one pair, as in the two-row model. We denote by $P_i(\xi, t) = \tilde{P}_{i,A}(\xi, t)$ the density of molecular motors at position $\xi$ and time $t$ attached to the $i$-th filament and who are in state $A$, walking on the $(i-1)$-th filament. Following again section 1.2, for $j \in \{1, \dots, N\}, t > 0$ and $\xi \in [0, \ell]$, we define $Q_j$ as $Q_j(\xi, t) = P_j(\xi + \Delta_j(t), t)$ and we finally obtain the system at any $t > 0$ and $\xi \in [0, \ell]$

$$ \left\{ \begin{array}{ll} \partial_t Q_j + v_j \partial_\xi Q_j = -(\omega_A + \omega_B)Q_j + \omega_B/\ell, & \text{for } j \in \{1, \dots, N\}, \\ \\ \eta(v_i - v_{i+1}) = \int_0^\ell (Q_i - Q_{i+1})\partial_\xi \Delta W - k(\Delta_i - \Delta_{i+1}), & \text{for } i \in \{1, \dots, N-1\}, \\ \\ \displaystyle\sum_{i=1}^N \dot{\Delta}_i(t) = 0 \text{ and } \dot{x}_N = 0, & \end{array} \right. \qquad (10) $$

where the first $N$ equations are transport equations for the motors probability density and the $N-1$ following equations are the force balances on the $i$-th filament with the motors and springs around it, for $i = 1 \dots N-1$.

Like in previous models, the equilibrium state is, for all $1 \le i \le N-1$ and $1 \le j \le N$

$$ \left\{ \begin{array}{l} v_{eq,i} = 0, \\ x_{eq,i} = 0, \\ Q_{eq,j} = \frac{\omega_B}{\ell(\omega_A + \omega_B)}. \end{array} \right. $$

Notice that again, the probability equilibrium is the same as in the previous models and, as in the two-row setting, the equilibrium position is zero.

## 2 Theoretical results

In this section we state and prove two theorems, giving the existence and uniqueness of solutions and existence of a Hopf bifurcation for two-row model (8). We then give some theoretical consideration on the general $N$-row model (10).

Let us define $C_{\#}^1([0, \ell])$ as the set of restrictions to $[0, \ell]$ of $C^1$ and $\ell$-periodic functions over $\mathbb{R}$, then the first theorem reads as follows:

**Theorem 2.1 (Existence and uniqueness).** Let us fix $\ell > 0$. Assume $\omega_A(\xi; \Omega)$ and $\omega_B(\xi; \Omega)$ as in (2). Moreover, assume $\omega_B$ and $\Delta W$ to be at least $C_{\#}^1([0, \ell])$. If the initial data $P(\xi, 0)$ and $Q(\xi, 0)$ are $C_{\#}^1([0, \ell])$ and $x(0) = 0$, then the system of equations (8), admits a unique solution $P$, $Q \in C^1([0, \ell] \times \mathbb{R}_+)$, with $\xi \mapsto P(\xi, \cdot)$ and $\xi \mapsto Q(\xi, \cdot)$ in $C_{\#}^1([0, \ell])$, and $x \in C^1(\mathbb{R}_+)$.

For the second result we chose $\Delta W(\xi) = U \cos(\frac{2\pi\xi}{\ell})$, as in [6]. Moreover, since the two motor heads are identical, then $\omega_B(\xi) = \omega_A(\xi + \ell/2)$. Taking the
transition rates' periodicity into account, we use the following Fourier expansion

$$ \omega_B(\xi; \Omega) = \frac{a_0(\Omega)}{2} + \sum_{n=2k+1, k \ge 0} \left( a_n(\Omega) \cos \frac{2n\pi\xi}{\ell} + b_n(\Omega) \sin \frac{2n\pi\xi}{\ell} \right), \quad (11) $$

with $a_n(\Omega)$ and $b_n(\Omega)$ real coefficients. The even coefficient disappear due to the fact that $\omega_B(\xi) = \omega_A(\xi + \ell/2)$ and $\omega_A(\xi) + \omega_B(\xi) = a_0$.

**Theorem 2.2 (Hopf bifurcation).** In addition to the hypothesis above, assume that $\Omega \ge 0$, $a_0(\Omega) > 0$, $a_1(\Omega) \ne 0$ and that $a_i(\Omega)$ with $i = 0, 1$ are at least $C^1(\mathbb{R}_+)$. Furthermore, let us define

$$ \tau(\Omega) := -\frac{1}{4} \left( 2a_0(\Omega) + \frac{\zeta\ell}{\pi} + 2\lambda \frac{a_1(\Omega)}{a_0(\Omega)\ell} \right), \quad (12) $$

with $\zeta = 2\pi k/\eta\ell$ and $\lambda = 2\pi^2 U/\eta\ell$. Suppose there exists $\Omega_0 \in \mathbb{R}_+$ such that $\tau(\Omega_0) = 0$ and $\tau'(\Omega_0) > 0$, then the solutions $P(\xi, t; \Omega)$, $Q(\xi, t; \Omega)$ and $x(t; \Omega)$ of the system (8) show a super-critical Hopf bifurcation in time near the bifurcation value $\Omega_0$ and near the fixed point

$$ (P_{eq}(\xi; \Omega), Q_{eq}(\xi; \Omega), x_{eq}(\Omega)) = \left( \frac{\omega_B(\xi; \Omega)}{a_0(\Omega)\ell}, \frac{\omega_B(\xi; \Omega)}{a_0(\Omega)\ell}, 0 \right). $$

In particular, the fixed point is asymptotically stable for $\Omega < \Omega_0$, and unstable for $\Omega > \Omega_0$. Moreover there exists an asymptotically stable periodic orbit of frequency $\omega_0$ given by and frequency

$$ \omega_0 = \sqrt{\frac{\zeta \ell a_0(\Omega_0)}{2\pi}} \quad (13) $$

for $\Omega > \Omega_0$. In that case, the displacement $x(t)$ has an amplitude given by

$$ \rho = \sqrt{\frac{(\Omega - \Omega_0)\tau'(\Omega_0)}{\tilde{\tau}}} + o\left(\sqrt{\Omega - \Omega_0}\right), \quad (14) $$

$$ \text{where } \tilde{\tau} = -\frac{3\pi\zeta}{4\ell} \left( \frac{\pi a_0(\Omega_0) + \ell\zeta}{\pi a_0(\Omega_0) + 2\ell\zeta} \right). $$

The methods introduced to prove the theorems can be also used to demonstrate similar results for the one-row model, and might prove useful to consider the general N-row structure. These two results rigorously prove the existence of a solution for system (8) which bifurcates in time for a given value of the ATP concentration. The relative sliding $x(t; \Omega)$ between microtubules shows oscillatory motion when $\Omega > \Omega_0$, depending on the coefficients $a_0, a_1$ in (11). This happens as a result of the binding and unbinding coordinated activity of molecular motors, whose dynamics is represented by $P(\xi, t; \Omega)$ and $Q(\xi, t; \Omega)$. The behavior of the bifurcation will be investigated at the end of the following section in terms of the way the coefficients $a_0, a_1$ depend on the ATP concentration $\Omega$. Notice that the other Fourier coefficients $a_n$, with $n > 1$, and $b_n$, with $n \ge 1$, don't play any role in the onset of bifurcation but influence the probabilities' amplitude of oscillation.
## 2.1 Proof of Theorem 2.1 (existence and uniqueness)

The proof follows a classical scheme in two steps: we first show local existence of the solution in time, and then extend it to $\mathbb{R}_+$.

**Step 1: local existence.** Let $T > 0$ be a time to be chosen afterwards and define the map $\psi: C^0([0, T]) \to C^0([0, T])$ as $\psi[x(\cdot)](t) = \int_0^t F[x(\cdot)](s) ds$, where

$$F[x(\cdot)](s) = \frac{1}{2\eta} \int_0^\ell (P_x(\xi, s) - Q_x(\xi, s)) \partial_\xi \Delta W(\xi) d\xi - \frac{k}{\eta} x(s), \quad (15)$$

and the functions $P_x(\xi, t)$ and $Q_x(\xi, t)$ are defined as the solutions of (6) and $v(t) = \dot{x}(t)$. Notice that $P_x$ and $Q_x$ are explicitly given in their integral form by

$$P_x(\xi, t) = e^{-a_0 t} P(\xi - x(t) + x(0), 0) + \frac{e^{-a_0 t}}{\ell} \int_0^t e^{a_0 s} \omega_B(\xi + x(s) - x(t)) ds, \quad (16)$$

and

$$Q_x(\xi, t) = e^{-a_0 t} Q(\xi + x(t) - x(0), 0) + \frac{e^{-a_0 t}}{\ell} \int_0^t e^{a_0 s} \omega_B(\xi - x(s) + x(t)) ds, \quad (17)$$

which remain well-defined when $x$ is only in $C^0([0, T])$.

Notice that $(\tilde{x}, P_{\tilde{x}}, Q_{\tilde{x}})$ is a solution of the system (6, 7) if and only if $\tilde{x}(t)$ is a fixed point of $\psi$, i.e. $\psi[\tilde{x}(\cdot)] = \tilde{x}(\cdot)$. In order to proceed, we prove that $\psi: C^0([0, T]) \to C^0([0, T])$ is a strict contraction for $T$ sufficiently small.

Let us take two functions $x_1$ and $x_2$ in $C^0([0, T])$, with $x_1(0) = x_2(0) = 0$ and with the corresponding $P_{x_1}$, $Q_{x_1}$ and $P_{x_2}$, $Q_{x_2}$ defined by the integral formulations (16, 17). The initial conditions are identical $P_{x_1}(\xi, 0) = P_{x_2}(\xi, 0) = P(\xi)$ and $Q_{x_1}(\xi, 0) = Q_{x_2}(\xi, 0) = Q(\xi)$. We want to estimate the quantity

$$\psi[x_1(\cdot)](t) - \psi[x_2(\cdot)](t) = \int_0^t F[x_1(\cdot)](s) - F[x_2(\cdot)](s) ds. \quad (18)$$

Using (15), we get

$$F[x_1(\cdot)](s) - F[x_2(\cdot)](s) = \frac{1}{2\eta} \left( \begin{aligned} & \int_0^\ell (P_{x_1}(\xi, s) - P_{x_2}(\xi, s)) \\ & - (Q_{x_1}(\xi, s) - Q_{x_2}(\xi, s)) \end{aligned} \right) \partial_\xi \Delta W(\xi) d\xi - \frac{k}{\eta} (x_1(s) - x_2(s)). \quad (19)$$

Now, we have from (16)

$$P_{x_1}(\xi, t) - P_{x_2}(\xi, t) = e^{-a_0 t} (P_{x_1}(\xi - x_1(t), 0) - P_{x_2}(\xi - x_2(t), 0)) \\ + \frac{e^{-a_0 t}}{\ell} \int_0^t e^{a_0 s} (\omega_B(\xi + x_1(s) - x_1(t)) - \omega_B(\xi + x_2(s) - x_2(t))) ds. \quad (20)$$
Thus, we have the estimate

$$
\begin{align*}
\| P_{x_1}(\xi, t) - P_{x_2}(\xi, t) \|_{L_{\xi}^{\infty} L_{t}^{\infty}}
&\le \| \partial_{\xi} P \|_{L_{\xi}^{\infty}} \| x_1(t) - x_2(t) \|_{L_{t}^{\infty}} \\
&\quad + \frac{1}{\ell} \left\| \int_{0}^{t} e^{a_0(s-t)} \| \partial_{\xi} \omega_B \|_{L_{\xi}^{\infty}} |x_1(s) - x_1(t) - x_2(s) + x_2(t)| ds \right\|_{L_{t}^{\infty}} \\
&\le \left( \| \partial_{\xi} P \|_{L_{\xi}^{\infty}} + \frac{2}{a_0 \ell} \| \partial_{\xi} \omega_B \|_{L_{\xi}^{\infty}} \right) \| x_1(t) - x_2(t) \|_{L_{t}^{\infty}}.
\end{align*}
$$

The same goes for the difference $Q_{x_1}(\xi, t) - Q_{x_2}(\xi, t)$, and we deduce

$$
\| P_{x_1}(\xi, t) - P_{x_2}(\xi, t) - Q_{x_1}(\xi, t) + Q_{x_2}(\xi, t) \|_{L_{\xi}^{\infty} L_{t}^{\infty}} \le C_1 \| x_1(t) - x_2(t) \|_{L_{t}^{\infty}}. \quad (21)
$$

with $C_1 = 4 \max \{ \|\partial_\xi Q\|_{L_\xi^\infty}, \|\partial_\xi P\|_{L_\xi^\infty}, \frac{2}{a_0\ell} \|\partial_\xi \omega_B\|_{L_\xi^\infty} \}$.
We now deduce from (19)

$$
\|F[x_1(\cdot)] - F[x_2(\cdot)]\|_{L_t^\infty} \le \left( \frac{C_1 \ell}{2\eta} \|\partial_\xi \Delta W\|_{L_\xi^\infty} + \frac{k}{\eta} \right) \|x_1(t) - x_2(t)\|_{L_t^\infty}, \quad (22)
$$

and obtain

$$
\| \psi[x_1(\cdot)] - \psi[x_2(\cdot)] \|_{L_t^\infty} \le C_2 T \| x_1 - x_2 \|_{L_t^\infty}, \quad (23)
$$

with $C_2 = \frac{C_1 \ell}{2 \eta} \| \partial_\xi \Delta W \|_{L_\xi^\infty} + k/\eta$.

Taking $T = \frac{1}{2C_2}$, this proves that $\psi$ is a contraction, as claimed. Then, there exists a unique fixed point $\tilde{x}(\cdot) \in C^0([0, T])$ which satisfies (7).

**Step 2: global solutions.** The previous construction can be extended as long as $P_x$ and $Q_x$ remain bounded in $C^1([0, \ell])$ with respect to $\xi$ and in $L^\infty(\mathbb{R}_+)$ with respect to $t$, as shown by the formulas for $C_1$ and $C_2$. But from (16, 17) we have

$$
\| P_x(\cdot, t) \|_{C_{\xi}^{1}} \leq e^{-a_0 t} \| P \|_{C_{\xi}^{1}} + (1 - e^{-a_0 t}) \frac{1}{a_0 \ell} \| \omega_B \|_{C_{\xi}^{1}},
$$

which shows that

$$
\| P_x \|_{C_{\xi}^{1}, L_{t}^{\infty}} \leq \max \left( \| P \|_{C_{\xi}^{1}}, \frac{1}{a_0 \ell} \| \omega_B \|_{C_{\xi}^{1}} \right),
$$

and the same goes for $Q_x$. We thus obtain that there is a unique global solution
$(x, P_x, Q_x)$ to (6, 7) for all time $t \ge 0$.

In fact, $x \in C^1(\mathbb{R}_+)$ and $P_x, Q_x \in C^1([0, \ell] \times \mathbb{R}_+)$, with $\xi \mapsto P(\xi, \cdot)$ and $\xi \mapsto Q(\xi, \cdot)$ in $C^1([0, \ell])$, as shown by the following bootstrap argument: From equations (16, 17), we know that $P_x$ and $Q_x$ are continuous in time, which enables us to deduce that $F[x(\cdot)] \in C^0(\mathbb{R}_+)$. Therefore, $\psi[x(\cdot)] \in C^1(\mathbb{R}_+)$. But, since $x = \psi[x(\cdot)]$, we deduce that $x \in C^1(\mathbb{R}_+)$. Re-using equations (16, 17), we may then deduce that $P_x$ and $Q_x$ are in fact $C^1$ in time (and they actually have the minimal regularity of the initial conditions $P, Q$ and of $\omega_B$).
## 2.2 Proof of Theorem 2.2 (Hopf bifurcation)

The functions $\xi \mapsto P(\xi, t)$ and $\xi \mapsto Q(\xi, t)$ are periodic with period $l$, and they can be then expanded in Fourier series

$$ \left\{ \begin{aligned} P(\xi, t) &= \frac{p_0(t)}{2} + \sum_{n>0} \left( p_n^c(t) \cos \frac{2n\pi\xi}{l} + p_n^s(t) \sin \frac{2n\pi\xi}{l} \right), \\ Q(\xi, t) &= \frac{q_0(t)}{2} + \sum_{n>0} \left( q_n^c(t) \cos \frac{2n\pi\xi}{l} + q_n^s(t) \sin \frac{2n\pi\xi}{l} \right). \end{aligned} \right. \quad (24) $$

We insert this expansion and the one for the transition rates (11) into the system (8). By matching same order terms, we get an infinite number of ordinary differential equations for the coefficients of $P$ and $Q$. Namely, for $n = 0$ we get two decoupled equations for $p_0$ and $q_0$:

$$ \dot{p}_0(t) = -a_0(\Omega)p_0(t) + a_0(\Omega)/2\ell, \quad \dot{q}_0(t) = -a_0(\Omega)q_0(t) + a_0(\Omega)/2\ell. \quad (25) $$

For $n \neq 0$ we obtain:

$$ \left\{ \begin{aligned} \dot{p}_n^c(t) + \frac{2\pi n}{\ell}v(t)p_n^s(t) &= -a_0(\Omega)p_n^c(t) + a_n(\Omega)/\ell, \\ \dot{p}_n^s(t) - \frac{2\pi n}{\ell}v(t)p_n^c(t) &= -a_0(\Omega)p_n^s(t) + b_n(\Omega)/\ell, \\ \dot{q}_n^c(t) - \frac{2\pi n}{\ell}v(t)q_n^s(t) &= -a_0(\Omega)q_n^c(t) + a_n(\Omega)/\ell, \\ \dot{q}_n^s(t) + \frac{2\pi n}{\ell}v(t)q_n^s(t) &= -a_0(\Omega)q_n^s(t) + b_n(\Omega)/\ell, \end{aligned} \right. \quad (26) $$

together with the force balance equation

$$ 2\eta\dot{x}(t) + 2kx(t) + \pi U(p_1^s - q_1^s) = 0. \quad (27) $$

Note that the coupling between the probabilities evolution and the force equilibrium equation takes place only for the first order coefficients. We are going to prove the existence of a Hopf bifurcation by treating order $n = 0$ first, then $n = 1$ and lastly $n > 1$. Combining all three results will complete the proof.

**Step 1.** Zeroth order coefficients. It is clear that (25) gives

$$ p_0(t) = e^{-a_0 t} \left( p_0(0) - \frac{1}{2\ell} \right) + \frac{1}{2\ell}, \quad q_0(t) = e^{-a_0 t} \left( q_0(0) - \frac{1}{2\ell} \right) + \frac{1}{2\ell}, $$

and both converge to $\frac{1}{2\ell}$ exponentially fast.

**Step 2.** First order coefficients. We now prove the onset of oscillatory patterns for the first order coefficients $p_1^{c,s}(t)$, $q_1^{c,s}(t)$, and for the position $x(t)$.
From (26) and (27) we obtain a first-order five-dimensional ODE:

$$ egin{cases} \dot{p}_1^c = -a_0(\Omega)p_1^c + \left(\zeta x + \frac{\lambda}{2}(p_1^s - q_1^s)\right)p_1^s + \frac{a_1(\Omega)}{\ell}, \\ \dot{p}_1^s = -a_0(\Omega)p_1^s - \left(\zeta x + \frac{\lambda}{2}(p_1^s - q_1^s)\right)p_1^c + \frac{b_1(\Omega)}{\ell}, \\ \dot{q}_1^c = -a_0(\Omega)q_1^c - \left(\zeta x + \frac{\lambda}{2}(p_1^s - q_1^s)\right)q_1^s + \frac{a_1(\Omega)}{\ell}, \\ \dot{q}_1^s = -a_0(\Omega)q_1^s + \left(\zeta x + \frac{\lambda}{2}(p_1^s - q_1^s)\right)q_1^c + \frac{b_1(\Omega)}{\ell}, \\ \dot{x} = -\frac{\ell}{2\pi}\left(\zeta x + \frac{\lambda}{2}(p_1^s - q_1^s)\right), \end{cases} \quad (28) $$

where $\zeta = 2\pi k/\eta\ell$ and $\lambda = 2\pi^2 U/\eta\ell$.

We linearize the system around its equilibrium point $p_{eq}(\Omega) = (p_{1,eq}^c(\Omega), p_{1,eq}^s(\Omega), q_{1,eq}^c(\Omega), q_{1,eq}^s(\Omega), x_{eq}(\Omega))$, where

$$ p_{1,eq}^c = q_{1,eq}^c = \frac{a_1(\Omega)}{a_0(\Omega)\ell}, \quad p_{1,eq}^s = q_{1,eq}^s = \frac{b_1(\Omega)}{a_0(\Omega)\ell}, \quad x_{eq} = 0. $$

We observe that the Jacobian matrix has five eigenvalues: three of them are equal to $-a_0(\Omega) < 0$ while the other two are of the form $\tau(\Omega) \pm \sqrt{-\frac{\zeta\ell a_0(\Omega)}{2\pi} + \tau^2(\Omega)}$, where $\tau(\Omega)$ is given by (12). Since, by hypothesis, there exists a real and positive $\Omega = \Omega_0$ such that $\tau(\Omega_0) = 0$, then $-\zeta\ell a_0(\Omega_0) + 2\pi\tau^2(\Omega_0) < 0$, and we deduce that the pair of complex eigenvalues can be written as $\tau(\Omega) \pm i\omega(\Omega)$ in a neighborhood of $\Omega_0$, where

$$ \omega(\Omega) := \sqrt{\frac{\zeta\ell a_0(\Omega)}{2\pi} - \tau^2(\Omega)}, \quad (29) $$

and they cross the imaginary axis at $\Omega = \Omega_0$.

In the following proposition, we are going to study the non-linear behavior of the vector field solution to (28) by exploiting the center manifold theorem; the orbit structure near the fixed point and $\Omega_0$ is determined by the restriction of the non linear system to the center manifold. In particular, system (28) restricted to the center manifold will show a super-critical Hopf bifurcation near $p_{eq}(\Omega_0)$ and $\Omega_0$.

**Proposition 2.1 (First order coefficients).** With the hypothesis of Theorem 2.2, the non-linear system (28) has a supercritical Hopf bifurcation near $(p_{eq}(\Omega_0), \Omega_0)$.

**Proof. Change of variables.**

In the first part of the proof we restrict the dynamical system to the center manifold. To compute the latter, we bring the system to a more suitable formulation.

Let us first transform the fixed point of (28) to the origin. We define the new variables as

$$ \delta p_1^c = p_1^c - p_{1,eq}^c, \quad \delta p_1^s = p_1^s - p_{1,eq}^s, \quad \delta q_1^c = q_1^c - q_{1,eq}^c, \quad \delta q_1^s = q_1^s - q_{1,eq}^s. $$
We then use the following linear and invertible change of variables

$$r = \frac{a_1}{b_1} \delta p_1^c + \delta p_1^s, \quad s = \frac{a_1}{b_1} \delta q_1^c + \delta q_1^s, \quad z = \frac{1}{2}(\delta p_1^s + \delta q_1^s), \quad y = \frac{1}{2}(\delta p_1^s - \delta q_1^s),$$

and shorten the notation by taking $X = (r, s, z)^T$, and $Y = (y, x)^T$ to transform the system (28) into

$$\frac{d}{dt} \begin{pmatrix} X \\ Y \end{pmatrix} = \mathbb{M}(\Omega) \begin{pmatrix} X \\ Y \end{pmatrix} + \begin{pmatrix} G(X, Y) \\ F(X, Y) \end{pmatrix}, \quad (30)$$

where

$$\mathbb{M}(\Omega) = \begin{pmatrix} -a_0 \mathbb{I}_{3,3} & 0 \\ 0 & \mathbb{A}(\Omega) \end{pmatrix}$$

is a block diagonal matrix with

$$\mathbb{A}(\Omega) = \begin{pmatrix} -a_0 - \frac{\lambda}{\ell} \frac{a_1}{a_0} & -\frac{\zeta}{\ell} \frac{a_1}{a_0} \\ -\frac{\lambda \ell}{2\pi} & -\frac{\zeta \ell}{2\pi} \end{pmatrix}, \quad (31)$$

and the non linear part is defined as

$$G(X, Y, \Omega) = (\zeta x + \lambda y) \begin{pmatrix} \frac{a_1}{b_1}(y+z) + \frac{b_1}{a_1}(-r+y+z) \\ \frac{a_1}{b_1}(y-z) + \frac{b_1}{a_1}(s+y-z) \\ -\frac{b_1}{a_1}\frac{(r-s-2y)}{2} \end{pmatrix},$$

and

$$F(X, Y, \Omega) = (\zeta x + \lambda y) \begin{pmatrix} -\frac{b_1}{a_1} \frac{(r+s-2z)}{2} \\ 0 \end{pmatrix}, \quad (32)$$

where we remind the reader that $a_0, a_1$ and $b_1$ depend on $\Omega$.

**Computation of the center manifold.**

The center manifold can then be computed by using standard techniques (see [29], Chapter 20, Section 2). We start by rewriting the system in such a way that $\Omega_0$ is moved to the origin through the change of variable $\delta\Omega = \Omega - \Omega_0$. As it is classical, we treat $\delta\Omega$ as a variable of the system. This means that we add the equation $\delta\dot{\Omega} = 0$ to the dynamical system and that the non linear part of the system includes all the products $\delta\Omega r, \delta\Omega s, \delta\Omega z$ etc.

Since the terms in the matrix $\mathbb{M}$ are nonlinear in $\Omega$, we expand them near $\Omega_0$ as for instance

$$a_0(\Omega_0 + \delta\Omega) = a_0(\Omega_0) + a'_0(\Omega_0)\delta\Omega + O(\delta\Omega^2). \quad (33)$$

Using such expansions in the system (30), we obtain

$$\frac{d}{dt} \begin{pmatrix} X \\ Y \\ \delta\Omega \end{pmatrix} = \mathbb{M}(\Omega_0) \begin{pmatrix} X \\ Y \\ 0 \end{pmatrix} + \begin{pmatrix} g(X, Y, \delta\Omega) \\ f(X, Y, \delta\Omega) \\ 0 \end{pmatrix}, \quad (34)$$
with

$$g(X, Y, \delta\Omega) = -a'_{0}(\Omega_{0})\delta\Omega X + G(X, Y, \Omega_{0}) + O(|\delta\Omega|^{3} + |X|^{3})$$

and

$$f(X, Y, \delta\Omega) = - \begin{pmatrix} \left(a'_{0}(\Omega_{0}) + \frac{\lambda}{\ell}\left(\frac{a_{1}}{a_{0}}\right)'(\Omega_{0})\right)\delta\Omega y + \frac{\zeta}{\ell}\left(\frac{a_{1}}{a_{0}}\right)'(\Omega_{0})\delta\Omega x \\ 0 \\ + F(X, Y, \Omega_{0}) + O(|\delta\Omega|^{3} + |X|^{3}) \end{pmatrix}$$

We can define a center manifold as

$$W^c(0) = \{(X, Y, \delta\Omega) : X = \mathbf{h}(Y, \delta\Omega), |Y| < \varepsilon, |\delta\Omega| < \bar{\varepsilon}, \mathbf{h}(0, 0) = 0, \nabla \mathbf{h}(0, 0) = 0\},$$
for $\varepsilon$ and $\bar{\varepsilon}$ sufficiently small and $\mathbf{h} = (h_1, h_2, h_3)$ smooth enough.

In order for $\mathbf{h}$ to be the center manifold for the system (34), the graph of $\mathbf{h}(Y, \delta\Omega)$ has to be invariant under the dynamics generated by (34). Hence, by plugging $\mathbf{h}$ into the system and compute its derivative, we get the following quasilinear differential equation

$$\nabla_Y \mathbf{h} \cdot (\mathbb{A}(\Omega_0)Y + f(\mathbf{h}, Y, \delta\Omega)) = -a_0(\Omega_0)\mathbf{h} + g(\mathbf{h}, Y, \delta\Omega). \quad (35)$$

The solution $\mathbf{h}(Y, \delta\Omega)$ of (35) can be approximated with a power series expansion up to any desired degree of accuracy. In our case we expand them up to order two defining, for $i = 1, 2, 3$ and $Y = (y, x)^T$

$$h_i(Y, \delta\Omega) := a_{i1}y^2 + a_{i2}yx + a_{i3}y\delta\Omega + a_{i4}x^2 + a_{i5}x\delta\Omega + a_{i6}\delta\Omega^2 + O(|\delta\Omega|^3 + |X|^3). \quad (36)$$

To completely determine (locally) the center manifold, we have to compute the coefficients $a_{ij}$ knowing that the functions defined in (36) must solve (35). For the detailed computations, we refer the reader to Appendix A.1.

Restricted to the center manifold, the original system (30) has the following form:

$$\frac{d}{dt}Y = \mathbb{A}(\Omega)Y + f(\mathbf{h}(Y), Y) + O(|\delta\Omega|^3 + |X|^3) \quad (37)$$

where

$$f(\mathbf{h}(Y), Y) = \left( \begin{array}{c} -\frac{2\pi a_0(\Omega_0) (\zeta^3 x^3 + \lambda^3 y^3) + \zeta \ell (\zeta x + \lambda y)^3}{2a_0(\Omega_0)\lambda(\pi a_0(\Omega_0) + 2\zeta\ell)} \\ 0 \end{array} \right) + O(|\delta\Omega|^3 + |X|^3).$$

Observe that $\tau(\Omega) \pm i\omega(\Omega)$ is the pair of conjugated eigenvalues of $\mathbb{A}(\Omega)$ and therefore we are in the hypothesis of the existence of a Hopf bifurcation for a two-dimensional system.

**Normal form.**
The next step is to bring system (37) into its normal form, from which we deduce the type of Hopf bifurcation that the system is attaining.

In order to proceed, we apply a further change of coordinates, such that $A(\Omega)$ is transformed to its real Jordan form. Namely, we define the transformation matrix

$$ \mathbb{P} = \begin{pmatrix} P_{11} & P_{12} \\ 1 & 0 \end{pmatrix}, $$

where

$$ P_{11} = \frac{\pi}{\lambda\ell} \left( \frac{\lambda}{\ell} \frac{a_1}{a_0} + a_0 \right) - \frac{\zeta}{2\lambda} $$

and

$$ P_{12} = \frac{\sqrt{-4\pi\lambda\ell a_0 a_1 (2\pi a_0 + \zeta\ell) - \ell^2 a_0^2 (\zeta\ell - 2\pi a_0)^2 - 4\pi^2\lambda^2 a_1^2}}{2\lambda\ell^2 a_0}. $$

The matrix $\mathbb{P}$ defines the new coordinates $Y = \mathbb{P}\tilde{Y}$ thanks to which equation (37) can be expressed in its normal form

$$ \frac{d}{dt} \tilde{Y} = \begin{pmatrix} \tau(\Omega) & -\omega(\Omega) \\ \omega(\Omega) & \tau(\Omega) \end{pmatrix} \tilde{Y} + \mathbb{P}^{-1} f(\mathbf{h}(\mathbb{P}\tilde{Y}), \mathbb{P}\tilde{Y}) + O(\|\delta\Omega\|^3 + \|Y\|^3). \quad (38) $$

To compute key properties of the system, such as the amplitude of the limit cycles, we change coordinates to the polar ones $\tilde{Y}^T = \rho(\sin(\theta), \cos(\theta))$, and get

$$ \left\{ \begin{aligned} \dot{\rho}(t) &= \tau'(\Omega_0)\delta\Omega\rho(t) + \tilde{\tau}\rho^3(t) + O(\delta\Omega^2\rho, \delta\Omega\rho^3), \\ \dot{\theta}(t) &= \omega(\Omega_0) + \omega'(\Omega_0)\delta\Omega + \tilde{\omega}\rho^2(t) + O(\delta\Omega^2, \delta\Omega\rho^2). \end{aligned} \right. \quad (39) $$

that is the normal form of (37) in polar coordinates around $\Omega_0$.

The non linear part $\mathbb{P}^{-1}f(\mathbf{h}(\mathbb{P}\tilde{Y}), \mathbb{P}\tilde{Y})$ determines the constants $\tilde{\tau}$ and $\tilde{\omega}$. The first one is involved in the expression for the limit cycle amplitude, and we compute it using a well know formula (see Appendix A.2). We obtain

$$ \tilde{\tau} = -\frac{3\pi\zeta}{4\ell} \left( \frac{\pi a_0(\Omega_0) + \ell\zeta}{\pi a_0(\Omega_0) + 2\ell\zeta} \right). $$

Since $\tilde{\tau}$ is negative and

$$ \tau'(\Omega_0) = - \frac{a_0'(\Omega_0) (\ell a_0(\Omega_0)^2 - \lambda a_1(\Omega_0)) + \lambda a_0(\Omega_0) a_1'(\Omega_0)}{2\ell a_0(\Omega_0)^2}, $$

is positive by hypothesis, the reduced system (37), and hence, the whole system (30), shows a supercritical Hopf bifurcation near the bifurcation parameter $\Omega_0$.

In particular, for $\Omega$ sufficiently near and greater then $\Omega_0$ there exists an asymptotically stable periodic orbit with radius and frequency as in (14) and (13), respectively.

This finishes the proof of Proposition 2.1

□
**Step 3.** Higher order terms. Lastly, we are going to prove that, after large times, the solutions of (26) $p_n^{c,s}(t)$ and $q_n^{c,s}(t)$, with $n > 1$ are periodic with the same period as $x(t)$. We remark again that, since the dynamics of $x$ (27) is independent of $p_n^c, p_n^s, q_n^c, q_n^s$ for $n > 1$, we may assume that $x(t)$ is given and periodic, and split the system into pairs of coupled equations as stated in the following Proposition.

**Proposition 2.2 (Higher order coefficients).** For $n > 1$ consider the infinite system of four equations:

$$ \left\{ \begin{aligned} \dot{p}_n^c(t) + \frac{2\pi n}{\ell} \dot{x} p_n^s(t) &= -a_0(\Omega) p_n^c(t) + a_n(\Omega)/\ell, \\ \dot{p}_n^s(t) - \frac{2\pi n}{\ell} \dot{x} p_n^c(t) &= -a_0(\Omega) p_n^s(t) + b_n(\Omega)/\ell. \end{aligned} \right. \qquad (40) $$

and

$$ \left\{ \begin{aligned} \dot{q}_n^c(t) - \frac{2\pi n}{\ell} \dot{x} q_n^s(t) &= -a_0(\Omega) q_n^c(t) + a_n(\Omega)/\ell, \\ \dot{q}_n^s(t) + \frac{2\pi n}{\ell} \dot{x} q_n^c(t) &= -a_0(\Omega) q_n^s(t) + b_n(\Omega)/\ell. \end{aligned} \right. \qquad (41) $$

Suppose that, for $t \ge 0$, the function $t \mapsto x(t)$ is periodic of period $T$. Then, the solutions to the systems (40) and (41) may have two behaviors: they either converge in time to periodic solutions of the same period $T$ if $n$ is odd, or they go to zero for large times when $n$ is even.

**Proof.** We introduce $z_n = p_n^c + ip_n^s$. The first two equations of (40) then become:

$$ \dot{z}_n + \left( a_0(\Omega) - \frac{2in\pi}{\ell} \dot{x} \right) z_n = c_n \qquad (42) $$

where $c_n = \frac{a_n(\Omega) + ib_n(\Omega)}{\ell}$ is constant. We thus deduce the following expression for $z_n(t)$

$$ z_n(t) = c_n \int_0^t e^{-\int_u^t (a_0(\Omega) - i\bar{v}(s))ds} du + z_n(0)e^{-\int_0^t (a_0(\Omega) - i\bar{v}(s))ds}, \qquad (43) $$

where $\bar{v} = \frac{2n\pi}{\ell} \dot{x}$.

Let us first notice that the second term, $z_n(0)e^{-\int_0^t (a_0(\Omega) - i\bar{v}(s))ds}$ goes to zero when $t$ goes to infinity. This, together with the fact that $a_n = b_n = 0$ if $n$ is even (see (11)), gives $c_n = 0$ and the solutions $p_n^{c,s}$ go to zero for large times when $n$ is even.

Let us now focus on the case where $n$ is odd. We want to prove that, when $\bar{v}$ is periodic of period $T$, $z_n$ converges towards a $T$-periodic solution after a transitory regime. We again notice that the second term of (43) converges to 0, and the result will follow by studying only the first term.
We thus denote by $\tilde{z}_n(t) = \int_0^t e^{-\int_u^t (a_0 - i\bar{v}(s))ds} du$, and compute

$$ 
\begin{aligned}
\tilde{z}_n(T+t) - \tilde{z}_n(t) &= \int_0^{t+T} e^{-\int_u^{t+T} (a_0 - i\bar{v}(s))ds} du - \int_0^t e^{-\int_u^t (a_0 - i\bar{v}(s))ds} du \\
&= \int_0^T e^{-\int_u^{t+T} (a_0 - i\bar{v}(s))ds} du + \int_T^{t+T} e^{-\int_u^{t+T} (a_0 - i\bar{v}(s))ds} du \quad (44) \\
&\quad - \int_0^t e^{-\int_u^t (a_0 - i\bar{v}(s))ds} du.
\end{aligned}
$$

Using the fact that $\bar{v}$ is $T$-periodic, the last two terms cancel and we obtain

$$ 
\begin{aligned}
\tilde{z}_n(T+t) - \tilde{z}_n(t) &= \int_0^T e^{-\int_u^{t+T} (a_0 - i\bar{v}(s))ds} du \\
&= e^{-(t+T)a_0} \int_0^T e^{ua_0 + \int_u^{t+T} i\bar{v}(s)ds} du
\end{aligned}
$$

which goes to 0 when $t$ goes to infinity.

We therefore proved that $z_n$ converges to a $T$-periodic function. The same result holds for $(q_n^c, q_n^s)$ solution to (41). $\square$

Finally, if we have initial conditions that are close to the equilibrium point of system (8), and an amount of ATP $\Omega$, which is close to $\Omega_0$, we know exactly how the solution of the system evolves in time. The first order coefficients of (24) go exponentially fast to the constant term $1/\ell$, the first order coefficients starts to oscillate in a limit cycle, and the higher order terms either go to zero or have the same patterns as $p_1^{c,s}$ and $q_1^{c,s}$ with same period of oscillation. Globally, the solution of (8), shows a supercritical Hopf-bifurcation in time, close enough to the parameter $\Omega_0$.

**Comment on hypotheses in Theorem 2.2.** We have shown that if there exists a value $\Omega_0 > 0$ such that $\tau(\Omega_0) = 0$ and $\tau'(\Omega_0) > 0$, then the system undergoes a supercritical Hopf bifurcation. The coefficients $a_0$ and $a_1$ have to be suitably chosen to satisfy these hypotheses.

In particular, to satisfy the first condition, the coefficient $a_1(\Omega_0)$ must be negative since it is given by:

$$ a_1(\Omega_0) = -\frac{\ell}{2\lambda} \left[ 2a_0^2(\Omega_0) + \frac{\zeta\ell}{\pi} a_0(\Omega_0) \right] < 0, \quad (45) $$

with $a_0(\Omega_0)$ positive.

Concerning the second condition, note that it has a physical meaning since increasing ATP concentration $\Omega$ drives the filament configuration from a stable
to an oscillatory regime. Notice that, if $\tau'(\Omega_0) < 0$, the system still undergoes a supercritical Hopf bifurcation, but in this case, the stable periodic orbit emerges as $\Omega$ decreases, which contradicts the hypothesis that $\Omega$ is the ATP concentration. This second condition is equivalent to

$$a'_{1}(\Omega_{0}) > -\frac{\zeta\ell^2}{2\lambda\pi}a'_{0}(\Omega_{0}), \qquad (46)$$

where equation (45) has been used to substitute $a_1(\Omega_0)$ in the derivative $\tau'(\Omega_0)$.

For example, if we take $a_0(\Omega) = c > 0$, then $a_1(\Omega_0) = -\frac{\ell}{2\lambda}\left[2c^2 + \frac{\zeta\ell}{\pi}c\right] = C$. We can then make a simple choice for $a_1(\Omega)$, taking an affine function $a_1(\Omega) = B(\Omega - \Omega_0) + C$, where $B$ is a constant. Depending on the sign of $B$ we are in the hypothesis of Theorem 2.2 (if $B < 0$) or not (if $B > 0$).

## 2.3 N-row model: theoretical considerations

In order to have a general idea of the system (10) dynamics, as in (24), we expand in Fourier series $Q_1, \dots, Q_N$ and treat order by order the Fourier coefficients of the probabilities. Then, to investigate the existence of a Hopf bifurcation we consider its linearization, as it was done with the two-row model. For the first order coefficients, we obtain a $3N - 1$ dimensional ODE system with unknowns:

$$p_j^{c,s}, \quad j = 1, \dots, N, \quad x_k, k = 1, \dots, N - 1.$$

We then perform some linear change of variables, similarly to (2.2), in order to write the Jacobian of the system as a block matrix. The first block is a $N + 1$ diagonal matrix with real and negative entries $-a_0(\Omega)$. The rest of the linearized system is defined by the following equations

$$\begin{pmatrix} \dot{q}_k \\ \dot{w}_k \end{pmatrix} = \begin{pmatrix} -\frac{2\pi^2}{\ell\eta}Up_e - a_0 & -\frac{2\pi}{\ell\eta}kp_e \\ -\pi U/\eta & -k/\eta \end{pmatrix} \begin{pmatrix} q_k \\ w_k \end{pmatrix} \quad (47)$$

for $k = 1, \dots, N - 1$, where $p_e = a_1/(a_0\ell)$, $q_k = p_k^s - p_{k+1}^s$ for $k = 1, \dots, N - 1$, and

$$\begin{pmatrix} w_1 \\ \vdots \\ w_{N-1} \end{pmatrix} = \begin{pmatrix} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & & \ddots & & \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_{N-1} \end{pmatrix}. \quad (48)$$

It follows that the Jacobian has $N + 1$ real and negative eigenvalues $-a_0(\Omega)$, and $N - 1$ identical pairs of complex and conjugated ones $\mu_k(\Omega)$ for $k = 1, \dots, N - 1$ which come from (47). We observe that $\mu_k(\Omega) = \tau(\Omega) + i\omega(\Omega)$, with $\tau$ and $\omega$ defined as in Theorem 2.2 by equations (12) and (29). Then, if there exists $\Omega = \Omega_0$ such that $\tau(\Omega_0) = 0$ and $\tau'(\Omega_0) > 0$, the eigenvalues
cross the imaginary axis all at the same bifurcation parameter. Thus, there is a suggestion of a bifurcation in the dynamics at $\Omega_0$, indicated by the linear part of the complete dynamical system. We will not investigate the theoretical aspects of this model further in this paper, as the Central Manifold Theorem used in previous section, does not apply to such systems. However, as shown in the next section, numerical simulations still suggest a potential pattern in the oscillations and hint that there is indeed a bifurcation.

### 3 Numerics

We now describe the numerical scheme used throughout the remaining of the paper. We only write it for $N = 2$, but it can easily be extended to $N$-row models, for $N > 2$ or restricted to the 1-row model.

We write a first-order upwind scheme for the densities $P$ and $Q$ where, after each step, we update the velocity $v$. Let us take $\Delta x$ such that $\ell = J\Delta x$ for $J \in \mathbb{N}^*$, and $\Delta t$ the time step. We define $P_j^n := P(j\Delta x, n\Delta t)$ and $Q_j^n := Q(j\Delta x, n\Delta t)$ for $n \ge 1$ and $1 \le j \le J$. Being $\ell$-periodic, we may extend $P_j^n$ for all $j \in \mathbb{Z}$ by setting $P_{j+J}^n = P_j^n$ for all $j \in \mathbb{Z}$ and all $n \in \mathbb{N}$; the same can be done for $Q_j^n$. At a fixed time step $t = n\Delta t$, $n \ge 1$, the velocity $v^n$ of the central tubule pair being known, we compute, for $j \in \mathbb{Z}$

$$ 
\left\{
\begin{aligned}
S_j^n &= -a_0 P_j^n + \omega_B(j\Delta x)/\ell, \\
T_j^n &= -a_0 Q_j^n + \omega_B(j\Delta x)/\ell, \\
P_j^{n+1} &= \begin{cases} \left(1 - v^n \frac{\Delta t}{\Delta x}\right) + v^n \frac{\Delta t}{\Delta x} P_{j-1}^n + \Delta t S_j^n, & \text{if } v^n > 0, \\ \left(1 + v^n \frac{\Delta t}{\Delta x}\right) - v^n \frac{\Delta t}{\Delta x} P_{j+1}^n + \Delta t S_j^n, & \text{if } v^n < 0, \end{cases} \\
Q_j^{n+1} &= \begin{cases} \left(1 + v^n \frac{\Delta t}{\Delta x}\right) - v^n \frac{\Delta t}{\Delta x} Q_{j-1}^n + \Delta t T_j^n, & \text{if } v^n > 0, \\ \left(1 - v^n \frac{\Delta t}{\Delta x}\right) + v^n \frac{\Delta t}{\Delta x} Q_{j+1}^n + \Delta t T_j^n, & \text{if } v^n < 0, \end{cases}
\end{aligned}
\right.
\quad (49) $$

and $(x^{n+1}, v^{n+1})$ are

$$ 
\left\{
\begin{aligned}
\eta v^{n+1} &= \Delta x \sum_{j=1}^{J} ((P_j^{n+1} - Q_j^{n+1}) \partial_{\xi} \Delta W(j\Delta x)) - kx^n, \\
x^{n+1} &= x^n + v^{n+1} \Delta t.
\end{aligned}
\right.
\quad (50) $$

Notice that, in (49), only one of the two updates is used for all $j$. Moreover, due to the stability of the upwind scheme, we need to check at each iteration that

$$ \left| v^n \frac{\Delta t}{\Delta x} \right| < 1. $$
## 3.1 Choice for parameters

In order to implement the numerical simulations, we have to characterize $\omega_A$ and $\omega_B$ by fixing their coefficients. For simplicity, we impose that $a_n = b_n = 0$ for $n > 1$ since, as observed in Theorem 2.2, they have no influence over the displacement's oscillations. We choose $b_1 \neq 0$ to make the transition rates not-symmetric, as in [12]. We opted for particularly simple polynomial functions, defining

$$a_0(\Omega) = c\Omega^{1/2}, \quad a_1(\Omega) = b_1(\Omega) = d\Omega,$$

assuming $c$ and $\Omega_0$ are chosen. We analytically recover $d$ by imposing $\tau(\Omega_0) = 0$

$$d = d(\Omega_0) = -\frac{c\ell}{2\pi\lambda}\left(2\pi c + \frac{\zeta\ell}{\sqrt{\Omega_0}}\right) \quad (51)$$

in such a way that the condition $\tau(\Omega_0) = 0$ is satisfied. By substituting $d$ into $\tau'(\Omega)$ we get the other condition $\tau'(\Omega_0) = \frac{\zeta\ell}{8\pi\Omega_0} > 0$.

Using the parameter values presented in Table 1, we get $d = d(\Omega_0) = -56.4588 s^{-1}$. Note that the parameters shown in Table 1, match F. Jülicher's [6]. We have taken for $k_B T$ the value at 36°C (human body temperature).

<table><caption>**Table 1:** Values of the parameters used in numerical simulations.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>ℓ</td><td>10 nm</td></tr><tr><td>k<sub>B</sub>T</td><td>4.2668 &middot; 10<sup>-3</sup> nN &middot; nm</td></tr><tr><td>η</td><td>1.0 &middot; 10<sup>-7</sup> kg/s</td></tr><tr><td>k</td><td>9.5 &middot; 10<sup>-5</sup> kg/s<sup>2</sup></td></tr><tr><td>U</td><td>10k<sub>B</sub>T</td></tr><tr><td>a<sub>0</sub><sup>0</sup></td><td>1.0 &middot; 10<sup>3</sup> s<sup>-1</sup></td></tr><tr><td>Ω<sub>0</sub></td><td>15k<sub>B</sub>T</td></tr></tbody></table>

In the following, all simulations are carried out defining $\Delta W(\xi) = U \cos(2\pi\xi/\ell)$.

## 3.2 2-row: numerical simulations

Using the upwind scheme (49-50) and the parameters chosen above, we run simulations for the two-row model. We compare its behavior to the one of the one-row model (1) for which the simulations are carried out using an obvious modification of the scheme.

When the ATP concentration $\Omega$ is lower than $\Omega_0$, the central pair does not move from its zero equilibrium position. On the contrary, when the ATP concentration is high enough, we observe an oscillatory displacement between microtubules and in the probabilities. As shown in Figure 5, the behavior of the two-row model with respect to the ATP concentration is comparable with the one for the one-row model. The main difference between the one-row model and the two-row one concerns the equilibrium position around which the solution
$x(t)$ oscillates, as expected. This is clearly illustrated in Figure 5(b), where the displacement takes place around zero, which is coherent with the symmetry of the physical problem.

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Two-row model ((a),(b)) and one-row model ((c),(d)): (a) Probability density $P_1(\xi, t)$ of the top layer of being in state 1 over time in a periodicity cell, in steady-state regime. The probability density $P_2(\xi, t)$ of the bottom layer of being in state 1 is similar to $P_1(\xi, t)$. (b) Relative tubule shift $x$ over time. (c) Probability density $P(\xi, t)$ of the moving layer in the one-row model. (d) Relative tubule shift $x$ over time in the one-row model.</figcaption></figure>

### 3.2.1 Amplitude of the oscillations

In this section, we study the evolution of the amplitude when the ATP concentration varies. In particular, we compare the *theoretical amplitude* given by the asymptotic formula (14) with the ones observed for the numerical solutions of the ODE system (28), solved using MATLAB's `ode45`, and of the full PDE model (8), computed with the upwind scheme. In view of the shape of the solution, to determine the *numerical amplitude* for both the ODE and the PDE, we calculate half of the difference between the maximum and minimum values of the displacement $x(t)$ along time, getting rid of the first 1000 time steps.

The parameter $\delta$ indicates the relative distance from the instability $\Omega_0$; therefore, the point where the simulations is performed is $\Omega(\delta) = \Omega_0(1 + \delta)$. For this choice of parameters and neglecting higher-order terms, the theoretical
amplitude (14) simplifies to:

$$ \rho(\delta) = \rho(\Omega(\delta)) = \sqrt{\delta \frac{\ell^2}{6\pi^2} \left( \frac{\pi a_0^0 \Omega_0^{1/2} + 2\ell\zeta}{\pi a_0^0 \Omega_0^{1/2} + \ell\zeta} \right)}. \quad (52) $$

In Figure 6, simulations only start after $\delta = 0.05$. This is related to the fact that the closer we are to $\Omega_0$, the slower the solution enters into its limit cycle.

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Two-row model. (a) Displacement amplitude vs. relative distance to instability: the red curve shows the amplitude from the ODE (28), the black curve from the PDE (8) and the theoretical amplitude from (52) in blue. (b) Relative error results between the theoretical amplitude and the PDE (in red), and the ODE and the PDE (in black).</figcaption></figure>

In Figure 6(a), we observe the amplitude of $x(t)$ computed with the PDE, the ODE, and the analytical result (52). As expected, the latter loses its prediction capability as we go away from the instability; and subsequently, the relative error between the formula (52) and the one measured from the PDE increases, as we can see in Figure 6(b). This is due to the contribution of the non linear terms, which play a role in determining the amplitude of oscillations even when the system is close to $\Omega_0$.

Note that the two-row model does not only center the equilibrium point for the oscillations around zero, but it also influences important physical quantities such as the amplitude of oscillations. This is emphasized in Figure 7 where the theoretical expansions for both models are plotted (the expansion (14) for the two-row model can be adapted to the one-row model).

### 3.2.2 Parameters sensitivity

In this section, we examine the sensitivity of the two-row model ODE (28) to the physical parameters $k$, $\eta$ and $\ell$, spring stiffness, viscosity and length of the slice, respectively. To assess their impact on the system, we analyze the oscillation
<figure><img src="image_8.png" /><figcaption>**Figure 7:** Theoretical amplitude for the two-row model computed with formula (14) in black and its equivalent for the one-row model, in red, against the distance from instability $\delta$.</figcaption></figure>

amplitude and frequency by varying one parameter while keeping the others fixed. The results, shown in Figure 8, are obtained using MATLAB's ode45 solver at $\Omega = \Omega_0(1 + \delta)$ with $\delta = 0.05$.

Our findings indicate that changes in $\ell$ mostly affect the oscillation amplitude, that goes from 0.3 to 0.8, while the frequency remains nearly constant, decreasing from 7.99 Hz to 7.98 Hz. Conversely, variations in the viscosity coefficient $\eta$ and the elastic coefficient $k$ change the frequency, with minimal impact on amplitude. Specifically, as $k$ increases, the amplitude shifts from 0.54 to 0.56, while the frequency rises from 60 Hz to 100 Hz. In contrast, increasing $\eta$ reduces the amplitude from 0.56 to 0.54 and lowers the frequency from 120 Hz to 60 Hz.

## 3.3 N-row model: Numerical results

In this section, the N-layer model is tested for $N = 8$, to match the 8 motor rows of an axoneme with a bridge and nine microtubule doublets, as in Figure 3. We use the same parameter values as for the previous systems, as specified in Table 1. The initial conditions are such that $Q_i(t = 0) = Q_{eq,i}$ for all $i$, and we destabilize the $x_i$s randomly in $[-\bar{x}, \bar{x}]$, where $\bar{x} = 0.01$ nm.

Using the same notation as before, i.e. $\Omega = (1+\delta)\Omega_0$, we measure the relative distance to the bifurcation point $\Omega_0$ using $\delta$.

Figure 9(a) shows the absence of oscillations before instability ($\Omega = 0.9\Omega_0$, i.e. $\delta = -0.1$), exactly as expected. Even though the system is put out of equilibrium by the initial conditions, since some tubule shifts are nonzero, they all quickly go back to zero and stop moving.

In Figure 9(b), we look at the system past the instability point, for $\Omega = 1.1\Omega_0$ i.e. $\delta = 0.1$. Here, the system clearly reaches a steady-state oscillating regime after a short amount of time. All oscillations also remain centered in zero. As expected from the linear study in section 2.2, all $\Delta_i$s have the same amplitude and oscillation frequency, but their phase difference varies depending on initial conditions. In the particular case shown in Figure 9(b), one group of tubule pairs is synchronized with pair 7, and the other group with pair 8. In fact, when looking at all layers separately, one can see that odd layers (respectively even layers)
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8:** Amplitude (first column) and frequency of oscillation (second column) of the two-row model ODE (28) by varying the physical parameters ($k, \eta, \ell$). The system was captured at the critical value of $\Omega = \Omega_0(1.05)$. In red $k, \eta, \ell$ are the model parameter shown in Table 1, while in blue all decreased and increased by 5% to perform sensitivity analysis.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 9:** N-row model (N = 8). Relative tubule displacement (in nm) with respect to time (in seconds) for random initial conditions (a) before bifurcation point, for $\delta = -0.1$ ($\Omega = 0.9\Omega_0$) and (b) and (c) after bifurcation point, for $\delta = 0.1$ ($\Omega = 1.1\Omega_0$), for two sets of random initial conditions.</figcaption></figure>

oscillate in sync. The other possible outcome when running the simulation with random initial conditions was groups {1, 2, 7, 8} and {3, 4, 5, 6} oscillating together, as shown in Figure 9(c). This influence of the initial conditions is understandable since there is no external force taking the whole axonemal structure and filament into account, resulting in limited coupling between layers at this level. The layers thus tend to keep their original phase difference. The synchronization shown in Figures 9(b) and 9(c) was predicted by Howard et al. [13], and proves the importance of having a bridge around which the system alternates between positive and negative displacement. In the general case, since there is no theory behind the bifurcation of this system of $3N - 1$ equations, formally understanding the coupling between phase difference and initial conditions remains an open question.

A similar problem has recently been numerically studied in Kuramoto oscillators [1]. The Kuramoto model is fairly different from ours but exhibits synchronization of coupled oscillators. Some first steps, only modeling individual motors along a single row, have been presented in [8].
# 4 Conclusion and outlook

In this paper, we propose a model that takes into account the cylindrical structure of an isolated axonome. Compared to the existing literature, this model generalizes previous model by [16] that possess only one row of molecular motors to the case of $N$, where $N$ represents the number of rows of motors walking between two microtubule doublets that are arranged in a circular configuration.

The first result on our model can be appreciated starting from the $N = 2$ case. The cylindrical structure answers the problem of the axonemal symmetrization proposed in [6], previously solved by choosing symmetric transition rates. Instead, in of this paper, general transition rates have been used all along. More in general, even in the $N$-row model, we observed that the equilibrium displacement with null external force was zero, independently from the potentials and the transition rates. This means that the desired symmetry of the axoneme is always preserved: at equilibrium, with no external forces, each microtubule has the same role and there is no initial shifting in the structure. We underline that this is made possible in our model without imposing specific symmetry for the potentials.

The problem appears as a integral partial differential equation for which we provide a mathematical framework to show its well posedness.

In the case where $N = 2$, we give the necessary and sufficient conditions for the system to get an Hopf bifurcation depending on the ATP concentration $\Omega$. Using the center manifold theory we have also shown that this bifurcation is supercritical. Physically it means that the molecular motors spontaneously oscillate under varying ATP concentrations around the axoneme.

We analyze the cases of $N = 2$ and $N = 8$ using both theoretical methods and numerical simulations. For the case $N = 8$, under appropriate initial conditions, we note the emergence of two out-of-phase groups among the microtubules, oscillating around the zero position with opposite shifting directions.

A natural extension to our work would be to take into account the presence of the central microtubule pair and, therefore, of the radial spokes coming from outer doublets towards the axonemal center. Moreover, this model can be thought of as a building block that, coupled with mechanics of an elastic flagellum, gives the feedback that generates oscillatory patterns following in the footsteps of [6, 7, 21, 28, 36] for the planar case and of [23] for the three dimensional case.

# Funding

This work was supported by a public grant as part of the Investissement d'avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH.
# A Details on Theorem 2.2

## A.1 Center manifold computation

From now on we define $a_{10}(\Omega) = a_1(\Omega)/a_0(\Omega)$ for simplicity.

Equation (35) is a system of three equations. We focus here on the first one, recovering the coefficients $a_{1j}$ of $h_1 = h_1(x, y, \delta\Omega)$ with $j = 1, \dots, 6$; the same can be done for the second and third equation, solving for the power series of $h_2$ and $h_3$.

The first equation reads

$$ (2a_{11}y + a_{12}x + a_{13}\delta\Omega) \left( \begin{aligned} & -a_0(\Omega_0)y - \frac{\lambda}{\ell}a_{10}(\Omega_0)y - \frac{\zeta}{\ell}a_{10}(\Omega_0)x \\ & -a_0'(\Omega_0)\delta\Omega y - \frac{\lambda}{\ell}a_{10}'(\Omega_0)\delta\Omega y - \frac{\zeta}{\ell}a_{10}'(\Omega_0)\delta\Omega x \\ & - \frac{b_1(\Omega_0)(h_1 + h_2 - 2h_3)}{2a_1(\Omega_0)}(\zeta x + \lambda y) \\ & + (a_{12}y + 2a_{14}x + a_{15}\delta\Omega) \left( -\frac{\beta\ell}{2\pi}y - \frac{\alpha\ell}{2\pi}x \right) \\ & + a_0(\Omega_0)h_1 + a_0'(\Omega_0)\delta\Omega h_1 \\ & - (\zeta x + \lambda y) \left( \frac{a_1^2(\Omega_0)(y + h_3) + b_1^2(\Omega_0)(-h_1 + y + h_3)}{a_1(\Omega_0)b_1(\Omega_0)} \right) \end{aligned} \right) = 0. $$

Since we only need terms $x, y$ and $\delta\Omega$ up to order two, we get rid of all the third order terms, getting

$$ (2a_{11}y + a_{12}x + a_{13}\delta\Omega) \left( \begin{aligned} & -a_0(\Omega_0)y - \frac{\lambda}{\ell}a_{10}(\Omega_0)y - \frac{\zeta}{\ell}a_{10}(\Omega_0)x \\ & + (a_{12}y + 2a_{14}x + a_{15}\delta\Omega) \left( -\frac{\beta\ell}{2\pi}y - \frac{\alpha\ell}{2\pi}x \right) \\ & + a_0(\Omega_0)h_1 - y(\zeta x + \lambda y) \left( \frac{a_1^2(\Omega_0) + b_1^2(\Omega_0)}{a_1(\Omega_0)b_1(\Omega_0)} \right) \end{aligned} \right) = 0. \quad (53) $$

We then proceed by matching the terms $x, y$ and $\delta\Omega$ with same order, and solve the resulting system of equations for $a_{1j}$ with $j = 1, \dots, 6$. Notice we could have kept only the zero order term in the Taylor expansion, since all the higher order terms disappear during this approximation.

From the first two equations we get the same sets of coefficients, namely we
obtain for $i = 1, 2$

$$a_{i1} = - \frac{\pi \lambda^2}{\ell(2\zeta\ell + \pi a_0(\Omega_0))} \frac{a_1^2(\Omega_0) + b_1^2(\Omega_0)}{a_0^2(\Omega_0)b_1(\Omega_0)}$$

$$a_{i2} = \frac{\zeta(\zeta\ell - \pi a_0(\Omega_0))}{(2\zeta + \pi a_0(\Omega_0))} \frac{a_1^2(\Omega_0) + b_1^2(\Omega_0)}{a_0(\Omega_0)b_1(\Omega_0)}$$

$$a_{i3} = 0, \quad a_{i4} = - \frac{\pi \zeta^2}{\ell(2\zeta\ell + \pi a_0(\Omega_0))} \frac{a_1^2(\Omega_0) + b_1^2(\Omega_0)}{a_0^2(\Omega_0)b_1(\Omega_0)}$$

$$a_{i5} = 0, \quad a_{i6} = 0.$$

Instead, for $i = 3$ we get

$$a_{31} = - \frac{\pi \lambda^2}{\ell(2\zeta\ell + \pi a_0(\Omega_0))} \frac{b_1(\Omega_0)}{a_0^2(\Omega_0)}$$

$$a_{32} = \frac{\zeta(\zeta\ell - \pi a_0(\Omega_0))}{(2\zeta + \pi a_0(\Omega_0))} \frac{b_1(\Omega_0)}{a_0(\Omega_0)a_1(\Omega_0)}$$

$$a_{33} = 0, \quad a_{34} = - \frac{\pi \zeta^2}{\ell(2\zeta\ell + \pi a_0(\Omega_0))} \frac{b_1(\Omega_0)}{a_0^2(\Omega_0)}$$

$$a_{35} = 0, \quad a_{36} = 0.$$

## A.2 Computation for $\tilde{\tau}$

We now consider system (38) at the bifurcation point

$$\frac{d}{dt}\tilde{Y} = \begin{pmatrix} 0 & -\omega(\Omega_0) \\ \omega(\Omega_0) & 0 \end{pmatrix}\tilde{Y} + \begin{pmatrix} 0 \\ \tilde{f}^t(x, y) \end{pmatrix} + O(|\delta\Omega|^3 + |Y|^3). \quad (54)$$

where

$$\begin{pmatrix} 0 \\ \tilde{f}^t(x, y) \end{pmatrix} = \mathbb{P}^{-1} f(\mathbf{h}(\mathbb{P}\tilde{Y}), \mathbb{P}\tilde{Y})$$

is the transformed non linear term.

The formula for $\tilde{\tau}$ can be found in [29] (Chapter 20, Section 2) and is the following

$$\tilde{\tau} = - \frac{1}{16 \omega(\Omega_0)} \left( \partial_{yy}^2 \tilde{f}^t + \partial_{xx}^2 \tilde{f}^t \right) \partial_{yx}^2 \tilde{f}^t + \frac{1}{16} \left( \partial_{yyy}^3 \tilde{f}^t + \partial_{xxy}^3 \tilde{f}^t \right),$$

where all the derivatives are evaluated at the bifurcation point, i.e (0, 0, $\Omega_0$).

## References

[1] A. Bayani, S. Jafari, H. Azarnoush, F. Nazarimehr, S. Boccaletti, and M. Perc. Explosive synchronization dependence on initial conditions: The minimal kuramoto model. *Chaos Soliton. Fract.*, 169:113243, 2023.
[2] P. V. Bayly and S. K. Dutcher. Steady dynein forces induce flutter instability and propagating waves in mathematical models of flagella. J. R. Soc. Interface, 13(123):20160523, 2016.

[3] C. J. Brokaw. Flagellar Movement: A Sliding Filament Model: An explanation is suggested for the spontaneous propagation of bending waves by flagella. Science, 178(4060):455-462, 1972.

[4] C. J. Brokaw. Bending patterns of ATP-reactivated sea urchin sperm flagella following high salt extraction for removal of outer dynein arms. Cell. Motil. Cytoskel., 42(2):125-133, 1999.

[5] C. J. Brokaw. Computer simulation of flagellar movement X: Doublet pair splitting and bend propagation modeled using stochastic dynein kinetics. Cytoskeleton, 71(4):273-284, 2014.

[6] S. Camalet and F. Jülicher. Generic aspects of axonemal beating. New J. Phys., 2(1):324, 2000.

[7] J. Cass and H. Bloomfield-Gadêlha. The reaction-diffusion basis of animated patterns in eukaryotic flagella. Nat. Commun., 14, 09 2023.

[8] G. Costantini and A. Puglisi. Thermodynamic precision of a chain of motors: the difference between phase and noise correlation. J. Stat. Mech.: Theory Exp., 2024(2):024003, feb 2024.

[9] E. Gaffney, H. Gadêlha, D. Smith, J. Blake, and J. Kirkman-Brown. Mammalian sperm motility: Observation and theory. Annu. Rev. Fluid Mech., 43(1):501-528, 2011.

[10] E. A. Gaffney, K. Ishimoto, and B. J. Walker. Modelling motility: The mathematics of spermatozoa. Front. Cell Dev. Biol., 9, 2021.

[11] J. Gray and G. J. Hancock. The propulsion of sea-urchin spermatozoa. J. Exp. Biol., 32, 1955.

[12] T. Guérin, J. Prost, and J.-F. Joanny. Dynamical behavior of molecular motor assemblies in the rigid and crossbridge models. Eur. Phys. J. E, 34:1-21, 2011.

[13] J. Howard, A. Chasteen, X. Ouyang, V. F. Geyer, and P. Sartori. Predicting the locations of force-generating dyneins in beating cilia and flagella. Front. Cell Dev. Biol., 10, 2022.

[14] T. Hu and P. V. Bayly. Finite element models of flagella with sliding radial spokes and interdoublet links exhibit propagating waves under steady dynein loading. Cytoskeleton, 75(5):185-200, 2018.

[15] F. Jülicher, A. Ajdari, and J. Prost. Modeling molecular motors. Rev. Mod. Phys., 69:1269-1282, 1997.
[16] F. Jülicher and J. Prost. Cooperative molecular motors. Phys. Rev. Lett., 75:2618-2621, 1995.

[17] F. Jülicher and J. Prost. Spontaneous oscillations of collective molecular motors. Phys. Rev. Lett., 78:4510-4513, 1997.

[18] F. Jülicher and J. Prost. Molecular Motors: From Individual to Collective Behavior. Prog. Theor. Phys. Supp., 130:9-16, 1998.

[19] C. B. Lindemann and K. A. Lesich. The many modes of flagellar and ciliary beating: Insights from a physical analysis. Cytoskeleton, 78(2):36-51, 2021.

[20] K. E. Machin. Wave Propagation along Flagella. J. Exp. Biol., 35(4):796-806, 1958.

[21] D. Oriola, H. Gadêlha, and J. Casademunt. Nonlinear amplitude dynamics in flagellar beating. R. Soc. Open Sci., 4(3):160698, 2017.

[22] E. M. Purcell. Life at low Reynolds number. Am. J. Phys., 45(1):3-11, 01 1977.

[23] P. Sartori, V. F. Geyer, J. Howard, and F. Jülicher. Curvature regulation of the ciliary beat through axonemal twist. Phys. Rev. E, 94:042426, Oct 2016.

[24] P. Sartori, V. F. Geyer, A. Scholich, F. Jülicher, and J. Howard. Dynamic curvature regulation accounts for the symmetric and asymmetric beats of Chlamydomonas flagella. eLife, 5:e13258, 2016.

[25] K. E. Summers and I. R. Gibbons. Adenosine triphosphate-induced sliding of tubules in trypsin-treated flagella of sea-urchin sperm. Proc. Nat. Acad. Sci. USA, 68(12):3092-3096, 1971.

[26] G. I. Taylor. Analysis of the swimming of microscopic organisms. Proc. R. Soc. Lond. A, 209:447-461, 1951.

[27] M. F. Velho Rodrigues, M. Lisicki, and E. Lauga. The bank of swimming organisms at the micron scale (boso-micro). PLoS One, 16(6):1-80, 06 2021.

[28] B. J. Walker, S. Phuyal, K. Ishimoto, C.-K. Tung, and E. A. Gaffney. Computer-assisted beat-pattern analysis and the flagellar waveforms of bovine sper-matozoa. R. Soc. Open Sci., 7(6), 2020.

[29] S. Wiggins. Springer, Texts in Applied Mathematics. 2003.

[30] L. G. Woodhams, Y. Shen, and P. V. Bayly. Generation of ciliary beating by steady dynein activity: the effects of inter-filament coupling in multi-filament models. J. R. Soc. Interface, 19(192):20220264, 2022.

[31] T. Yagi and R. Kamiya. Novel mode of hyper-oscillation in the paralyzed axoneme of a chlamydomonas mutant lacking the central-pair microtubules. Cell Motil., 31(3):207-214, 1995.
[32] H. Soares, B. Kamiya., S. Nolasco, L. Viseu Melo and J. Goncalves Cilia Distal Domain: Diversity in Evolutionarily Conserved Structures. *Cells*, 14;8(2):160, 2019.

[33] I. Riedel-Kruse, A. Hilfinger, J. Howard and F. Jülicher How Molecular Motors Shape The Flagellar Beat *HFSP journal*, Sep;1(3):192-208, 2007.

[34] G. Lattanzi and A. Maritan Force dependent transition rates in chemical kinetics models for motor proteins. *J. Chem. Phys.*, 117 (22): 10339-10349, 2002.

[35] G. Lattanzi, A. Maritan, Coarse grained models: the kinetics of motor proteins, *Computational Materials Science*, Volume 30, Issues 1-2, pages 172-179, 2004.

[36] I. Anello, F. Alouges and A. DeSimone Beating of eukaryotic flagella via Hopf bifurcation of a system of stalled molecular motors *Arxiv preprint*, 2024.