Some mathematical models for flagellar activation mechanisms
François Alouges Irene Anello Antonio DeSimone Aline Lefebvre-Lepot Jessie Levillain


Bias correction of posterior means using MCMC outputs
    Yukito Iba
    September 9, 2024
=====================================================




§ ABSTRACT 

    Several micro-swimmers in nature use active slender appendages such as cilia or flagella to swim inside fluids. The characteristic periodic motion of these cilia and flagella allows them to swim in viscous fluids at low Reynolds numbers. In order to generate this motion, some form of internal activity along the appendage is necessary.
This activity is created by an internal structure extending all along the flagellum, called the axoneme. It is composed of nine pairs of microtubules arranged on the surface of a cylinder. Protein macro-molecules, called molecular motors, are located between microtubule pairs and create local sliding between two neighboring pairs. The motors' work is powered by a chemical component called ATP. 
A “two-row model" is introduced, which is a 

simplified model of the axoneme with two motor rows in total, using a system of partial differential equations to describe the motors' behavior. The main theorems of this paper then prove existence and uniqueness of the solutions of this model is well-posed, and that the system undergoes a Hopf bifurcation depending on the quantity of ATP available to the motors. Near the instability, one can also compute the oscillation amplitude.
The two-row model is then extended to a system with an arbitrary number N of motor rows, and
all models are also studied through numerical simulations, in which one can observe periodic oscillations 

past a certain ATP value, without any additional form of regulation in the system.

This paper focuses on studying a model for molecular motors responsible for the bending of the axoneme in the flagella of microorganisms. The model is a coupled system of partial differential equations inspired by <cit.>, incorporating two rows of molecular motors between microtubules filaments. Existence and uniqueness of a solution is proved, together with the presence of a supercritical Hopf bifurcation. Additionally, numerical simulations are provided to illustrate the theoretical results. A brief study on the generalization to N-rows is also included.









 §.§ Division of paragraphs


    
  * Irene: 
    
        
        
        
        
        
  * 15/12: General form of the theorem (old file), 4.3.1
        
  * have to do : 3.5 (only figures), find a link between obs on b1 e comparison one-row two-row (notes to bring at the meeting)
        
  * 15/01 Rewrite theorem and appendix
    
    
  * Jessie:
    
        
  * 1-4/12: I've modified the intro
        
  * I need to write Appendix A at some point
        
  * 4/12: I wrote paragraph 4.1
        
  * 7/12: I wrote 4.2. I've added numerical results for 4.4.2 and 4.5. I'm writing descriptions
    
    
  * Remaining paragraphs:
    
        
        
  * 3.6.2 tug of war 
        
  * Simulations Appendix A 
        
  * 3.5 (only figures)
        
  * 4.3.2
        
  * all of 4.4
        
        
  * conclusion
        
  * abstract
    


where did we say how oscillation amplitude was computed?
In 3.7.1



 §.§ Molecular motors bibliography

cited:
<cit.>, <cit.>,  <cit.>, <cit.>, <cit.> <cit.>, <cit.>,   <cit.>, <cit.>. <cit.>, <cit.>,
General review: <cit.>.
<cit.>

Review <cit.>. 

Curvature regulation models
Lindemann's geometric clutch model
<cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>. 

Steady dynein loading: <cit.>, <cit.>, <cit.>.
inhibition/excitation <cit.>

Filament structure, simulations
<cit.>
theory and simulation by Brokaw: <cit.>, <cit.>,<cit.>.
Sliding filament theory (muscle contraction, A. Huxley, R. Niedergerke, H. Huxley, J. Hanson). Review by J. M. Squire: <cit.>.
Optimal swimming (Pironneau): <cit.>.

Motor protein types and models: <cit.>.

Deriving molecular models from Langevin equations
Read papers by H. Wang, on mathematical theory for molecular motors:  <cit.>.

Predicting location of force-generating dyneins: <cit.>.



 §.§ not cited yet


numerically solving Fokker-Planck equations: <cit.>, <cit.>

Spontaneous cilia synchronization
Spontaneous cilia coordination: <cit.> 
determining an optimal number of motors <cit.>

Misc?
<cit.>, 
tracking kinesin jumps: <cit.>




§ INTRODUCTION



Cilia and flagella are active 

slender organelles employed by 

eukaryotic cells to swim through viscous fluid <cit.>. They exploit hydrodynamic friction to induce 

self-propulsion and show a characteristic non-reciprocal periodic motion in their dynamics <cit.>. This enables them to move without relying on inertial forces, i.e.,

at a low Reynolds number <cit.>. 
Pioneering works ( <cit.>) 
established that some form of internal activity along the flagellum is necessary in order to have sustained oscillations. 


    Pioneering work by Taylor <cit.>, Gray and Hancock <cit.>, and
Machin <cit.> 
established that actuation at one end is not sufficient to establish a traveling wave of bending similar to the one experimentally observed, and that 
some form of internal activity along the flagellum is necessary to sustain the oscillations. 
Numerical simulations with distributed activity were then conducted by Brokaw <cit.>, by matching parameters values in order to reproduce a flagellum's beating pattern <cit.>.



This internal activity is present in all cilia and flagella due to their common internal structure, the axoneme. The axoneme consists of nine pairs of protofilaments called microtubules. It is an active structure due to the molecular motors located between neighboring microtubule pairs. These molecular motors are attached to one microtubule pair and walk along the adjacent pair, creating local sliding between contiguous microtubule pairs.



This mechanism, which is fueled by the conversion of chemical energy (ATP) into mechanical work, creates bending along the filament <cit.>. 







       






The three-dimensional structure of the axoneme can be modeled as two planar filaments, with motors attached to the top filament exerting force on the fixed bottom filament (<cit.>,<cit.>,<cit.>). Since we are interested in the emergence of mesoscopic order in the axoneme from the uncoordinated activity of many independent molecular motors, we consider a portion of the filaments of arc-length of the order of 10 ℓ, where ℓ is a sub-micron length scale defined by the periodic microtubule structure (specifically, the typical inter-dynein spacing is 24 nm, see <cit.>). Over this length scale, the filaments can be considered as rigid. 

Throughout this paper, we will call the one-row model the system composed of a single row of motors between two filaments. In this work, we extend the one-row model by introducing a second row of molecular motors, as in Figure <ref>(b). We call this new framework, which symmetrizes the previous model, the two-row model. We also present a more comprehensive N-row model, where N corresponds to the number of motor rows.



    The two-row model is a simplified version of the axoneme, which forms a cylinder with only two microtubule pairs linked by two motor rows, as shown in Figure <ref>(a). It shall also be seen as a first step to the real model, which we are going to study via an N-row model in section <ref>, where N corresponds to the quantity of motor rows (see Figure <ref>). Since we are interested in the microscopic description of the axoneme, we assume rigid filaments all the N-row models.
   








This complex three dimensional structure has first been modeled in <cit.> as a pair of filaments that bends in a plane; throughout this paper we will call this model the one-row model. This kind of symplified model of the axoneme, further discussed in <cit.>, implies a great simplification on the study of the dynamics of the system, while still 
exhibiting
the typical oscillatory behaviour of cilia and flagella. 

In particular, it is possible to 
couple explicitly the dynamics of active filaments with the stochastic equations that regulate the molecular activity of the dyneins and still keep the governing equations at a manageable level of complexity. Our focus is going to be on this last part, and we study an amended
simplified model filament with additional motors rows while still using the stochastic model presented in <cit.> and <cit.>. Here motors are attached to one of  two filaments 
and walk on the other one. This model has no additional components regulating sliding, and nevertheless proves that dyneins can synchronize and thus auto-regulate themselves to create a form of alternating sliding on both sides of the axoneme. This matches, in particular, with the curvature regulation hypothesis called steady dynein loading <cit.> in which axonemal curvature is created without any inhibition mechanism transmitted along the system. 

Other curvature regulation models, that use additional mechanical signaling to provide some form of inhibition have been proposed in the literature <cit.>.





To study the two-row model, we generalize the stochastic model proposed in <cit.>. Each motor is anchored to one of two filament pairs, pair 1 or pair 2, and has two identical heads, e.g. head A and head B, each of them being either bound or unbound to the opposite filament. The motor has two different chemical states, state 1, where head A is bound and head B is unbound, or state 2, where the opposite holds. In this model, both heads cannot be unbound (or bound) at the same time. Between the pair of filaments, passive elastic and viscous elements resist the motion; they are modeled by the positive constants k and η, respectively.


 



Each state i = 1, 2 is defined by its potential energy W_i(ξ) at position ξ∈ℝ. The transition rates, ω_i(ξ)= ω_ij(ξ) with i,j = 1,2 and i ≠ j, represent the probability per unit time for a motor to switch between state i to state j. Moreover, since the filament is a periodic structure, W_i and ω_i are periodic functions with period ℓ.
 As ATP provides the energy for the motor to change states, the transition rates depend on its concentration Ω. A detailed explanation of this is available in <cit.>.


In the limit of an infinite system of motors, we introduce P_i(ξ,t) as the probability density that a motor anchored to tubule pair i is in state 1 at position ξ and time t.
To obtain the probability density that a motor anchored to tubule pair i is in state 2, one must compute 1/ℓ - P_i. 
The shift between the moving and the non-moving filament is measured by x(t), with velocity v(t) = d/dtx(t).

   
Let P(ξ, t) = P_1(ξ,t) and Q(ξ,t) = P_2(ξ + x(t),t).
Then, we obtain (see section <ref>) the following system of equations, for ξ∈ [0, ℓ] and t>0:


    {[                                                                            ∂ P/∂ t(ξ,t) + v(t) ∂ P/∂ξ(ξ, t)  = -( ω_1(ξ) +ω_2(ξ)) P(ξ, t) + ω_2(ξ)/ℓ,;                                                                              ∂ Q/∂ t(ξ,t) - v(t) ∂ Q/∂ξ(ξ,t)  = -( ω_1(ξ) +ω_2(ξ)) Q(ξ,t) + ω_2(ξ)/ℓ,;                                                                                           v(t) = 1/2η( ∫_0^ℓ (P(ξ,t)- Q(ξ,t))∂_ξΔ W (ξ)dξ - 2kx(t) ),;                                                                                                                  P(0,t) = P(ℓ,t),    Q(0,t) = Q(ℓ,t), ].

where ω_i = ω_i(ξ; Ω) and Δ W (ξ) = W_2(ξ) - W_1(ξ) are ℓ-periodic.
As one can see in Figure <ref>(b), pair 0 and pair 2 are the same and they are both fixed. 


The goal of this paper is twofold. On the one hand, we rigorously prove an existence and uniqueness result for the solution of (<ref>). Furthermore, we show the existence of a Hopf bifurcation for that same system when reasonable formulas are taken for potentials W_i and transition rates ω_i, depending on the ATP concentration via Ω. 

Moreover, we provide the reader with numerical simulations of (<ref>), where an upwind scheme is employed to solve both transport equations. The numerical method is then extended to the more realistic N-row model.

Before stating the theoretical results, we point out that we assume (see <cit.>),

    ω_1(ξ; Ω) +ω_2(ξ; Ω) = a_0(Ω),

where a_0(Ω)>0. Let us define C_#^1([0,ℓ]) as the set of restrictions to [0,ℓ] of C^1 and periodic functions over ℝ, then the first theorem reads as follows:
 
Let us fix ℓ>0. Assume ω_1(ξ;Ω) and ω_2(ξ;Ω) as in (<ref>). Moreover, assume ω_2 and Δ W to be at least C_#^1([0,ℓ]).
If the initial data P(ξ, 0) and Q(ξ, 0) are C_#^1([0,ℓ]) and x(0)=0, then the system of equations (<ref>), admits a unique solution P,   Q ∈ C^1([0,ℓ]×ℝ_+), with ξ↦ P(ξ,·) and ξ↦ Q(ξ,·) in C_#^1([0,ℓ]), and x ∈ C^1(ℝ_+).


For the second result we chose Δ W (ξ) = U cos(2 πξ/ℓ), as in <cit.>. Moreover, since the two motor heads are identical, then ω_2(ξ) = ω_1(ξ + ℓ/2). Taking the transition rates' periodicity into account, we use the following Fourier expansion

    ω_2(ξ; Ω) = a_0(Ω)/2 + ∑_n = 2k + 1,  k ≥ 0( a_n(Ω) cos2n πξ/ℓ + b_n(Ω)sin2n πξ/ℓ),

with a_n(Ω) and b_n(Ω) real coefficients. 
 In addition to the hypothesis above, assume that a_0(Ω)=a_0^0Ω^α_0, a_1(Ω)=a_1^0Ω ^α_1, b_1(Ω)=b_1^0Ω^α_1 for a_1^0,b_1^0 ∈ℝ and α_0,α_1,a_0^0, Ω∈ℝ_+. Furthermore, let us define 

    τ (Ω) := -1/4(2a_0^0Ω^α_0+ ζℓ/π+  2 λa_1^0/a_0^0 ℓΩ^α_1-α_0),

with ζ = 2 π k / ηℓ and λ = 2 π^2 U / ηℓ.
Suppose there exists Ω_0 ∈ℝ_+ such that τ(Ω_0)=0 and τ'(Ω_0)>0, then the solutions P(ξ,t;Ω), Q(ξ,t;Ω) and x(t;Ω) of the system (<ref>) show a super-critical Hopf bifurcation in time near the bifurcation value Ω_0 and near the fixed point 
    
    (P_eq(ξ;Ω),Q_eq(ξ;Ω),x_eq(Ω)) = (ω_2(ξ;Ω)/a_0(Ω)l,ω_2(ξ;Ω)/a_0(Ω)l,0).

    In particular, the fixed point is asymptotically stable for Ω < Ω_0, and unstable for Ω > Ω_0. Moreover there exists an asymptotically stable periodic orbit for Ω > Ω_0 with radius 
 
   


    ρ = √(-(Ω-Ω_0)τ'(Ω_0)/τ̃) +o(√(Ω-Ω_0)),

where τ̃ = - 3 πζ/4 ℓ(π a_0^0Ω_0^α_0 + ℓζ/π a_0^0Ω_0^α_0 + 2 ℓζ), τ̃<0.


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.


The next section presents the two-row model and proves Theorems <ref> and <ref>, with numerical simulations confirming these results. The final section explores new questions by extending the model to an N-row structure.



§ INTRODUCTION2

Cilia and flagella are slender, active organelles that enable eukaryotic cells to swim through viscous fluids by exploiting hydrodynamic friction. They exhibit periodic motion, allowing movement without relying on inertial forces at low Reynolds numbers. All cilia and flagella share a common internal structure, the axoneme, consisting of nine pairs of microtubules.

The axoneme's activity is driven by molecular motors between neighboring microtubule pairs. These motors, attached to one microtubule pair, walk along the adjacent pair, creating local sliding by converting chemical energy (ATP) into mechanical work, causing the filament to bend.

This study focuses on a microscopic segment of the axoneme, considering the filaments rigid. The axoneme, that is therefore a complex three-dimensional structure, is modeled as a pair of filaments sliding relative to each other due to motor activity, referred to as the one-row model.

We introduce a two-row model of the axoneme, depicted in Figure 1, by adding another row of motors and a more comprehensive N-row model, where N represents the number of motor rows. The two-row model is studied using the two-state model proposed in [10]. Each motor is attached to one of the two filaments (pair 1 or pair 2) and can be either in state 1 (bound) or state 2 (unbound), with state-dependent potential energy W_i(ξ). Transition rates ω_1(ξ) and ω_2(ξ)  indicate the probability of transitioning between states, dependent on ATP concentration Ω, as detailed in [9].

In the infinite motor limit, Pi(ξ, t) denotes the probability of a motor on tubule pair i being in state 1. The probability of being in state 2 is 1/ℓ - Pi. The moving filament's position is x(t) with velocity v(t) = dx/dt. Defining P (ξ, t) = P1(ξ, t) and Q(ξ, t) = P2(ξ + x(t), t), we derive the following system (see section 3).




§ NUMERICAL SIMULATIONS OF THE ONE-ROW MODEL

In this section, we describe the numerical scheme used throughout the whole paper. It is used here to compute the solution to the one-row model, but can easily be extended to N-row models, for N>1.
We also show that well-known results for the 1-row model are recovered using our numerical scheme.






    
    















    




 §.§ One-Row model

In the one-row model the motors are attached to the top filament only and they exerts the force on the bottom and non moving filament. 

We introduce densities P^i(ξ)≥ 0 with i= 1,2 which give the probability to find a motor at position ξ in state i. These densities are not independent but obey P^1(ξ) + P^2(ξ) = 1/l, reflecting the fact that a motor is always either in state 1 or in state 2. The top filament is at position x(t), with velocity v(t) = d/dtx(t).
The one-row model, given in <cit.>, expresses this motion as a transport equation on P=P^1, namely:

    {[ ∂_t P(ξ,t) + v(t)∂_ξ P(ξ,t) = -( ω_1 +ω_2) P(ξ,t) + ω_2/ℓ;             v(t) = 1/η(∫_0^ℓ dξ P(ξ,t)∂_ξΔ W(ξ) - kx(t)), ].
 
where ω_i = ω_i(ξ; Ω) and Δ W (ξ) = W_2(ξ) - W_1(ξ). 

This system has a non moving and stationary solution:

    P_eq(ξ; Ω) = ω_2/ℓ(ω_1+ω_2),   v_eq = 0,   x_eq = 1/k∫_0^ldξ P_eq(ξ; Ω)∂_ξΔ W(ξ).


Depending on the chosen expressions for transitions rates ω_i and potentials W_i, a threshold value Ω=Ω_0 may exist, above which this solution becomes unstable and the system oscillates.




 §.§ Numerical method

For numerical simulations of the whole 1-row PDE system (<ref>), we write a first-order upwind scheme for the density P where, after each step, we update the velocity v. Let us take Δ x such that ℓ = J Δ x for J ∈ℕ^*, and Δ t the time step. We define P_j^n := P(j Δ x, n Δ t), for n≥ 1 and 1 ≤ j ≤ J. Being ℓ-periodic, we may extend P_j^n for all j ∈ℤ by setting P_j+J^n=P_j^n for all j ∈ℤ and all n ∈ℕ. At a fixed time step t= n Δ t, n≥ 1, the velocity v^n of the central tubule pair being known,  we compute, for j ∈ℤ

    {[                           S^n_j = -a_0 P^n_j + ω_2(jΔ x)/ℓ,                                                            ; P^n+1_j = ( 1 - v^nΔ t/Δ x) + v^nΔ t/Δ xP^n_j-1 + Δ t S^n_j                                                 if  v^n >0,; P^n+1_j = ( 1 + v^nΔ t/Δ x) - v^nΔ t/Δ xP^n_j+1 + Δ t S^n_j                                                 if  v^n <0. ].

and (x^n+1,  v^n+1) are

    {[ η v^n+1 = Δ x ∑_j=1^J (P^n+1_j ∂_ξΔ W(jΔ x)) - kx^n,;                              x^n+1 = x^n + v^n+1Δ t. ].

Notice that, in (<ref>), 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 | v^n Δ t/Δ x| <1.



 §.§ Numerical results


All simulations are carried out defining the model functions that satisfy the hypotheses from theorem <ref>: Δ W(ξ) = U cos (2 πξ / ℓ), a_0(Ω) =a_0^0 √(Ω), a_1(Ω)=b_1(Ω) =a_1^0Ω and a_n(Ω)=b_n(Ω) = 0, if n>1, with α_0 = 1/2 and α_1=1. All parameters values used for all simulations (unless specified otherwise) are shown in Table 1, matching F. Jülicher's <cit.>. We have taken for k_BT the value at 36^∘C (human body temperature). We compute a_1^0 such that the threshold of value of the system remains Ω_0, which is given.

These choices will be justified in section <ref>.


As written above, in <cit.>, the transition rates were taken symmetric to reflect the symmetry of the three dimensional body. Namely, ω_is are symmetric  in the sense of <cit.> if |ω_i(ξ+l/2)|=|ω_i(-ξ+l/2)|.
To point out what happens if the system is not symmetric we choose b_1 ≠ 0, as in <cit.>.


 

The system needs a specific ATP concentration Ω_0 to start moving. This can be observed in Figure <ref>. Here, we plot the graphs of P(ξ,t) in the periodicity cell [0,ℓ] for the space variable ξ. 
When Ω<Ω_0, the system quickly goes to its equilibrium state. The filament does not move and maintains a non-zero displacement, x_0 ≠ 0, as shown in Figures <ref>(a), <ref>(c) and  <ref>(e). In this regime, the quantity P(ξ,t), which is the probability of finding a motor with the right head attached to the filament and the left head detached from it, is constant in time. Lastly, due to our choices for the transition rates, the natural equilibrium configuration of the motor's heads is not independent from the position in space.



As soon as the system with elastic links attains enough chemical energy (namely, Ω > Ω_0), the position of the filament, together with the motor density P, oscillate around their equilibrium positions, as in Figures <ref>(d) and <ref>(f).  

Notice that, since the transition rates are not symmetric, the equilibrium displacement for the one-row model is different from zero (Figures <ref>(c) and <ref>(d)). 
    
It is also known that the dynamical behavior of the filament depends on the presence of elastic elements. They are defined through the elastic force -kx(t), with k ≠ 0. If enough chemical ATP is given to the system, namely Ω > Ω_0, then the filament is expected to deviate from its equilibrium position. If no elastic elements are present, i.e. k = 0, then the motors walk and the filament slides, with v(t) ≠ 0, in one arbitrary direction. In this case, | x | keeps growing until both filaments are not connected to one another (through the motors) anymore. Instead, if k ≠ 0, the filament starts to oscillate. 


    



§ THE TWO-ROW MODEL



 §.§ Motivation

In this section we focus on the two-row model, from both theoretical and numerical perspectives. 



The two-row model serves the purpose of symmetrizing the one-row model.
In the one-row model <cit.>, 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 <cit.>, solving it by choosing symmetrical potentials, as in our case, and symmetrical transition rates, namely |ω_i(ξ+l/2)|=|ω_i(-ξ+l/2)|, with ξ∈ [0,ℓ]. Instead, in this work, we model non-symmetrical transition rates and introduce a second row of molecular motors, producing a zero equilibrium x_eq=0 at rest.



The two-row model is a simplified version of the axoneme that acts like its two-dimensional projection (<cit.>), forming a cylinder with only two microtubule pairs linked by two motor rows, as shown in Figure <ref>(a). Moreover, this model is a first step to the N-row model presented in section <ref>, where N corresponds to the quantity of motor rows. 



Despite lacking additional sliding regulation components, the model demonstrates that dyneins can synchronize and self-regulate to create alternating sliding on both sides of the axoneme. This matches, in particular, with the curvature regulation hypothesis called steady dynein loading <cit.> in which axonemal curvature is created without any inhibition mechanism transmitted along the system.




 §.§ Model structure

We now derive in detail the system (<ref>), which is depicted in Figure <ref>(b). We consider three pairs of filaments, numbered from 0 to 2. The pairs 0 and 2 are identical, and fixed, while the central pair 1, may move. We call x(t) its displacement, and v(t)= d/dtx(t) its velocity. We define P_1(ξ,t) and P_2(ξ,t) as the probabilities for the motors to be in state one at position ξ 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, as in (<ref>), and obtain two transport equations for ξ∈ [0, ℓ] and t>0:

    {[                                                                       ∂ P_1/∂ t(ξ,t) + v(t) ∂ P_1/∂ξ(ξ,t) = -a_0(Ω) P_1(ξ,t) + ω_2(ξ; Ω)/ℓ,;                                                                                       ∂ P_2/∂ t(ξ,t)  = -a_0(Ω)P_2(ξ,t) + ω_2(ξ-x(t); Ω)/ℓ. ].


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

    f_ext-2η v - 2kx +f_mot = 0,

where f_ext is the external force applied on the central filament, and f_mot is the active force exerted by the motors, which is given by
    
    f_mot(t) = ∫_0^ℓ (P_1(ξ,t)∂_ξΔ W (ξ) - P_2(ξ,t)∂_ξΔ W (ξ - x(t)))dξ.

Here, k x and η v represent the elastic and viscous resistances, respectively.
For the purpose of our work, we do not consider any external forcing, which means
f_ext = 0.
Let P(ξ, t) = P_1(ξ,t) and Q(ξ,t) = P_2(ξ + x(t),t). Then, (<ref>) becomes 

    {[                                                               ∂ P/∂ t(ξ,t) + v(t) ∂ P/∂ξ(ξ,t)  = -a_0(Ω) P(ξ,t) + ω_2(ξ)/ℓ,;                                                               ∂ Q/∂ t(ξ,t) - v(t) ∂ Q/∂ξ(ξ,t)  = -a_0(Ω) Q(ξ,t) + ω_2(ξ)/ℓ. ].

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

    f_mot(t) = ∫_0^ℓ (P(ξ,t)- Q(ξ,t))∂_ξΔ W (ξ)dξ ,

and obtain, from (<ref>),

    v(t) = ẋ(t) = 1/2η( ∫_0^ℓ (P(ξ,t)- Q(ξ,t))∂_ξΔ W (ξ)dξ - 2kx(t) ).


Finally, we recover the system of equations (<ref>), which will then be studied both theoretically and numerically, before generalizing it to an arbitrary number of layers N.




 §.§ Proof of Theorem <ref> (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 ℝ_+.


Step 1: local existence.
Let T>0 be a time to be chosen afterwards and define the map ψ: C^0([0,T]) → C^0([0,T]) as ψ[x(·)](t) = ∫_0^t F[x(·)](s) ds, where 

    F[x(·)](s) = 1/2η∫_0^ℓ (P_x(ξ,s)- Q_x(ξ,s))∂_ξΔ W (ξ)dξ - k/η x(s),

and the functions P_x(ξ,t) and Q_x(ξ,t) are defined as the solutions of (<ref>) and v(t) = ẋ(t). Notice that P_x and Q_x are explicitly given in their integral form by

    P_x(ξ,t)= e^-a_0 tP (ξ - x(t)+ x(0), 0) +e^-a_0 t/ℓ∫_0^t e^a_0 sω_2(ξ + x(s)-x(t)) ds,

and 

    Q_x(ξ,t)=  e^-a_0 tQ (ξ +x(t)-x(0), 0) +e^-a_0 t/ℓ∫_0^t e^a_0 sω_2(ξ - x(s)+x(t)) ds,

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


Notice that (x̃, P_x̃, Q_x̃) is a solution of the system (<ref>, <ref>) if and only if x̃(t) is a fixed point of ψ, i.e. ψ[x̃(·)] = x̃(·). In order to proceed, we prove that ψ: C^0([0,T]) → 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 (<ref>, <ref>). The initial conditions are identical P_x_1(ξ, 0) = P_x_2(ξ, 0)=P(ξ) and Q_x_1(ξ, 0) = Q_x_2(ξ, 0)=Q(ξ). We want to estimate the quantity 

    ψ[x_1(·)](t)-ψ[x_2(·)](t) = ∫_0^t F[x_1(·)](s)-F[x_2(·)](s) ds.

Using (<ref>), we get

    F[x_1(·)](s)   -F[x_2(·)](s) = 1/2η( ∫_0^ℓ(P_x_1(ξ, s) - P_x_2(ξ, s)  
        - ( Q_x_1(ξ, s) - Q_x_2(ξ, s) ))∂_ξΔ W (ξ)dξ)- k/η (x_1(s)-x_2(s)).

Now, we have from (<ref>)

    P_x_1(ξ,t)     - P_x_2(ξ,t)=  e^-a_0 t(P_x_1 (ξ - x_1(t), 0) - P_x_2 (ξ - x_2(t), 0))
       +e^-a_0 t/ℓ∫_0^t e^a_0 s(ω_2(ξ + x_1(s)-x_1(t))-ω_2(ξ + x_2(s)-x_2(t)) ds.



We have the estimate

    P_x_1(t,ξ) - P_x_2(t,ξ)_L_t^∞L_ξ^∞
       ≤∂_ξ P_L_ξ^∞x_1(t)-x_2(t)_L_t^∞
            +1/ℓ∫_0^t e^a_0 (s-t)∂_ξω_2_L_ξ^∞ | x_1(s)-x_1(t)- x_2(s)+x_2(t))| ds _L_t^∞
       ≤( ∂_ξ P_L_ξ^∞ + 2/a_0 ℓ∂_ξω_2_L_ξ^∞) x_1(t)-x_2(t)_L_t^∞,

The same goes for the difference Q_x_1(ξ,t) - Q_x_2(ξ,t), and we deduce

    P_x_1(t,ξ) - P_x_2(t,ξ)- Q_x_1(t,ξ) + Q_x_2(t,ξ)_L_t^∞L_ξ^∞≤  C_1 x_1(t)-x_2(t)_L_t^∞.

with C_1 = 4 max{∂_ξQ_L_ξ^∞,∂_ξP_L_ξ^∞,2/a_0 ℓ∂_ξω_2_L_ξ^∞}. 

We now deduce from (<ref>)

    F[x_1(·)]-F[x_2(·)]_L_t^∞≤( C_1ℓ/2η∂_ξΔ W_L_ξ^∞+ k/η) x_1(t)-x_2(t)_L_t^∞,

and obtain

    ψ[x_1(·)]-ψ[x_2(·)]_L_t^∞≤ C_2 T x_1-x_2_L_t^∞,

with C_2 = C_1 ℓ/2η∂_ξΔ W_L_ξ^∞ +k/η.

Taking T = 1/ 2C_2, this proves that ψ is a contraction, as claimed. Then, there exists a unique fixed point x̃ (·) ∈ C^0([0,T]) which satisfies (<ref>). 


Step 2: global solutions.
The previous construction can be extended as long as P_x and Q_x remain bounded in C^1([0,ℓ]) with respect to ξ and in L^∞(ℝ_+) with respect to t, as shown by the formulas for C_1 and C_2. But from (<ref>, <ref>) we have

    ‖ P_x(·, t) ‖_C_ξ^1≤ e^-a_0t‖ P ‖_C_ξ^1 + (1-e^-a_0t)1/a_0 ℓ‖ω_2‖_C_ξ^1,

which shows that 

    ‖ P_x ‖_C_ξ^1, L_t^∞≤max( ‖ P ‖_C_ξ^1,1/a_0 ℓ‖ω_2‖_C_ξ^1),

and the same goes for Q_x.
We thus obtain that there is a unique global solution (x,  P_x,   Q_x) to (<ref>, <ref>) for all time t≥ 0.

In fact, x∈ C^1(ℝ_+) and P_x,   Q_x ∈  C_#^1([0,ℓ] ×ℝ_+), with ξ↦ P(ξ,·) and ξ↦ Q(ξ,·) in C_#^1([0,ℓ]), as shown by the following bootstrap argument:
From equations (<ref>, <ref>), we know that P_x and Q_x are continuous in time, which enables us to deduce that F[x(·)] ∈ C^0(ℝ_+). Therefore, ψ[x(·)] ∈ C^1(ℝ_+). But, since x = ψ[x(·)], we deduce that x ∈ C^1(ℝ_+). Re-using equations (<ref>, <ref>), 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 ω_2).



 §.§ Proof of Theorem <ref> (Hopf bifurcation)


The functions ξ↦ P(ξ,t) and ξ↦ Q(ξ,t) are periodic with period l, and they can be then expanded in Fourier series

    {[ P(ξ,t) = p_0(t)/2 + ∑_n>0( p^c_n(t) cos2n πξ/ℓ + p^s_n(t)sin2n πξ/ℓ),;  Q(ξ,t) = q_0(t)/2+ ∑_n>0( q^c_n (t)cos2n πξ/ℓ + q^s_n(t)sin2n πξ/ℓ). ].

We insert this expansion and the one for the transition rates (<ref>) into the system (<ref>). 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:

    ṗ_̇0̇(t) = - a_0(Ω) p_0(t) + a_0(Ω)/ 2 ℓ,  q̇_̇0̇(t) = - a_0(Ω) q_0(t) + a_0(Ω)/ 2ℓ.

For n ≠ 0 we obtain:

    {[  ṗ^̇ċ_̇ṅ(t) + 2 π n /ℓ v(t) p^s_n(t)= - a_0(Ω) p^c_n(t) + a_n(Ω)/ ℓ ,;  ṗ^̇ṡ_̇ṅ(t) - 2 π n /ℓ v(t) p^c_n(t)= - a_0(Ω) p^s_n(t) + b_n(Ω)/ ℓ ,; q̇^̇ċ_̇ṅ(t) - 2 π n /ℓ v(t) q^s_n(t)= - a_0(Ω) q^c_n(t) + a_n(Ω)/ ℓ ,; q̇^̇ṡ_̇ṅ(t) + 2 π n /ℓ v(t) q^c_n(t)= - a_0(Ω) q^s_n(t) + b_n(Ω)/ ℓ , ].

together with the force balance equation

    2 ηẋ(t) + 2 k x(t) + π U (p_1^s-q_1^s) = 0 .

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 (<ref>) gives 
    
    p_0(t)=e^-a_0t(p_0(0)- 1/2l) + 1/2ℓ ,   q_0(t)=e^-a_0t(q_0(0)- 1/2l) + 1/2ℓ ,

and both converge to 1/2ℓ 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 (<ref>) and (<ref>) we obtain a first-order five-dimensional ODE:

    {[   ṗ^̇ċ_̇1̇ = - a_0(Ω) p_1^c +   (ζ x +   λ/2 (p_1^s-q_1^s))p^s_1 + a_1(Ω)/ℓ,;  ṗ^̇ṡ_̇1̇ = - a_0(Ω) p_1^s - (ζ  x +   λ/2 (p_1^s-q_1^s)) p^c_1  + b_1(Ω)/ℓ,;  q̇^̇ċ_̇1̇=  - a_0(Ω) q_1^c - (ζ x +   λ/2 (p_1^s-q_1^s)) q^s_1  + a_1(Ω)/ℓ,; q̇^̇ṡ_̇1̇ = - a_0(Ω) q_1^s + (ζ  x +   λ/2 (p_1^s-q_1^s)) q^c_1  + b_1(Ω)/ℓ,;                                  ẋ = -ℓ/2 π(   ζ  x +   λ/2 (p_1^s-q_1^s)) , ].

where ζ = 2 π k / ηℓ and λ = 2 π^2 U / ηℓ.

We linearize the system around its equilibrium point p_eq(Ω) = (p_1,eq^c(Ω),p_1,eq^s(Ω),q_1,eq^c(Ω),q_1,eq^s(Ω), x_eq(Ω)), where

    p_1,eq^c = q_1,eq^c = a_1(Ω) / a_0(Ω) ℓ ,    p_1,eq^s = q_1,eq^s = b_1(Ω) / a_0(Ω) ℓ ,    x_eq=0 .
    
        p_1,eq^c = q_1,eq^c = a_1^0 /a_0^0 ℓΩ^α_1-α_0,  p_1,eq^s = q_1,eq^s = b_1^0/a_0^0 ℓΩ^α_1-α_0,  x_eq=0 ,

We observe that the Jacobian matrix has five eigenvalues: three of them are equal to 
-a_0(Ω)<0 while the other two are of the form τ(Ω) ±1/2 √(π)√(-2  ζℓ a_0(Ω) +  π τ^2(Ω)) , where τ(Ω) is given by (<ref>).

Since, by hypothesis, there exists a real and positive Ω=Ω_0 such that τ(Ω_0)=0, then -2  ζℓ a_0(Ω_0) +  π τ^2(Ω_0)<0, and we deduce that the pair of complex eigenvalues can be written as τ(Ω) ± i ω(Ω) where  

    ω(Ω) := -1/2 √(π)√(2  ζℓ a_0(Ω) -  π τ^2(Ω)),

and they cross the imaginary axis at Ω=Ω_0.

In the following proposition, we are going to study the non-linear behavior of the vector field solution to (<ref>) by exploiting the center manifold theorem; the orbit structure near the fixed point and Ω_0 is determined by the restriction of the non linear system to the center manifold. In particular, system (<ref>) restricted to the center manifold will show a super-critical Hopf bifurcation near p_eq(Ω_0) and Ω_0.

 
    With the hypothesis of Theorem <ref>, the non-linear system (<ref>) has a supercritical Hopf bifurcation near (p_eq(Ω_0),Ω_0).





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 (<ref>) to the origin. We define the new variables as

    δ p_1^c = p_1^c - p_1,eq^c, δ p_1^s = p_1^s - p_1,eq^s, δ q_1^c = q_1^c - q_1,eq^c,  δ q_1^s = q_1^s - q_1,eq^s.



    

 We then use the following linear and invertible change of variables

    [ r = a_1/b_1δ p_1^c + δ p_1^s,  s = a_1/b_1δ q_1^c +δ q_1^s,   z = 1/2(δ p_1^s + δ q_1^s),   y = 1/2(δ p_1^s - δ q_1^s), ]
 
and shorten the notation by taking X=(r,s,z)^T, and Y=(y,x)^T to transform the system (<ref>) into
















    d/dt([ X; Y ]) = 𝕄(Ω) ([ X; Y ]) + ([ G(X,Y); F(X,Y) ]),

where

    𝕄(Ω)=(
    [ 𝔹(Ω)    0;    0 𝔸(Ω) ])












is a block diagonal matrix. Since we assumed that a_0(Ω)=a_0^0Ω^α_0, a_1(Ω)=a_1^0Ω ^α_1, b_1(Ω)=b_1^0Ω^α_1 for a_1^0 ∈ℝ and α_0,α_1,a_0^0, Ω∈ℝ_+, then for the linear part we have 

    𝔹(Ω)=-a_0^0 Ω^α_0Id_3,3 ,  𝔸(Ω)= (
    [ -a_0^0 Ω^α_0 - λ/ℓa_1^0/a_0^0Ω^α_1-α_0              - a_1^0/a_0^0ζ/ℓΩ^α_1-α_0; ;                                -λℓ/2 π                                -ζℓ/2 π;                                        ]),

and the non linear part is defined as

    G(X,Y)=(ζ   x+λ  y)([ a_1^0/b_1^0 (y+z)+ b_1^0/a_1^0(-r+y+z); a_1^0/b_1^0 (y-z)+ b_1^0/a_1^0 (s+y-z);               -b_1^0/a_1^0 (r-s-2 y)/2 ]) ,

and

    F(X,Y)=(ζ  x+λ  y)([ -b_1^0/a_1^0(r+s-2 z)/2;                       0 ]).


Computation of the center manifold.

The center manifold can then be computed by using standard techniques (see <cit.>, Chapter 20, Section 2). We start by rewriting the system in such a way that Ω_0 is moved to the origin through the change of variable δΩ = Ω - Ω_0. As it is classical, we treat δΩ as a variable of the system. This means that we add the equation δ̇Ω̇=0 to the dynamical system and that the non linear part of the system includes all the products δΩ  r, δΩ  s, δΩ  z etc. 

Since the terms in the matrix 𝕄 are nonlinear in Ω, we expand them as

    (Ω_0 + δΩ)^α = Ω_0^α + c(α) δΩ  + O(δΩ ^2),

where c(α) = αΩ_0^α-1.

We insert the expansion (<ref>) into the system (<ref>), neglecting terms of order two in δΩ and getting

    d/dt([  X;  Y; δΩ ]) = 𝕄(Ω_0)([ X; Y; 0 ]) + ([ g(X,Y,δΩ); f(X,Y,δΩ);         0 ]),

with 
    g(X,Y,δΩ)= -a_0^0c(α_0)δΩ X +G(X,Y)


    


    ([                           - a_0^0c(α_0) δΩ  r + (ζ  x+λ  y) ( (a_1^0)^2 (y+z)+(b_1^0)^2 (-r+y+z)/a_1^0b_1^0);                              -a_0^0 c(α_0) δΩ  s +(ζ  x+λ  y) ((a_1^0)^2 (y-z)+(b_1^0)^2 (s+y-z)/a_1^0b_1^0);                                                      -a_0^0c(α_0) δΩ  z -b_1^0 (r-s-2 y)/2 a_1^0 (ζ  x+λ  y); -a_0^0 c(α_0) δΩ  y- λ/ℓa_1^0/a_0^0c(α_1-α_0) δΩ  y -ζ/ℓ c(α_1-α_0) δΩ  x -b_1^0(r+s-2 z)/2 a_1^0(ζ  x+λ  y);                                                                                                            0 ]).


We call 

    𝔹_0:=(
    [ -a_0^0 Ω_0^α_0              0              0;              0 -a_0^0 Ω_0^α_0              0;              0              0 -a_0^0 Ω_0^α_0;                ]),  𝔸_0:= (
    [ -a_0^0 Ω_0^α_0 - λ/ℓa_1^0/a_0^0Ω_0^α_1-α_0                - a_1^0/a_0^0ζ/ℓΩ_0^α_1-α_0; ;                                    -λℓ/2 π                                    -ζℓ/2 π;                                            ]),

the two sub-matrices of 𝕄.
The non linear part is split into 

    g(r,s,z,   y,x,δΩ):=
       ([ - a_0^0c(α_0) δΩ  r + (ζ  x+λ  y) ( (a_1^0)^2 (y+z)+(b_1^0)^2 (-r+y+z)/a_1^0b_1^0);   -a_0^0 c(α_0) δΩ   s +(ζ  x+λ  y) ((a_1^0)^2 (y-z)+(b_1^0)^2 (s+y-z)/a_1^0b_1^0);                          -a_0^0 c(α_0) δΩ   z -b_1^0 (r-s-2 y)/2 a_1^0 (ζ  x+λ  y);                                                                                    ])

and

and

    f(X,Y,δΩ) = -([ (a_0^0 c(α_0) y+λ/ℓa_1^0/a_0^0 c(α_1-α_0))δΩ  y +ζ/ℓ c(α_1-α_0)δΩ  x; ;                                                                    0 ])+ F(X,Y).


We can define a center manifold as 

    W^c(0) = {(X,Y,δΩ): X= 𝐡(Y,δΩ),
        |Y|<, |δΩ |<, 𝐡(0,0)=0, ∇𝐡(0,0)=0},

for  and  sufficiently small and 𝐡=(h_1,h_2,h_3) smooth enough.

In order for 𝐡 to be the center manifold for the system (<ref>), the graph of 𝐡(Y,δΩ) has to be invariant under the dynamics generated by (<ref>). Hence, by plugging 𝐡 into the system and compute its derivative, we get the following quasilinear differential equation 

    ∇_Y𝐡·(𝔸(Ω_0)Y+ f(𝐡,Y,δΩ))=𝔹(Ω_0)𝐡 + g(𝐡,Y,δΩ).


The solution 𝐡(Y,δΩ) of (<ref>) 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 

    h_i(Y,δΩ):= a_i1y^2+a_i2 y x+a_i3 y δΩ  +a_i4 x^2+a_i5 x δΩ  +a_i6δΩ  ^2 + O(Y^3,δΩ ^3).

To completely determine (locally) the center manifold, we have to compute the coefficients a_ij knowing that the functions defined in (<ref>) must solve (<ref>). For the detailed computations, we refer the reader to Appendix <ref>.

Restricted to the center manifold, the original system (<ref>) has the following form:

    d/dtY = 𝔸(Ω) Y + f(𝐡(Y),Y)+O(Y^3,δΩ ^3) .


A complete formula for f(𝐡(Y),Y)) in terms of Y=(y,x)^T follows from our computations and reads

    f(   𝐡(Y),Y)=
       ([ 2 π a_0^0ℓ (λ  y+ζ  x) (πa_1^0Ω_0^α_1-α_0(λ ^2 y^2+ζ ^2 x^2)+ζ a_0^0ℓ (π a_0^0Ω_0^α_0-ζℓ) y x)/(2 πa_1^0λΩ_0^α_1-α_0+ζ a_0^0ℓ^2) (2 πa_1^0λΩ_0^α_1-α_0+a_0^0 ℓ (π a_0^0Ω_0^α_0-ζℓ));                                                                                                                                                                                   0 ]).
 

Observe that τ(Ω)  ± i ω(Ω) is the pair of conjugated eigenvalues of 𝔸(Ω) 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 (<ref>) 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 𝔸(Ω) is transformed to its real Jordan form. Namely, we define the transformation matrix

    ℙ=[ P_11      P_12;    1         0 ],

where 

    P_11=π/λℓ(a_1^0λΩ_0^α_1-α_0/a_0^0 ℓ+a_0^0 Ω_0^α_0)-ζ/2 λ
 and 
    P_12 = -1/λℓ√(πζℓ(a_0^0 Ω_0^α_0-a_1^0λΩ_0^α_1-α_0/a_0^0 ℓ)- (πa_1^0λΩ_0^α_1-α_0/a_0^0 ℓ+π a_0^0Ω_0^α_0)^2-ζ ^2 ℓ^2/4).

The matrix ℙ defines the new coordinates Y= ℙỸ thanks to which equation (<ref>) can be expressed in its normal form

    d/dtỸ = [  τ(Ω) -ω(Ω);  ω(Ω)  τ(Ω) ]Ỹ + ℙ^-1f(𝐡(ℙỸ),ℙỸ)+O(Y^3,δΩ ^3).


To compute key properties of the system, such as the amplitude of the limit cycles, we change coordinates to the polar ones Ỹ^T= ρ (sin(θ),cos(θ)), and get

    {[     ρ̇(t) = τ'(Ω_0) δΩρ(t)+ τ̃ρ^3(t) + O(δΩ^2ρ,δΩρ^3),; θ̇(t) =ω(Ω_0) + ω'(Ω_0) δΩ +  ω̃ρ^2(t) +O(δΩ^2,δΩρ^2). ].,

that is the normal form of (<ref>) in polar coordinates around Ω_0.

The non linear part ℙ^-1f(𝐡(ℙỸ),ℙỸ) determines the constants τ̃ and ω̃. The first one is involved in the expression for the limit cycle amplitude, and we compute it using a well know formula (see <cit.>, Chapter 20, Section 2). We obtain

    τ̃ = - 3 πζ/4 ℓ(π a_0^0Ω_0^α_0 + ℓζ/π a_0^0Ω_0^α_0 + 2 ℓζ)



    Additionally, it is straightforward to compute

    τ'(Ω_0) = -1/2(a_1^0λα_1-α_0Ω_0^α_1-α_0 -1/a_0^0 l+a_0^0).




Since τ̃ is negative and

    τ'(Ω_0) = -1/2(a_1^0λα_1-α_0Ω ^α_1-α_0-1/a_0^0 ℓ + α_0 a_0^0 Ω ^α_0-1),

is positive by hypothesis, the reduced system (<ref>), and hence, the whole system (<ref>), shows a supercritical Hopf bifurcation near the bifurcation parameter Ω_0.

In particular, for Ω sufficiently near and greater then Ω_0 there exists an asymptotically stable periodic orbit with radius as in (<ref>).

    

This finishes the proof of Proposition <ref>


Step 3. Higher order terms. Lastly, we are going to prove that, after large times, the solutions of (<ref>) 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 (<ref>) 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.

   For n > 1 consider the infinite system of four equations:

    {[ ṗ^̇ċ_̇ṅ(t) + 2 π n/ℓẋ p^s_n(t)= - a_0(Ω) p^c_n(t) + a_n(Ω)/ ℓ,; ṗ^̇ṡ_̇ṅ(t) - 2 π n/ℓẋ p^c_n(t)= - a_0(Ω) p^s_n(t) + b_n(Ω)/ ℓ. ].

and

    {[ q̇^c_n(t) - 2 π n/ℓẋ q^s_n(t)= - a_0(Ω) q^c_n(t) + a_n(Ω)/ ℓ,; q̇^s_n(t) + 2 π n/ℓẋ q^c_n(t)= - a_0(Ω) q^s_n(t) + b_n(Ω)/ ℓ. ].

Suppose that, for t ≥ 0, the function t ↦ x(t) is periodic of period T. Then, the solutions to the systems (<ref>) and (<ref>) 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. 


We introduce z_n = p_n^c + ip_n^s.
The first two equations of (<ref>) then become: 

    ż_n + (a_0(Ω) -2inπ/ℓẋ)z_n = c_n

where c_n = a_n(Ω) + ib_n(Ω)/ℓ is constant. We thus deduce the following expression for z_n(t)

    z_n(t)  =  c_n ∫_0^te^-∫_u^t(a_0(Ω) - iv̅(s))dsdu + z_n(0)e^-∫_0^t (a_0(Ω)-iv̅(s))ds ,

where v̅ = 2nπ/ℓẋ.

Let us first notice that the second term, z_n(0)e^-∫_0^t (a_0(Ω)-iv̅(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 (<ref>)), 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 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 
(<ref>) converges to 0, and the result will follow by studying only the first term.

We thus denote by z̃_n(t)= ∫_0^te^-∫_u^t(a_0 - iv̅(s))dsdu, and compute

    [                                                         z̃_n(T+t)- z̃_n(t)   =  ∫_0^t+Te^-∫_u^t+T(a_0 - iv̅(s))dsdu - ∫_0^te^-∫_u^t(a_0 - iv̅(s))dsdu;                                                                            =  ∫_0^Te^-∫_u^t+T(a_0 - iv̅(s))dsdu + ∫_T^t+Te^-∫_u^t+T(a_0 - iv̅(s))dsdu;                                         -∫_0^te^-∫_u^t(a_0 - iv̅(s))dsdu . ]

Using the fact that v̅ is T-periodic, the last two terms cancel and we obtain

    z̃_n(T+t)- z̃_n(t)   =   ∫_0^Te^-∫_u^t+T(a_0 - iv̅(s))dsdu
       =    e^-(t+T)a_0∫_0^Te^ua_0 + ∫_u^t+Tiv̅(s)dsdu

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^c_n,q^s_n) solution to (<ref>).


Finally, if we have initial conditions that are close to the equilibrium point of system (<ref>), and an amount of ATP Ω, which is close to Ω_0, we know exactly how the solution of the system evolves in time. The first order coefficients of (<ref>) go quickly to the constant term 1/l, 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 (<ref>), shows a supercritical Hopf-bifurcation in time, close enough to the parameter Ω_0.


    We remark that, up until the computation of the center manifold, the coefficients a_0, a_1 and b_1 may remain general: there is no need for the power laws to obtain three negative eigenvalues and two complex eigenvalues in system (<ref>). To facilitate the application of center manifold techniques, we have opted to explicitly express the dependence of a_0, a_1, and b_1 on Ω. The selection of power laws is somewhat intuitive, as it extends previous findings (<cit.>,<cit.>) and proves to be advantageous for our analysis.




    
    
  * add numerics and cite: https://link.springer.com/article/10.1140/epje/i2011-11060-5, <cit.>


It is worth mentioning the role of b_1 in the dynamics of the two-row system, since we introduce it to point out that x_eq for the two-row model is always zero, while for the one-row model this is not the case. 

    By looking at system (<ref>), we can observe that the coefficient d appears only in the non linear part ℕ𝕃, implying that the effect of b_1 can be appreciated only with the complete non linear simulations.

First notice that b_1^0 does not appear in 𝕄. Moreover, after the restriction to the center manifold, the non linear part f(𝐡(Y),Y) is independent on b_1^0 too, at least up to third order expansion. Since ẋ(t) only depends on Y, the effect of b_1(Ω) on the displacement solution and its velocity is negligible.

We also remark how the two systems (<ref>) and (<ref>) depends on x(t) and therefore won't depend on b_1(Ω), but rather on a_0(Ω), on a_n( Ω) with n>0 odd, and on b_n(Ω) with n>1 odd. On the other hand, we expect the probabilities P and Q to be sensitive to b_1(Ω), through their first order coefficients. 




We point out that both Theorem <ref> and <ref> can be applied to system (<ref>), leading to the existence of an oscillatory solution for the one-row model, and completing the work done in <cit.>. 

The Hopf bifurcation theorem gives a formal proof of the results carried out in <cit.>, together with new insights on the non linear dynamics which have been often studied  exploiting the non linear system theory, considering the displacement x(t) as input and the external force f_ext(t) as output.




 §.§ Numerical method

We now describe the numerical scheme used throughout the remaining of the paper.

It is used here to compute the solution to the two-row model (<ref>), but 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 Δ x such that ℓ = J Δ x for J ∈ℕ^*, and Δ t the time step. We define P_j^n := P(j Δ x, n Δ t) and Q_j^n := Q(j Δ x, n Δ t) for n≥ 1 and 1 ≤ j ≤ J. Being ℓ-periodic, we may extend P_j^n for all j ∈ℤ by setting P_j+J^n=P_j^n for all j ∈ℤ and all n ∈ℕ; the same can be done for Q_j^n. At a fixed time step t= n Δ t, n≥ 1, the velocity v^n of the central tubule pair being known,  we compute, for j ∈ℤ

    {[                                                                                                                                                                                            S^n_j = -a_0 P^n_j + ω_2(jΔ x)/ℓ,                                                                                                                                                                                                                             ; ;                                                                                                                                                                                            T^n_j = -a_0 Q^n_j + ω_2(jΔ x)/ℓ,                                                                                                                                                                                                                             ; ; P^n+1_j = {[ ( 1 - v^nΔ t/Δ x) + v^nΔ t/Δ xP^n_j-1 + Δ t S^n_j,                                        if  v^n >0,; ( 1 + v^nΔ t/Δ x) - v^nΔ t/Δ xP^n_j+1 + Δ t S^n_j,                                        if  v^n <0, ].; ; Q^n+1_j = {[ ( 1 + v^nΔ t/Δ x) - v^nΔ t/Δ xQ^n_j-1 + Δ t T^n_j,                                        if  v^n >0,; ( 1 - v^nΔ t/Δ x) + v^nΔ t/Δ xQ^n_j+1 + Δ t T^n_j,                                        if  v^n <0, ]. ].

and (x^n+1,  v^n+1) are

    {[ η v^n+1 = Δ x ∑_j=1^J ((P^n+1_j - Q^n+1_j)∂_ξΔ W(jΔ x)) - kx^n,;                                         x^n+1 = x^n + v^n+1Δ t. ].

Notice that, in (<ref>), 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 | v^n Δ t/Δ x| <1.
 



  §.§.§ Conditions on α_0 and α_1


In order to implement the numerical simulations, we have to characterize ω_1 and ω_2 by fixing their coefficients. We first set the values for a_0^0 and Ω_0. For simplicity, we impose that a_n=b_n=0 for n>1 since, as observed in Theorem <ref>, they have no significant influence over the system's oscillations. We now want to find some a_1^0,  α_0,  α_1, that verify  τ(Ω_0)=0 from equation (<ref>) and τ'(Ω_0)>0. Starting with α_0 and α_1, we search for sufficient conditions such that τ'(Ω_0)>0. This inequality reads

    1/2a_0^0α_0 Ω_0^α_0-1 - (α_0 -α_1)λ a_1^0/2a_0^0ℓΩ_0^α_1-α_0-1 < 0,

Substituting a_1^0 into it we get 

    a_0^0(2 α_0-α_1) Ω_0^α_0 +(α_0 -α_1) ℓζ/2 π < 0.

The condition (<ref>) is fulfilled with different choices for α_0, α_1. A particularly simple solution that we will use hereafter is α_0 = 1/2 and α_1 = 1.
Then, to recover a_1^0=a_1^0(Ω_0), we solve τ(Ω_0)=0, from equation (<ref>), which leads to 

    a_1^0= - a_0^0 ℓ/2 πλΩ_0^α_0-α_1(2a_0^0Ω_0^α_0π + ℓζ).

Finally, using the parameter values presented in Table 1, we compute a_1^0 at the instability and get a_1^0(Ω_0) = -56.4588   s^-1.



    

  * Starting from a non zero shift, explain tug of war theory <cit.>
page 10 of Sartori's PhD: <https://www.pks.mpg.de/fileadmin/user_upload/MPIPKS/group_pages/BiologicalPhysics/dissertations/PabloSartori2015.pdf>

  * Conclude on inhibition and curvature regulation





 §.§ Numerical simulations


Motivated by the previous computations, all simulations are carried out defining Δ W(ξ) = U cos (2 πξ / ℓ), a_0(Ω) =a_0^0 √(Ω), a_1(Ω)=b_1(Ω) =a_1^0Ω and a_n(Ω)=b_n(Ω) = 0, if n>1, with α_0 = 1/2 and α_1=1. It is fundamental to notice that the choice to define b_1 ≠ 0, makes the transition rates not-symmetric, as in <cit.>.


All parameters values used for all simulations are shown in Table 1, matching F. Jülicher's <cit.>. We have taken for k_BT the value at 36^∘C (human body temperature). 



As shown in Figure <ref>, the behavior of the two-row model with respect to the ATP concentration is comparable with the one for the one-row model.
When the ATP concentration Ω is lower than Ω_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.



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 <ref>(b), where the displacement takes place around zero. In this case, the equilibrium position does not create any curvature, meaning it has the potential to bend equally in both directions. 

The simulations for the one-row model where carried out using the same upwind numerical scheme. In this case there is only one probability density P(ξ,t). It is linked to the only moving filament, whose displacement is defined as x(t). We used the same notation as the literature, see <cit.> for example.


  



  §.§.§ Illustrating Theorem <ref>

In the following section, we are going to test the results of Theorem <ref> against numerical simulations of the original partial differential equations system (<ref>) and its first order approximation as the ordinary differential equations system (<ref>). 

 All constants are chosen as described in Table 1, and in Section <ref>. The parameter δ indicates the relative distance from the instability Ω_0; therefore, the point where the simulations is performed is Ω_0(1+δ).



Motivated by the previous computations, we numerically solve the ODE (<ref>) and the PDE (<ref>) for Ω = Ω_0(1+δ) with δ > 0, focusing on the amplitude of oscillations for the displacement.

To determine the numerical amplitude for both the ODE and the PDE, we calculate the difference between the maximum and minimum values of the solution after some time t̅, when the numerical solution has reached its limit cycle. This difference is then divided by two.

We then compare the numerical amplitude with the theoretical amplitude defined by (<ref>). For our choice of parameters and neglecting higher-order terms, the theoretical amplitude simplifies to:

    ρ(Ω)= √(δℓ^2/6 π^2(π a_0^0Ω_0^1/2 + 2 ℓζ/π a_0^0Ω_0^1/2 + ℓζ)).





In Figure <ref>, simulations only start after δ=0.05. This is related to the fact that the closer we are to Ω_0, the slower the solution enters into its limit cycle. For this reason, we should have taken a time t̅=t(δ) after which evaluating the amplitude of x(t). In this case, when δ goes to zero, t̅ (δ)→∞. In the simulations this was not reasonable, so we chose a time t̅ independently from δ, but large enough so that the solution enters in its steady-state regime for all δs. 

In Figure <ref>(a), we observe the amplitude of x(t) computed with the PDE, the ODE, and the analytical result (<ref>). As expected, the last loses its prediction capability the more we go away from the instability; and subsequently, the relative error between the truncated formula (<ref>) and the one measured from the PDE increases, as we can see in Figure <ref>(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 Ω_0. 











    




    
    

    



    
   
    
    







  §.§.§ Comment on the amplitude


The expansion (<ref>) of the amplitude ρ can also be adapted to the ODE formulation of the one-row model. In Figure <ref>, 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. 







§ TOWARDS MODELING THE AXONEME: N-LAYER MODEL WITH FIXED EXTREMITIES

Starting from two rows of molecular motors, we extend the model from section <ref> to N ≥ 2 rows of molecular motors, without any explicit form of inhibition in the system. 

The system is then composed of N+1 microtubule 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. 

To take the bridge of an axoneme into account (as shown in Figure <ref>) in our N-layer model, filament pairs 0 and N are considered to be the same pair. They have no shift between each other, as shown in the split view Figure <ref>.
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 ℓ.



 §.§ Full mathematical model


For i ∈{ 0, …, N }, we denote by X_i the horizontal shift of the i-th tubule, as shown in Figure <ref>.
We measure displacement with respect to the 0-th filament: X_0 = X_N = 0. For 1 ≤ i ≤ N, we introduce Δ_i = X_i - X_i-1, the relative displacement between filaments i and i-1. It follows that ∑_i=1^NΔ_i(t) =0. The shifting speed of the i-th filament is defined by V_i = Ẋ_i - Ẋ_i-1 = Δ̇_i, and we have

    ∑_i=1^N V_i = ∑_i=1^NΔ̇_i(t)=0.



Once again, the variable ξ∈ [0, ℓ] represents the local variable along one pair, as in the two-row model. We denote by Q_i(ξ, t) the density of molecular motors at position ξ and time t attached to the i-th filament and who are in state 1, walking on the (i-1)-th filament. 


As for the two-row model, in section <ref>, (Q_i,  Δ_i) are solution to,  for ξ∈ [0, ℓ] and t>0

    {∂_t Q_j(ξ,t)      = -(ω_1(ξ - Δ_j(t), t) + ω_2(ξ - Δ_j(t))) Q_j(ξ, t) 
        + ω_2(ξ - Δ_j(t)))/ℓ,   for j ∈{ 1, …, N },
    η (V_ i(t)-V_i+1(t))     =  ∫_0^ℓ (Q_i(ξ,t)∂_ξΔ W(ξ) - Q_i+1(ξ,t)∂_ξΔ W(ξ - Δ_i(t))) dξ
          - k(Δ_i(t) - Δ_i+1(t)),   for i ∈{ 1, …, N-1 }.
      .
  
Where the first equation is a transport equation for the motors probability density. Shifting speed is measured by looking at the current filament, as well as the one above. We thus obtain the second equation in system (<ref>) by writing the force balance on the i-th filament with the motors and springs around it. We then have N-1 force-balance equations that determine the filaments' motion, in which we impose the external forces to be zero.

Following again section <ref>, for j ∈{ 1, …, N }, t>0 and ξ∈ [0, ℓ], we define P_j as P_j(ξ,t) = Q_j(ξ + Δ_j(t),t).
We then observe that ∂_t P_j(ξ + Δ_j, t ) = V_j(t) ∂_ξ Q_j(ξ,t) + ∂_t Q_j(ξ,t), and we finally obtain the system at any t>0 and  ξ∈ [0, ℓ]

    {∂_t P_j + V_j ∂_ξ P_j = -(ω_1 + ω_2) P_j + ω_2/ℓ,       for j ∈{ 1, …, N },
    η (V_ i-V_i+1) = ∫_0^ℓ (P_i - P_i+1) ∂_ξΔ W - k(Δ_i - Δ_i+1),       for i ∈{ 1, …, N-1 },
    ∑_i=1^NΔ̇_i(t)=0   and Ẋ_n = 0.
      .

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

    {[                       V_i = 0,;                       X_i = 0,; P_j = P^0 = ω_2/ℓ (ω_1 + ω_2). ].





 §.§ Theoretical considerations
 
In this situation we want to investigate the existence of a Hopf bifurcation in time, as it was done with the two-row model. To do so, we consider a linearized system arising from the N-row model (<ref>), 
in order to have a general idea of the system's dynamics, as in (<ref>). We expand in Fourier series P_1, …, P_N and treat order by order the Fourier coefficients of the probabilities. 
For the first order coefficients, we obtain a 3N-1 dimensional ODE system with unknowns: 
    p_j^c,s,   j=1,…,N,     X_k,  k=1,…,N-1.

We then perform some linear change of variables, similarly to (<ref>), 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(Ω). The rest of the linearized system is defined by the following equations 

    ([ q̇_k;  ẇ_k;      ])=
        [ - 2 π^2/ℓη U p_e - a_0         - 2 π/ℓη k p_e;               - π U /η                - k / η;                        ]([ q_k; w_k;     ])

for k=1,…,N-1, where  p_e = a_1/(a_0 ℓ), q_k = p_k^s - p_k+1^s for k=1,…,N-1, and 

    ([   w_1;     ⋮; w_N-1 ]) =
         [  2 -1         ; -1  2 -1      ;        ⋱      ;       -1  2 -1;          -1  2 ]([   X_1;     ⋮; X_N-1 ]).


It follows that the Jacobian has N+1 real and negative eigenvalues -a_0(Ω), and N-1 identical pairs of complex and conjugated ones μ_k(Ω) for k=1,…, N-1 which come from (<ref>). We observe that μ_k(Ω)=τ(Ω) + i ω(Ω), with τ and ω defined as in Theorem <ref> by equations (<ref>) and (<ref>). 
Then, if there exists Ω=Ω_0 such that τ(Ω_0)=0 and τ'(Ω_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 Ω_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. 









 §.§ 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 <ref>. We use the same parameter values as for the previous systems, as specified in Table 1.
The initial conditions are such that P_i(t=0) = P^0 for all i, and we destabilize the X_is randomly in [-x_0, x_0], where x^0 = 0.01   nm.



As we are looking for a setting in which the axoneme generates curvature along the flagellum, we are looking for a setting in which the system alternates between shifting directions in an organized manner. We initially start by destabilizing each half of the system in opposite directions:

    {[              X_i = ix^0,      if  1 ≤ i ≤⌊ N/2 ⌋,;          X_i = (N-i)x^0, if ⌊ N/2 ⌋ +1 ≤ i ≤ N-1,;                 X_N = 0.                          ].

    where x^0 = 0.01   nm. 



Using the same notations as before, 
i.e. Ω = (1+δ) Ω_0, we measure distance the to the bifurcation point Ω_0 using δ.






Figure <ref>(a) shows the absence of oscillations before instability (Ω = 0.9Ω_0), 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 <ref>(b), we look at the system past the instability point, for Ω = 1.1 Ω_0. 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 <ref>, all Δ_is have the same amplitude and oscillation frequency, but their phase difference varies depending on initial conditions. In the particular case shown in Figure <ref>(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) 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  <ref>(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.







Then, using another set of initial conditions, where exactly one half of the system is pushed in one direction, and the other half in the other one, one can observe that the phase difference is maintained through time and each layer in one half of the axoneme oscillates in sync with the other ones from the same half, as shown in Figure <ref>. These initial conditions are denoted as "middle" in the next paragraphs, and implemented as

    {[              X_i = ix_0,      if  1 ≤ i ≤⌊ N/2 ⌋,;          X_i = (N-i)x_0, if ⌊ N/2 ⌋ +1 ≤ i ≤ N-1,;                 X_N = 0.                          ].



The synchronization shown in Figure <ref>(c) was predicted by Howard et al. <cit.>, 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 <cit.>. The Kuramoto model is fairly different from ours, as each oscillator has its own velocity, whereas all of our layers have one common velocity value. Our system of equations cannot thus be easily reduced to a Kuramoto model, but the behaviors of the two systems in terms of synchronization and dependence on initial conditions are very much alike. Some first steps, only modelling individual motors along a single row, have been presented in <cit.>.




















  §.§.§ Is there an optimal number of layers?


Figure <ref> shows the oscillation amplitude depending on the number of layers. For both relative and absolute shift, the maximum amplitude is reach for N between 8 and 12. In particular, the relative shift is clearly maximal for N=8, i.e. for a 9-axoneme model. Relative shift is directly responsible for curvature, showing that the number of motor rows in an axoneme is the optimal value to maximize bending.



either we plot a better graph or we can't really conclude anything from it
    






  §.§.§ Influence of b_1 on the oscillations -  Note: doesn't work with middle initial conditions


In this subsection, the N-layer model is simulated with a different transition rate. We impose b_1(Ω) = 0, i.e. ω_2(ξ; Ω) = a_0(Ω)/2 + a_1(Ω) cos2 πξ/ℓ. 

P0 changes with b1. test b1 =  alpha a1
explore initial conditions
amplitude graph?



Figure <ref> shows tubule oscillation with this new ω_2, after the instability point, with all other parameters unchanged. As expected from the amplitude estimation from the proof of Theorem <ref>, the coefficient b_1 has no influence over oscillation amplitude. Oscillation frequency remains unchanged as well. However, every tubule is slightly out of phase with respect to its neighbor. This time, all oscillations are centered in zero but there are three groups of Δs oscillating together instead of three. 

We would like to find some parameters with a physical meaning, for which we can clearly identify two halves of the N-axoneme oscillating in opposite directions, with a π/2 phase lag and centered in 0.





 §.§ Changing potentials

In this section, we highlight the influence of our choice for the potentials W_i and transition rates ω_i. To do so, let us consider once again the one-row model.
In Jülicher's first models <cit.>, transitions rates are directly derived from chemical considerations, by looking at reaction speed in ATP hydrolysis. In the case of a two headed motor, the reactions can be written as:

    {[ M_1 + ATP k_1k_-1⇌ M_2 + ADP +P; M_1 + ADP + Pk_2k_-2⇌ M_2 + ATP ].


Which leads to the following reaction speed:


    d[M_1]/dt    =  -k_1 [M_1] + k_-1 [M_2] -k_2 [M_1] + k_-2[M_2] 
       =  -δ(x) e^W_1(x) + μ_ATP/kT [M_1] + δ(x) e^W_2(x) + μ_ADP + μ_P/kT [M_2] 
        -δ(x+ℓ/2) e^W_1(x) + μ_ADP + μ_P/kT [M_1] + δ(x+ℓ/2) e^W_2(x) + μ_ATP/kT [M_2] 
       =  -ω_1[M_1]+ω_2 [M_2],

where

    {[ ω_1(ξ)=  [δ (ξ)e^μ_ATP/kT+δ(ξ+ℓ/2)] e^W_1(ξ) / k_B T;  ω_2(ξ)=  [δ(ξ)+δ(ξ+ℓ/2)e^μ_ATP/kT] e^W_2(ξ) / k_B T ].

which directly and explicitly depend on both ATP and physical potentials.

mu adp and mu p? go back to original computation

Here, δ is a Gaussian centered around the minimum of W_1, which helps favor transition at a specific point in the periodicity cell.
For W_1, we chose an arbitrary periodic potential that was highly nonlinear in ξ, so as to ensure it was not symmetric around ℓ/2:


    {[ W_1(ξ) = U cos(2 π/ℓ (0.6 ξ^3 + 0.4ξ));                    W_2(ξ) = W_1(ξ+ℓ/2);                                        ].


We also introduce, as previously done by <cit.>, the local deviation from the chemical equilibrium of the system Ω(ξ), to model the dependence of the system on ATP concentration, depending on our parameter Ω:


    Ω(ξ) = ω_1(ξ)/ω_2(ξ) - exp[(W_1(ξ)-W_2(ξ))/ k_B T],


where k_B is the Boltzmann constant and T is the temperature of the system. 
We use the amplitude of this function, defining it as Ω := sup_ξ∈ [0,ℓ]| Ω(ξ)|.

In practice, we will only work with Δ W and its usual assumptions, as well as the ones we originally made on ω_2 in the introduction. Moreover, using both our previous conditions and our chemical expressions (<ref>), we can arrive to a series of equations on both Ω and α, which could eventually prove that there exists a solution that combines all our assumptions. The system is explained in further detail in <ref>.



We can thus clearly observe in Figure <ref> our choice has had an influence on the solution, as displacement x becomes highly nonlinear, and of much larger periodicity in time.





§ CONCLUSION AND OUTLOOK





In this paper, we presented a new model for the axoneme, the cytoskeleton of cilia and flagella, focusing on the presence of spontaneous oscillations and potential synchronization of microtubule doublets under varying ATP concentrations. We called this model N-row model, where N represents the number of rows of motors walking between two microtubule doublets that are arranged in a circular configuration. Molecular motors are put into motion thanks to a chemical reaction involving ATP, thus inducing microtubule sliding. This specific structural arrangement results in all microtubule doublets having a zero equilibrium position. The dynamic is governed by N transport equations, coming from a stochastic two-state model, and N-1 force balance equations, collectively forming a PDE system with 2N-1 equations. We analyzed this model using both theoretical methods and numerical simulations.
Our first studies focused on the N=2 configuration, revealing the existence and uniqueness of the displacement solution. Moreover, under sufficient chemical energy, the middle microtubule doublet exhibits spontaneous oscillations, corresponding to the existence of a super-critical Hopf-bifurcation in the dynamical system. Theoretical results are corroborated by numerical simulations. Subsequently, we conduct numerical analysis of the N-row model, focusing particularly on the case N=8. Here too, we observe that the system starts to oscillate once the ATP concentration surpasses a critical threshold. Remarkably, 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.

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 <cit.>, previously solved by choosing symmetric transition rates. Instead, in of this paper, asymmetric 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.

From a more theoretical point of view, both Theorem <ref> (well-posedness) and Theorem <ref> (existence of a Hopf bifurcation) give a complete mathematical description on the two-row model. The same results can be proven for the one-row model as well, formalizing the work done in <cit.>.
With the N-row model with fixed extremes, we set up a framework that is useful to gain insights on the microscopic structure of the axoneme. In particular, we comment on the specific case N=8, which is the closest representation of the real axoneme in terms of number of microtubule doublets. With the aid of the numerical scheme presented in the two-row model section, we observed the onset of oscillations at a critical value of ATP. To perturb the system we tried several random initial conditions. We then noticed that imposing an opposite shift on each half of the model immediately led to robust and realistic patterns, potentially leading to planar flagellar beating patterns, as discussed in <cit.>.






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 <cit.> and <cit.> for the planar case and of <cit.> 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.




§ CENTER MANIFOLD COMPUTATION

In practice, equation (<ref>) is a system of three equations. With the first equation below we aim to find the coefficients a_1j, with j=1,…,6

    (2 a_11 y + a_12  x + a_13δΩ  )    (  -a_0^0 Ω_0^α_0y - λ/ℓa_1^0/a_0^0Ω_0^α_1-α_0y  - a_1^0/a_0^0ζ/ℓΩ_0^α_1-α_0x.
       .-a_0^0 c(α_0)δΩ    y
            - λ/ℓa_1^0/a_0^0 c(α_1-α_0)δΩ   y -ζ/ℓ c(α_1-α_0)δΩ   x .
       - . b_1^0(r+s-2 z)/2 a_1^0(ζ  x+λ  y) ) 
        + (a_12 y  + 2 a_14  x + a_15δΩ  )  (-βℓ/2 π y -αℓ/2 π x) 
        + a_0^0Ω_0^α_0 h_1 +  a_0^0c(α_0)δΩ    h_1
       - (ζ  x+λ  y) ( (a_1^0)^2 (y+z)+(b_1^0)^2 (-r+y+z)/a_1^0b_1^0) = 0.

Since we only need terms x, y and δΩ up to order two, we get rid of all the third order terms, getting

    (2 a_11 y + a_12  x + a_13δΩ  )    (  -a_0^0 Ω_0^α_0y - λ/ℓa_1^0/a_0^0Ω_0^α_1-α_0y  - a_1^0/a_0^0ζ/ℓΩ_0^α_1-α_0x )
        + (a_12 y  + 2 a_14  x + a_15δΩ  )  (-βℓ/2 π y -αℓ/2 π x) 
        + a_0^0 Ω_0^α_0 h_1  - (ζ  x+λ  y)  ( (a_1^0)^2 y+(b_1^0)^2 y/a_1^0b_1^0) = 0.


We then proceed by matching the terms x, y and δΩ with same order, and solve the resulting system of equations for a_1j. Notice we could have kept only the zero order term in the Taylor expansion, since all the higher order terms disappear during this approximation.

We do the same for the second and third equation, solving for the power series of h_2 and h_3.

For i=1,2, we obtain

    a_i1= -2 π ^2  λ ^2a_0^0ℓ((a_1^0)^2+(b_1^0)^2) Ω_0^α_1-α_0/b_1^0(2 πa_1^0λΩ_0^α_1-α_0+ζ a_0^0ℓ^2) (2 πa_1^0λΩ_0^α_1-α_0+a_0^0 ℓ(π a_0^0Ω_0^α_0-ζℓ))
        a_i2= 2 πζ  (a_0^0)^2 ℓ^2 ((a_1^0)^2+(b_1^0)^2) (ζ  l-π a_0^0Ω_0^α_0)/a_1^0b_1^0(2 πa_1^0λΩ_0^α_1-α_0+ζ a_0^0ℓ^2) (2 πa_1^0λΩ_0^α_1-α_0+a_0^0 ℓ(π a_0^0Ω_0^α_0-ζℓ))
       a_i3= 0,    a_i4= -2 π ^2 ζ ^2 a_0^0ℓ((a_1^0)^2+(b_1^0)^2) Ω_0^α_1-α_0/b_1^0(2 πa_1^0λΩ_0^α_1-α_0+ζ a_0^0ℓ^2) (2 πa_1^0λΩ_0^α_1-α_0+a^0_0 ℓ(π a_0^0Ω_0^α_0-ζℓ))
       a_i5= 0,    a_i6= 0.

For i=3 we get

    a_31= -2 π ^2 λ ^2a_0^0 b_1^0ℓΩ_0^α_1-α_0/(2 πa_1^0λΩ_0^α_1-α_0+ζ a_0^0ℓ^2) (2 πa_1^0λΩ_0^α_1-α_0+a_0^0 ℓ(π a_0^0Ω_0^α_0-ζℓ))
        a_32= 2 πζ  c^2b_1^0ℓ^2 (π a_0^0Ω_0^α_0-ζℓ)/a_1^0(2 πa_1^0λΩ_0^α_1-α_0+ζ a_0^0ℓ^2) (2 πa_1^0λΩ_0^α_1-α_0+c ℓ(π a_0^0Ω_0^α_0-ζℓ))
       a_33= 0,    a_34= -2 π ^2 ζ ^2a_0^0 b_1^0ℓΩ_0^α_1-α_0/(2 πa_1^0λΩ_0^α_1-α_0+ζ a_0^0ℓ^2) (2 πa_1^0λΩ_0^α_1-α_0+ a_0^0 ℓ(π a_0^0Ω_0^α_0-ζℓ))
       a_35= 0,    a_36= 0.






10

BAYANI2023
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.

bayly_steady_2016
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.

brokaw_flagellar_1972
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.

brokaw_bending_1999
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.

brokaw_computer_2014
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.

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

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

Costantini_2024
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.

Gaffney2011
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.

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

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

guerin2011dynamical
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.

Howard2022
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.

hu_finite_2018
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.

Julicher1997Modeling
F. Jülicher, A. Ajdari, and J. Prost.
Modeling molecular motors.
Rev. Mod. Phys., 69:1269–1282, 1997.

Julicher1995
F. Jülicher and J. Prost.
Cooperative molecular motors.
Phys. Rev. Lett., 75:2618–2621, 1995.

Julicher1997
F. Jülicher and J. Prost.
Spontaneous oscillations of collective molecular motors.
Phys. Rev. Lett., 78:4510–4513, 1997.

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

lindemann_many_2021
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.

machin_wave_1958
K. E. Machin.
Wave Propagation along Flagella.
J. Exp. Biol., 35(4):796–806, 1958.

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

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

sartori2016twist
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.

Sartori2016
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.

Summers
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.

Taylor1951
G. I. Taylor.
Analysis of the swimming of microscopic organisms.
Proc. R. Soc. Lond. A, 209:447–461, 1951.

Velho2021
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.

Walker2020
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 spermatozoa.
R. Soc. Open Sci., 7(6), 2020.

Wiggins
S. Wiggins.
Springer, Texts in Applied Mathematics.
2003.

woodhams_generation_nodate
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.

Yagi1995
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.



