# Tracking phase synchronization between flagella in the time-frequency domain resolves photophobic response

Lucas Fedelspiel, Françoise Argoul, Antoine Allard, Jorge Arrieta, Marco Polin

► To cite this version:

Lucas Fedelspiel, Françoise Argoul, Antoine Allard, Jorge Arrieta, Marco Polin. Tracking phase synchronization between flagella in the time-frequency domain resolves photophobic response. 2025. hal-05405071

HAL Id: hal-05405071

https://hal.science/hal-05405071v1

Preprint submitted on 8 Dec 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
# Tracking phase synchronization between flagella in the time-frequency domain resolves photophobic response

Lucas FEDERSPIEL, Françoise ARGOUL, Antoine ALLARD*

Univ. Bordeaux, CNRS, LOMA, UMR 5798, F-33400 Talence, France

Jorge ARRIETA

Departmento de Física, Universitat de les Illes Balears, 07122, Palma, Spain

Marco POLIN

Mediterranean Institute for Advanced Studies, IMEDEA, UIB-CSIC, Esporles, 07190, Spain

(Dated: October 27, 2025)

The unicellular microalga *Chlamydomonas reinhardtii* (CR) is well known for its bi-flagellated swimming in response to light stimuli. This work aims to study the resynchronization of CR flagella after a high light intensity stimulus, known as photoshock. The synchronization is estimated thanks to a quantity defined as the Phase Synchronization Index (PSI). The originality of this approach is to perform a time-frequency computation of a complex PSI based on continuous wavelet transform. Thanks to this analysis, we distinguish three swimming stages involving different frequency bands and phase shifts: normal breaststroke, escaping, and resynchronization. This approach also reveals the presence of signal harmonics that set the photoshock response, independently of cell variability. Our results suggest that CR modulates the balance between fundamental and harmonic beating modes, providing a mechanism for robust adaptation to sudden environmental stresses.

## I. INTRODUCTION

Oscillators are systems that exhibit repetitive behavior around a steady state. They are found across a wide range of fields, including mechanics (e.g., pendulums and vibrating strings), optics (e.g., electromagnetic oscillator), or active matter systems both inanimate and living [1]. In biology, rhythmic processes such as neuron spikes [2], cardio-respiratory activity [3], and circadian rhythms [4], display oscillatory dynamics across multiple spatial and temporal scales. While the behavior of one isolated oscillator is generally well described, two or more coupled oscillators exhibit much more complex responses [5, 6]. The concept of synchronization, where oscillators adjust their dynamics to cooperate, was first introduced in physics by Huygens in the 17th century, who observed phase or anti-phase coupling of two pendulum clocks [7].

In living systems, synchronization is essential for robust physiological function: it coordinates gait during locomotion, supports coherent brain activity, and drives fluid transport by ciliary carpets [8–10]. However, many biological oscillators exhibit nonstationary and nonlinear behavior. This means that their natural frequency can change dynamically in response to internal processes or environmental perturbations [11, 12]. In such systems, transient coupling between oscillators appears, making synchronization detection challenging. Addressing this challenge requires time-resolved tools that track not only the instantaneous frequency evolution of each oscillator, but also their relative phase dynamics, which govern synchronization.

Here, we introduce a time-frequency synchronization framework, based on continuous wavelet transform (CWT) capable of capturing transient and frequency-dependent coupling, enabling insight into nonstationary biological oscillators [13–18]. Within this framework, we introduce a time- and frequency-resolved Phase Synchronization Index (PSI), that quantifies transient coupling and reveals how synchronization evolves across multiple time scales. This method not only detects when synchronization is present, but also identifies how characteristic frequency bands shift during external perturbations. We applied this framework to experimental data on flagella synchronization in the green alga *Chlamydomonas reinhardtii* (CR), and complemented these results with a minimal phase-coupling model with time-varying intrinsic frequencies [19, 20]. The later allows us to probe how varying coupling strength can help interpret the observed experimental behaviors.

CR is a model organism to study motility, phototaxis or sensory responses [21]. CR typically swims using a breaststroke-like motion, in which its two flagella beat synchronously [22]. This coordination results from intracellular coupling, which is necessary [23] and sufficient [24] for the wild-type strain CC125. In the case of the CR strain *ptx1*, known to have a weaker intracellular coupling, hydrodynamic coupling impacts resynchronization [24, 25]. The ability of CR to thrive in various conditions makes it an ideal organism for exploring adaptive responses to changes in environmental cues [24, 26–30]. Under strong light stimulation, CR undergoes a photo-shock response, transiently reversing its swimming direction through a switch to a high-frequency, undulatory beating mode [11, 28, 31]. Although steady-state flagellar coordination in CR has been extensively studied, little is known about how synchronization dynamically

<footnote>* francoise.argoul@u-bordeaux.fr and antoine.allard@ubordeaux.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. **Experimental setup.** (a) *C. reinhardtii* cell trapped in a micropipette and bathed with tracer beads. Dashed squares indicate the regions of interest where velocity signals are extracted ($Y = 0, \pm 1$). Schematic lines represent a screenshot of the flagella before photoshock (black, -PS) and after photoshock (gray, +PS). Scale bar: 10 µm. (b-c) Velocity components along the *x*-axis (b) and *y*-axis (c), averaged in front of the cell ($Y = 0$), before photoshock (black) and after photoshock (gray).</figcaption></figure>

breaks and re-emerges following such perturbations.

In this work, we address this gap by studying flagellar resynchronization during photoshock. Since the oscillatory motion of flagella is performed at low Reynolds number, we infer oscillatory dynamics from the surrounding flow field instead of directly tracking flagellar shapes, providing robust estimates of beating frequency and phase [32, 33]. Using our wavelet-based PSI, we resolve the time-dependent synchronization of the two flagella. We identify three distinct dynamical stages: (i) pre-stimulus forward swimming with breaststroke beating around 40 Hz, (ii) a transient high-frequency mode (up to 80 Hz) associated with backward swimming, and (iii) gradual recovery into the original breaststroke rhythm. Beyond synchronization, our analysis uncovers the presence and evolution of harmonic components throughout the alga response to photoshock. Intriguingly, we find that the second harmonic of the breaststroke mode is continuously present, even during photoshock, and becomes dominant in the backward-swimming stage. The relative energy of the fundamental and harmonic modes shifts gradually, suggesting that CR does not simply switch between discrete beating patterns but rather modulates the balance of coexisting oscillatory modes. This continuous interplay between fundamental and higher harmonics may underlie the robustness and adaptability of flagellar coordination in changing environments.

## II. MATERIALS AND METHODS

### A. Experimental setup

*Chlamydomonas reinhardtii* strain CC125 (Chlamydomonas Resource Centre, USA) was cultured in liq-

uid medium under continuous agitation at 120 rpm to prevent sedimentation. A 14 h light/10 h dark cycle was maintained to synchronize cell division. Experiments were performed on a bright-field inverted microscope (Nikon Eclipse TE2000-U) modified with a 780 nm LED (Thorlabs, M780L3) for Köhler illumination, a 60X/1.00W water-immersion objective (Olympus), and a high-speed camera (Edgertronic SC1).

Experimental chambers were prepared with 100 µL of algal suspension at a concentration of $5 \times 10^4$ cells/mL, bathed with 1 µm diameter polystyrene beads at a volume fraction of 0.37%. Chambers were kept under infrared illumination for at least 30 minutes prior to measurements to stabilize the cells. Micropipettes were fabricated from borosilicate capillaries (1B100-6, WPI) using a laser puller (Sutter P-2000; parameters: Heat 996, Pull 90, Vel 13, Time 250, Pressure 560), and subsequently forged with a heated filament (MF200-H3, WPI) to refine the tip. For each experiment, a pipette filled with water was introduced into the chamber. Cells were immobilized by gentle suction, and oriented so that both flagella beat within the imaging plane (Fig. 1(a)), though the *cis* and *trans* flagella could not be reliably distinguished experimentally.

Thirteen independent cells were considered. For each, movies of 25 s duration were recorded at 1000 frames per second, with a spatial resolution of $496 \times 496$ px (1 px = 0.222 µm). Photophobic responses were triggered using a 470 nm LED (Thorlabs, M470L5) mounted laterally to the microscope and directed towards the chamber. The LED delivered 50 ms pulses of high-intensity light (5 V amplitude) every 5 s. Each cell was subjected to 3-5 consecutive photoshocks.

The fluid velocity field was measured using Ghost Particle Velocimetry (GPV) [34], with velocity vectors computed over $32 \times 32$ px interrogation windows (17 px overlap) along a rolling window of four frames. To analyze flagellar dynamics, we extracted velocity signals from three representative regions of interest (Fig. 1(a)). Two windows of $20 \mu m \times 10 \mu m$, located symmetrically on either side of the cell body ($Y = -1$ and $Y = 1$), at a distance of approximately 15 µm to predominantly capture the flow generated by the adjacent flagellum. A third window of $20 \mu m \times 30 \mu m$ was placed in front of the cell ($Y = 0$), providing a measure of the bulk swimming flow. For this central region, the recorded signal reflects the combined contribution of both flagella, without distinction between them. From each window, velocity signals were decomposed along two orthogonal directions: the *x*-axis, aligned with the main body axis of the cell, and the *y*-axis, perpendicular to it. Fig. 1(b-c) shows representative velocity traces from the bulk ($Y = 0$) before (black) and after (gray) a photoshock.
## B. Modeling of coupled oscillators

Harmonic oscillators can generally be described as $x(t) = A(t) \cos(\varphi(t))$ with $\varphi(t) = \omega t + \varphi_0$, $\varphi(t)$ the phase, $A(t)$ the amplitude and $\omega/2\pi$ the oscillator frequency [25]. For simplicity, we assume that the amplitude $A$ is constant and does not depend on the oscillator, and we concentrate on the coupling dynamics of oscillator phases. Note that this assumption implies that either the dynamics of the amplitude is uncoupled to that of the oscillator phase, or the amplitude varies very little and the system is not much dissipative. Such an hypothesis is obviously not fully rigorous in the context of living systems. However, for finite time intervals, dissipation effects could be considered as negligible.

In a very generic way, phase coupling dynamics can be written as an ordinary differential equation (ODE) system:

$$ \left\{ \begin{aligned} \dot{\varphi}_1 &= \omega_1 + C(\varphi_1 - \varphi_2) \\ \dot{\varphi}_2 &= \omega_2 + C(\varphi_2 - \varphi_1) \end{aligned} \right. . \quad (1) $$

The coupling function $C(\cdot)$ is $2\pi$ periodic, and can be written as an harmonic function: $C(\delta(t)) = B \cos(\delta(t) + \text{cst})$, with $\delta(t) = \varphi_1(t) - \varphi_2(t)$. We take $\text{cst} = 0$. The difference of these two equations yields the Adler equation $\dot{\delta} = \Delta\omega + 2B \sin(\delta)$, [19] where $\Delta\omega = \omega_1 - \omega_2$.

If we include independent additive Gaussian white noise terms ($\xi_1$ and $\xi_2$) in the phase dynamics, the ODE system (1) reads:

$$ \left\{ \begin{aligned} \dot{\varphi}_1(t) &= \omega_1 + B \sin(\delta(t)) + \xi_1(t) \\ \dot{\varphi}_2(t) &= \omega_2 - B \sin(\delta(t)) + \xi_2(t) \end{aligned} \right. . \quad (2) $$

The noise terms $\xi_1$ and $\xi_2$ follow $\langle \xi_i(t) \rangle = 0$ and $\langle \xi_i(t) \xi_j(t') \rangle = 2D\delta_{ij}\delta(t-t')$, where $\delta_{ij}$ is the Kronecker delta and $\delta(t-t')$ the Delta function. To model this ODE system (2), we use the Python module `sdeint` [35], with built-in noise amplitudes $G_1 = G_2 = 5$ for Fig. 2 and $G_1 = G_2 = 1$ for Fig. 3.

## C. Time frequency analysis of non-stationary signals

The simplest and most common method for spectral decomposition of signals is the Fourier transform [36]. Although effective for stationary signals, temporal windows must be introduced to restrict the analysis to finite time intervals when the spectral signature of the signal changes over time as described in [13].

Here, we use the Continuous Wavelet Transform (CWT), which provides a time-frequency representation of a signal, i.e. the temporal changes of both its amplitude and its phase [14]. For a real signal $S(t)$, the CWT yields the complex two-parameter function $W(a, b)$:

$$ W(a, b) = \frac{1}{|a|^{1/2}} \int_{-\infty}^{+\infty} S(t) \bar{\psi} \left( \frac{t-b}{a} \right) dt \quad (3) $$

<figure><img src="image_3.png" /><figcaption>FIG. 2. Simulations of two coupled oscillators with Gaussian noise, modeled by Eq. (2). (a) Strong negative coupling ($B/2\pi = -20$ Hz < 0, $|2B| > \Delta\omega$), with $\varphi_1(0) = 0$, $\varphi_2(0) = \pi$: leads to stable in-phase synchronization. (b) Weak negative coupling ($B/2\pi = -3$ Hz < 0, $|2B| < \Delta\omega$), with same initial phases: leads to desynchronization with transient coordination. (c) Strong positive coupling ($B/2\pi = 20$ Hz > 0, $|2B| > \Delta\omega$) with $\varphi_1(0) = \varphi_2(0) = 0$: leads to stable anti-phase synchronization. Left panels: $S_1 = \cos(\varphi_1)$ (gray) and $S_2 = \cos(\varphi_2)$ (black). Right panels: $\delta = \varphi_1 - \varphi_2$. Dashed and solid lines represent anti-phase and in-phase locking, respectively.</figcaption></figure>

with $\psi(t)$ the mother wavelet, $a$ the scaling factor, and $b$ the translation factor. $\bar{\psi}$ stands for the complex conjugate of $\psi$. For convenience, in the remainder of the paper we use the equivalent notation $W(t, f)$, where the scale-frequency conversion is given by $f = f^*/a$, with $f^*$ being the frequency at which the Fourier transform of the mother wavelet, $\Psi(f)$, peaks. Therefore, $W(t, f)$ and $W(a, b)$ represent the same quantity, expressed respectively in the time-frequency and time-scale domains. The Continuous Wavelet Transforms were performed using the python library `pycwt`. Since our simulated and experimental signals are close to sinusoidal, we use a Morlet mother wavelet, with a wave number $n_0$, made of a sine function convoluted with a Gaussian function (Supp. Fig. S1) [14, 15, 37]. What distinguishes the wavelet transform from other time-frequency transforms is its constant quality factor $Q$ across frequencies (Supp. Fig. S2). This makes wavelets especially suitable for analyzing oscillatory signals with varying frequency, such as the one from CR cell undergoing a photophobic response (Fig. 1(a,b)).
## D. Phase Synchronization Index: PSI

To analyze the occurrence of synchronization, we extend the already published Phase Synchronization Index (PSI) [16, 17] to its complex form. From discrete phase signals {$\varphi_1^1, \varphi_1^2, ..., \varphi_1^N$} and {$\varphi_2^1, \varphi_2^2, ..., \varphi_2^N$}, where N refers to the size of the dataset, the PSI, noted $\tilde{\Upsilon}$ reads:

$$ \tilde{\Upsilon} = \left\langle e^{i(n\varphi_1^k - m\varphi_2^k)} \right\rangle_k, \quad (4) $$

where k is a constant, and n, m are positive integers that correspond to different harmonic factors. $|\tilde{\Upsilon}|$ ranges from 0 (no synchronization) and 1 (full synchronization). In the following, we focus on the synchronization between fundamental modes: n = m = 1.

While this standard Phase Synchronization Index (PSI) analysis described above is powerful to study synchronization over time, it assumes that oscillations have stationary frequencies and is therefore limited when applied to transient or nonstationary signals. In contrast, our refined approach described below extends the PSI to the time-frequency domain, providing precise insight into how synchronization evolves during rapid transitions such as photoshock. We therefore introduce $\tilde{\Upsilon}_{\Psi}(t, f)$, the PSI defined in the time-frequency domain by including the CWT of each signal, $W_1(t, f)$ and $W_2(t, f)$, in Eq. (4), without the need of computing explicitly their phases, which yields:

$$ \tilde{\Upsilon}_{\Psi}(t, f) = \left\langle \frac{W_1(t, f)\overline{W_2}(t, f)}{|W_1(t, f)||W_2(t, f)|} \right\rangle_{\Delta T_{\Psi}} = \Upsilon_{\Psi} e^{i\Phi}. \quad (5) $$

where the temporal averaging $\langle \cdot \rangle_{\Delta T_{\Psi}}$ is performed over the wavelet window $[t - \frac{\Delta T_{\Psi}}{2}, t + \frac{\Delta T_{\Psi}}{2}]$, with $\Delta T_{\Psi} \propto n_0/f^*$ ensuring that the PSI is computed over a timescale appropriate to the analyzed frequency. In Supplementary Fig. S4, we compare the PSI computation on a fixed size window (1 s) to this adaptive size $\Delta T_{\Psi}$. Following this framework, the phase and the modulus of $\tilde{\Upsilon}_{\Psi}$ are variables of t and f: $\Phi(t, f)$ and $\Upsilon_{\Psi}(t, f)$, respectively.

The time-frequency complex PSI can be compared to the time-frequency complex coherence that quantifies the linear coupling between two signals [17, 38, 39]:

$$ \kappa_{1,2}(t, f) = \frac{\langle W_1(t, f) \overline{W_2}(t, f) \rangle_{\Delta T_{\Psi}}}{\sqrt{\langle |W_1(t, f)|^2 \rangle_{\Delta T_{\Psi}} \langle |W_2(t, f)|^2 \rangle_{\Delta T_{\Psi}}}} \quad (6) $$

Perfect coherence between two signals at a given frequency occurs when they maintain both a constant phase difference and a constant amplitude ratio over the time interval considered. Eq. (6) differs markedly from Eq. (5). Specifically, in coherence, the averaging is performed separately on the numerator and denominator before computing the ratio. As a result, the coherence computation includes both amplitude and phase information from the full signal. In contrast, $\tilde{\Upsilon}_{\Psi}$ is computed by normalizing

out amplitude at each time-frequency point prior to averaging. This makes it particularly suited for detecting phase relationships independently of amplitude fluctuations. Importantly, if the signal amplitudes remain constant over time, as implicitly assumed in the derivation of PSI, then $\tilde{\Upsilon}_{\Psi}$ becomes equivalent to the modulus of the complex coherence $\kappa_{1,2}$.

## III. RESULTS

### A. Coupled oscillators simulations

Fig. 2 illustrates the dynamics of a system composed of two coupled oscillators governed by Eq. (2), with fixed natural frequencies $\omega_1/2\pi = 30$ Hz, $\omega_2/2\pi = 20$ Hz. The coupling strength B and initial phases $\varphi_1(0), \varphi_2(0)$ vary between cases. Steady-state synchronization is reached when the phase difference becomes constant over time, i.e. $\dot{\delta} = \dot{\varphi}_1 - \dot{\varphi}_2 = 0$. Analytically, based on Adler's equation given above, this condition is met in the limit $t \to \infty$, which means that a sufficiently strong coupling exists $|2B| > |\Delta\omega|$. In such a strong coupling regime (Fig. 2(a)), the two oscillators quickly synchronize (left panel) with a constant phase difference (right panel) indicating in-phase synchronization ($\delta = 0$). In contrast, in a weak coupling condition the system does not maintain synchronization ($|2B| < |\Delta\omega|$, Fig. 2(b)). Although transient locking occurs, the oscillators remain predominantly desynchronized, with a quasi periodic behavior. Note that the sign of B determines the nature of the synchronization: the two oscillators can synchronize in phase if B is negative (Fig. 2(a)) or in anti-phase if B is positive ($\delta = \pi$, Fig. 2(c)). Furthermore, even in the synchronized time-region, the instantaneous frequencies of the oscillators vary. This occurs because the synchronization frequency $\omega_S = \dot{\varphi}_1 = \dot{\varphi}_2 = (\omega_1 + \omega_2)/2$ generally differs from the natural frequency of each uncoupled oscillator. These transient dynamics highlight the importance of time-frequency methods to detect synchronization in nonstationary systems.

To better illustrate nonstationary dynamics, in addition of a positive coupling, we now define a system of two chirp oscillators whose natural frequencies increase linearly over time (Fig. 3(a)). This supposedly mimics realistic biological or physical systems where oscillator properties evolve gradually due to internal modulation or external perturbations.

We pick an intermediate value for coupling ($B = 2\Delta\omega$) in Eq. (2) to favor longer synchronization states. The spectrograms of the two oscillators are shown in Fig. 3(b-c), where black and gray lines display the instantaneous frequency $f^*(t)$, obtained as the local maxima of the wavelet modulus, also called “ridges”, which verify $\partial_f|W(t, f^*)| = 0$ [40].

Obviously, this computation enables us to follow the linear increase of each oscillator frequency. While these two signals are mostly synchronized, few phase slips are
<figure><img src="image_4.png" /><figcaption>FIG. 3. System of two coupled oscillators with linearly increasing frequency. (a) Simulated signals generated from Eq. (2). Gray: cos($\varphi_{1}$); black: cos($\varphi_{2}$). Frequencies: $\omega_{1}/2\pi = 5$Hz, $\omega_{2}/2\pi = 4$Hz. Coupling strength: $B_{1}/2\pi = B_{2}/2\pi = 0.5$Hz (b) $|W_{1}(t, f)|$ for the first oscillator (gray line in (a)) using Morlet wavelet with $n_{0} = 6$. The gray line indicates the local frequency maxima $f_{1}^{*}(t)$. (c) $|W_{2}(t, f)|$ for the second oscillator (black line in (a)) using Morlet wavelet with $n_{0} = 6$. The black line indicates the local frequency maxima $f_{2}^{*}(t)$. (d) $\Upsilon_{\Psi}(t, f)$ computed with an adaptive time window $\Delta T_{\Psi} \propto n_{0}/f^{*}(n_{0} = 6)$. Frequency ridges $f_{1}^{*}(t)$ (gray) and $f_{2}^{*}(t)$ (black) are overlaid. (e) $\Upsilon_{\Psi}(t, f_{12}^{*}(t))$ (black), computed along the frequency trajectory $f_{12}^{*}(t) = (f_{1}^{*}(t) + f_{2}^{*}(t))/2$. The derivative of the phase difference between $S_{1}$ and $S_{2}$, $\dot{\delta}_{W}(t)$ (gray) was computed using a Gaussian derivative wavelet to filter out the signal noise (Supplementary Fig. S5) [41].</figcaption></figure>

visible on the time-domain signals (Fig. 3(a)) (from 3 s to 3.2 s for instance). The time-frequency maps confirm the presence of transient phase slips, displaying frequency mismatch (lines in Fig. 3(d)).

Fig. 3(d) shows the time-frequency map of the PSI $\Upsilon_{\Psi}(t, f)$ for the two coupled chirp oscillators. For reference, the instantaneous frequency ridges $f_{1}^{*}(t)$ (gray) and $f_{2}^{*}(t)$ (black), previously extracted in Fig. 3(b–c), are overlaid as guide lines. Regions of low $\Upsilon_{\Psi}$ appear as localized “islands” in this map, marking transient losses of synchronization that are clearly bounded by the main frequency ridges.

To quantify the synchronization dynamics more directly, we compute $\Upsilon_{\Psi}(t, f_{12}^{*}(t))$ along the averaged frequency trajectory $f_{12}^{*}(t) = [f_{1}^{*}(t) + f_{2}^{*}(t)]/2$ (black line, Fig. 3(e)). On the same plot, we display the temporal derivative of the phase difference between the two oscillators, $\dot{\delta}_{W}(t)$ (gray curve), which provides a complementary measure of phase locking, and therefore a validation of phase synchronization. Both quantities were computed directly in the wavelet domain: $\dot{\delta}_{W}(t)$ was obtained using a Gaussian derivative wavelet, which efficiently filters

out high-frequency noise while preserving the temporal localization of phase jumps (Supplementary Fig. S5).

The comparison between $\Upsilon_{\Psi}$ and $\dot{\delta}_{W}$ in Fig. 3(e) confirms that regions of strongest synchronization ($\Upsilon_{\Psi} \to 1$) correspond to the lowest phase-derivative values ($\dot{\delta}_{W} \to 0$). For comparison, Supplementary Fig. S4 shows $\Upsilon_{\Psi}$ computed with a constant time window; in this case, the “islands” of desynchronization vanish, and synchronization states are poorly resolved, demonstrating the advantage of the adaptive formulation. Furthermore, applying this framework to the coupled oscillator model of Fig. 2 reproduces the same qualitative behavior (Supplementary Fig. S3), confirming the robustness of the method. Together, these results show that the proposed time-frequency PSI method provides an accurate and temporally resolved characterization of synchronization dynamics, thereby motivating its application to biological oscillators.
## B. Flagellar resynchronization in response to a photophobic stimulus

The time response of CR to a photoshock unfolds in three distinct phases (Fig. 4(a)). First, before stimulation ($t < 0$ ms), the cell swims in a breaststroke motion in which the two flagella are synchronized, resembling two noisy oscillators with strong coupling, as described in Fig. 2(a,c). The velocities recorded along the x-axis are in phase, while the ones along the y-axis are in anti-phase. Then, immediately after the photoshock ($0.05 s < t < 0.6$ s), the beating switches to an undulatory mode, which drives backward swimming with smaller amplitude. In this regime, the two flagella oscillate in phase along the y-axis. It should be noted that during this phase, the signals detected at $Y = \pm 1$ reflect the combined flow generated by both flagella, which can no longer be distinguished individually. However, once a flagellum resumes its breaststroke-like motion, the signal from the nearest window recaptures its specific dynamics, while the opposite window no longer detects it. Lastly, for $t \gtrsim 0.6$ s, the breaststroke pattern is progressively restored, until full resynchronization.

To quantitatively characterize the dynamics of the transitions between these different stages, we apply our wavelet-based time-frequency analysis to the velocity signals (Fig. 4(b-c)). This approach allows us to clearly identify the three behavioral phases described above, while providing detailed insight into the instantaneous frequency evolution of the flagellar beating. Both x- and y-axis velocity components (left and right panels, respectively) show qualitative similar trends, while their relative amplitudes reflect the dominance of different swimming modes at each stage.

The wavelet maxima frequency (marked with lines on the spectrograms of Fig. 4(b-c)) dynamics follows a reproducible sequence. In the pre-stimulus phase, both flagella exhibit steady beating frequencies centered around 40 Hz (Fig. 4(b-c)). Following the photoshock, the beating frequency rises sharply to 70 Hz, consistent with the undulatory backward-swimming mode. Then, it gradually decreases to around 50 Hz over several hundred milliseconds. This transition is followed by a sharp drop into a low-frequency regime, approximately 30 Hz, marking the start of resynchronization. The beating frequency then progressively increases and stabilizes near the pre-stimulus breaststroke value, indicating the recovery of the original synchronized state.

To quantify phase synchronization more precisely, we apply our PSI framework to these flagellar signals, using Eq. (5). Instantaneous frequency maxima $f_{-1}^*$ and $f_{+1}^*$ were extracted from $|W_{-1}(t, f)|$ and $|W_{+1}(t, f)|$, respectively (lines in Fig. 4(b-c)). The subscripts $\pm 1$ refer to the averaging windows located on the two opposite sides of the cell body ($Y = \pm 1$; see Fig. 1(a)). $\Upsilon_\Psi(t, f_{\pm 1}^*)(t)$ is evaluated using an adaptive time window such as $f_{\pm 1}^*(t) = (f_{-1}^* + f_{+1}^*)/2$ (Fig. 4(d)). This analysis displays two desynchronization events: imme-

diately following the photoshock (from 0 to 0.2 s), and after the stimulus, roughly from 0.5 to 1 s. While the first is linked to an abrupt change in flagellar beating modes, the second reflects a subtler change in phase synchronization. This is further confirmed by the PSI phase $\Phi(t, f_{\pm 1}^*)$, computed as the phase of the complex $\tilde{\Upsilon}_\Psi$ (Equation (5)), which transitions along $y$ from in-phase ($\Phi \approx 0$) to anti-phase ($\Phi \approx \pi$) during the recovery (Fig. 4(e)). This gradual shift highlights reorganization of the underlying coupling dynamics between the two flagella. Finally, while Fig. 4 shows a representative case, we confirmed robustness across the dataset by averaging $\Upsilon_\Psi(t, f_{\pm 1}^*)$ over all recorded photoshock events (Fig. 5(a,b)), retaining only experiments with strong pre- and post-shock synchronization. More precisely, the PSI verifies $\langle\Upsilon_\Psi\rangle_t > 0.85$, where the temporal averaging $\langle\cdot\rangle_t$ excludes the photoshock dynamics: $t \notin [-1, +2]$ s.

Beyond tracking fundamental frequencies, the time-frequency maps in Fig. 4(b-c) also reveal higher-harmonic structure in the flagellar beating. Before the shock, two main frequency bands are present: a fundamental at $\sim 40$ Hz and its harmonic at $\sim 80$ Hz (Fig. 4(b,c) and 5(c,d)). The lower mode carries most of the energy, especially along the x-axis (see color scale in Fig. 5(c,d)). Immediately following the photoshock, energy redistribution occurs: along the x-axis both modes remain but with transient amplitude dynamics (Fig. 5(c)), while the low-frequency mode disappears from the y-axis, leaving the harmonic dominant (Fig. 5(d)). These modes evolve dynamically over the course of few hundreds of milliseconds until steady state is reached. These observations suggest a subtle mechanism in which the persistent presence of higher harmonic modes during breaststroke beating may provide a way for smooth switching into backward swimming under stress. This pre-existing spectral component may thus act as a “reserve” mode that facilitates fast and efficient photophobic responses.

## IV. DISCUSSION

Using a time-frequency framework combined with a wavelet-based Phase Synchronization Index (PSI), we investigated the photophobic response of *Chlamydomonas reinhardtii* cells held in place by a micropipette (Fig. 1). This approach enabled us to follow the evolution of flagellar beating across different swimming stages, with high temporal and spectral resolution, and to capture how synchronization breaks and re-emerges. Our analysis shows that flagellar coordination does not switch abruptly between discrete modes, but instead involves a gradual redistribution of energy across coexisting oscillatory components. In particular, we find that the fundamental breaststroke frequency temporarily vanishes during photoshock, while a higher-frequency mode becomes dominant, before both modes re-establish during resynchronization. The persistence of these harmonic modes, even during breaststroke beating, may provide a “spec-
<figure><img src="image_5.png" /><figcaption>FIG. 4. Flagellar resynchronization of CR following a photoshock. (a) Averaged fluid velocity signals along the x-axis (left) and y-axis (right), extracted from regions on each side of the alga: $Y = -1$ (gray) and $Y = +1$ (black), as shown in Fig. 1. (b-c) Time-frequency representation of the CWT signals, $log|W_Y(t, f)|$ (Morlet wavelet with $n_0 = 12$), for $Y = -1$ (b) and $Y = +1$ (c). Frequency ridges $f_{-1}^*(t)$ (gray) and $f_{+1}^*(t)$ (black) are also plotted. (d) Phase Synchronization Index, $|\tilde{\Upsilon}_{\Psi}|(t, f_{\pm1}^*(t))$, computed at the mean instantaneous frequency $f_{\pm1}^*(t) = [f_{-1}^*(t) + f_{+1}^*(t)]/2$. (e) Phase $\Phi(t)$ obtained from the complex PSI $\tilde{\Upsilon}_{\Psi}$. In all panels, the shaded gray region marks the duration of the photoshock.</figcaption></figure>

that facilitates rapid switching into backward swimming under stress. So far, our PSI analysis has focused on the fundamental modes ($n = m = 1$; Eq. 4), but the method can be readily extended to probe synchronization between higher-order components ($n, m > 1$). This opens the possibility of testing whether harmonics themselves are phase-locked, and how their relative weights evolve during photoshock. This richer spectral organization resonates with studies suggesting that multiple internal oscillatory states coexist in CR flagella beating [23, 25].

At a broader level, our measurements are based on flow fields generated by the two flagella rather than their direct motion. As a result, the recorded signals represent an effective measure of flagellar coordination but do not allow us to isolate the detailed contributions of each flagellum or identify the specific coupling pathways responsible for resynchronization. Direct imaging of flagellar motion, ideally resolving *cis* and *trans* beating patterns, would help bridge this gap and allow a more detailed comparison with models that explicitly include amplitude dynamics and hydrodynamic interactions [10, 22, 24, 25, 42]. In this context, applying our

time-frequency synchronization framework to mutants such as *ptx1*, which exhibit both in-phase and anti-phase beating modes, with the latter occurring at higher frequencies, could provide a valuable test case. Interestingly, this would perhaps lead us to quantify the energy loss or recovery of algae in the different swimming stages: does the algae select a swimming mode on the basis of a criterion of economy or a criterion of adaptation? Such analyses would enable direct investigation of how coupling strength and frequency asymmetry influence synchronization and gait selection.

Together, these results demonstrate the power of adaptive time-frequency methods to resolve transient synchronization dynamics. They also highlight new questions on the interplay between intrinsic oscillatory modes, amplitude regulation, and coupling mechanisms in eukaryotic flagella. While demonstrated here on flagellar dynamics, this approach may broadly be applied to other biological or physical systems where oscillatory coordination evolves in time.
V. ACKNOWLEDGMENTS

We thank Ahmad Badr for insightful discussions. The authors acknowledge support of the Leverhulme Trust (grant RPG-2018-345; A.A. and M.P.) ; the SMR Department and CNRS Tremplin@INP (A.A.) ; MICINN (M.P.: grant PID2019-104232BI00/AEI; J.A.: grant PID2022-143018NB-I00/AEI). M.P. acknowledges the fact that IMEDEA is an accredited ‘María de Maeztu Excellence Unit’ (grant CEX2021-001198, funded by MCIN/AEI/10.13039/501100011033).

A.A., J.A. and M.P. conceived the project. A.A. and J.A. designed and performed the experiments. L.F., A.A. and J.A. developed the analysis code. L.F. and F.A. designed and performed the signal analysis work. L.F., A.A. and F.A. wrote the manuscript with feedback from all authors.

The data and code supporting this work are avail-
able in Zenodo (https://doi.org/10.5281/zenodo.
17406887).

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

FIG. 5. Flagellar synchronization and frequency dy-
namics around photoshock. (a-b) PSI averaged over
44 photoshocks (13 algae) computed using the two signals
Y = -1 and Y = 1, along the x-axis (a) and the y-axis
(b). Black line: Median PSI. Gray shade: First quartile (Q1)
and Third quartile (Q3). (c-d) The red lines correspond to
the averaged wavelet frequency maxima of each band accross
all the algae, computed using signal in Y = 0. The gray
shades are the standard deviations. The color map is linked
to the amplitude and represents the average wavelet transform
squared modulus of the corresponding frequency maxima, in
log scale. Frequencies between t = -50 ms and t = 100 ms
are not shown because they are not robustly detected using
the wavelet transform.

[1] M. Riedl, I. Mayer, J. Merrin, M. Sixt, and B. Hof, "Synchronization in collectively moving inanimate and living active matter," *Nature Communications*, vol. 14, p. 5633, Sept. 2023. Publisher: Nature Publishing Group.

[2] E. M. Izhikevich, *Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting*. The MIT Press, July 2006.

[3] C. Schäfer, M. G. Rosenblum, H.-H. Abel, and J. Kurths, "Synchronization in the human cardiorespiratory system," *Physical Review E*, vol. 60, pp. 857-870, July 1999. Publisher: American Physical Society.

[4] U. Schibler, “Steven A. Brown and the synchronization of circadian rhythms by body temperature cycles,” *European Journal of Neuroscience*, vol. 60, no. 2, pp. 3891–3900, 2024. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/ejn.16431.

[5] S. H. Strogatz and I. Stewart, "Coupled Oscillators and Biological Synchronization," *Scientific American*, vol. 269, no. 6, pp. 102–109, 1993. Publisher: Scientific American, a division of Nature America, Inc.

[6] C. M. Denis and P. François, “Unclocklike Oscillators with Frequency Memory for the Entrainment of Biological Clocks,” *PRX Life*, vol. 2, p. 043008, Oct. 2024. Publisher: American Physical Society.

[7] C. Huygens, *Horologium oscillatorium sive De motu pendulorum ad horologia aptato demonstrationes geometricae*. 1673.

dulorum ad horologia aptato demonstrationes geometri-
cae. 1673.

[8] A. Pikovsky, Author, M. Rosenblum, J. Kurths, and R. C. Hilborn, “Synchronization: A Universal Concept in Nonlinear Science,” *American Journal of Physics*, vol. 70, p. 655, June 2002.

[9] B. Button, L.-H. Cai, C. Ehre, M. Kesimer, D. B. Hill, J. K. Sheehan, R. C. Boucher, and M. Rubinstein, “A Periciliary Brush Promotes the Lung Health by Separating the Mucus Layer from Airway Epithelia,” *Science*, vol. 337, pp. 937–941, Aug. 2012. Publisher: American Association for the Advancement of Science.

[10] S. Hu and F. Meng, “Multiflagellate Swimming Controlled by Hydrodynamic Interactions,” *Physical Review Letters*, vol. 132, p. 204002, May 2024. Publisher: American Physical Society.

[11] E. M. Holland, H. Harz, R. Uhl, and P. Hegemann, “Control of phobic behavioral responses by rhodopsin-induced photocurrents in Chlamydomonas,” *Biophysical Journal*, vol. 73, pp. 1395–1401, Sept. 1997.

[12] A. Guillet and F. Argoul, “Uncertainty and information in physiological signals: Explicit physical trade-off with log-normal wavelets,” *Journal of the Franklin Institute*, vol. 361, p. 107201, Dec. 2024.

[13] B. Torrésani, *Analyse continue par ondelettes*. Savoirs
actuels / Série Physique, CNRS Editions / EDP Sciences,
1995.

[14] C. Torrence and G. P. Compo, “A Practical Guide to Wavelet Analysis,” *Bulletin of the American Meteorological Society*, vol. 79, pp. 61–78, Jan. 1998.

[15] R. Kronland-Martinet, J. Morlet, and A. Grossmann, “Analysis of sound patterns through wavelet transforms,” *International Journal of Pattern Recognition and Artificial Intelligence*, vol. 01, pp. 273–302, Aug. 1987. Publisher: World Scientific Publishing Co.

[16] B. Schack and S. Weiss, “Quantification of phase synchronization phenomena and their importance for verbal memory processes,” *Biological Cybernetics*, vol. 92, pp. 275–287, Apr. 2005.

[17] Z. Liang, Y. Bai, Y. Ren, and X. Li, “Synchronization Measures in EEG Signals,” in *Signal Processing in Neuroscience* (X. Li, ed.), pp. 167–202, Singapore: Springer Singapore, 2016.

[18] V. V. Grubov, A. K. Kuc, S. A. Kurkin, D. A. Andrikov, N. Utyashev, V. A. Maksimenko, O. E. Karpov, and A. E. Hramov, “Harnessing Long-Range Temporal Correlations for Advanced Epilepsy Classification,” *PRX Life*, vol. 3, p. 013005, Jan. 2025.

[19] R. Adler, “A study of locking phenomena in oscillators,” *Proceedings of the IRE and Waves and Electrons*, vol. 34, pp. 351–357, 1946.

[20] Y. Kuramoto, *Chemical Oscillations, Waves, and Turbulence*. Springer Series in Synergetics, Berlin Heidelberg: Springer-Verlag, 1984. ZSCC: NoCitationData[s0].

[21] E. H. Harris, *The Chlamydomonas Sourcebook: Introduction to Chlamydomonas and Its Laboratory Use: Volume 1*. Academic Press, 1989. Google-Books-ID: xTjJGV5GWY0C.

[22] Y. Liu, R. Claydon, M. Polin, and D. R. Brumley, “Transitions in synchronization states of model cilia through basal-connection coupling,” *Journal of The Royal Society Interface*, vol. 15, p. 20180450, Oct. 2018. Publisher: Royal Society.

[23] K. Y. Wan and R. E. Goldstein, “Coordinated beating of algal flagella is mediated by basal coupling,” *Proceedings of the National Academy of Sciences*, vol. 113, May 2016.

[24] L. Zorrilla, A. Allard, K. Desai, and M. Polin, “The role of hydrodynamics in the synchronisation of {\it Chlamydomonas} flagella,” *Apr. 2025*. arXiv:2504.02680 [physics].

[25] B. Friedrich, “Hydrodynamic synchronization of flagellar oscillators,” *The European Physical Journal Special Topics*, vol. 225, pp. 2353–2368, Nov. 2016.

[26] R. E. Goldstein, M. Polin, and I. Tuval, “Noise and Synchronization in Pairs of Beating Eukaryotic Flagella,” *Physical Review Letters*, vol. 103, p. 168103, Oct. 2009. Publisher: American Physical Society.

[27] P. Hegemann and B. Bruck, “Light-induced stop response in Chlamydomonas reinhardtii: Occurrence and adaptation phenomena,” *Cell Motility*, vol. 14, no. 4, pp. 501–515, 1989.

[28] U. Rüffer and W. Nultsch, “Flagellar photoresponses of Chlamydomonas cells held on micropipettes: II. Change in flagellar beat pattern,” *Cell Motility*, vol. 18, no. 4, pp. 269–278, 1991.

[29] J. Arrieta, A. Barreira, M. Chioccioli, M. Polin, and I. Tuval, “Phototaxis beyond turning: persistent accumulation and response acclimation of the microalga Chlamydomonas reinhardtii,” *Scientific Reports*, vol. 7,

p. 3447, June 2017. Number: 1 Publisher: Nature Publishing Group.

[30] A. Ramamonjy, J. Dervaux, and P. Brunet, “Non-linear Phototaxis and Instabilities in Suspensions of Light-Seeking Algae,” *Physical Review Letters*, vol. 128, p. 258101, June 2022.

[31] N. Ueki and K.-i. Wakabayashi, “Dynein-mediated photobehavioral responses in *Chlamydomonas*,” in *Dyneins: Structure, Biology and Disease (Second Edition)* (S. M. King, ed.), pp. 368–385, Academic Press, Jan. 2018.

[32] C. Brennen and H. Winet, “Fluid Mechanics of Propulsion by Cilia and Flagella,” *Annual Review of Fluid Mechanics*, vol. 9, pp. 339–398, Jan. 1977. Publisher: Annual Reviews.

[33] A. K. Boggon, A. D. Hastewell, J. Dunkel, and K. Y. Wan, “Embodied behavioural complexity in a ciliated microorganism,” Aug. 2025. ISSN: 2692-8205 Pages: 2025.08.20.671109 Section: New Results.

[34] M. Riccomi, F. Alberini, E. Brunazzi, and D. Vigolo, “Ghost Particle Velocimetry as an alternative to µPIV for micro/milli-fluidic devices,” *Chem. Eng. Res. Des.*, vol. 133, pp. 183–194, May 2018. Publisher: Elsevier BV.

[35] M. J. Aburn, C. A. Holmes, J. A. Roberts, T. W. Boonstra, and M. Breakspear, “Critical Fluctuations in Cortical Models Near Instability,” *Frontiers in Physiology*, vol. 3, Aug. 2012. Publisher: Frontiers.

[36] F. Roddier, *Distributions et transformation de Fourier : A l'usage des physiciens et des ingénieurs*. New York NY: McGraw-Hill, 1984.

[37] M. X. Cohen, “A better way to define and describe Morlet wavelets for time-frequency analysis,” *NeuroImage*, vol. 199, pp. 81–86, Oct. 2019.

[38] A. Bruns, “Fourier-, Hilbert- and wavelet-based signal analysis: are they really different approaches?,” *Journal of Neuroscience Methods*, vol. 137, pp. 321–332, Aug. 2004.

[39] A. Guillet, A. Arneodo, and F. Argoul, “Tracking Rhythms Coherence From Polysomnographic Records: A Time-Frequency Approach,” *Frontiers in Applied Mathematics and Statistics*, vol. 7, Apr. 2021. Publisher: Frontiers.

[40] R. Carmona, W. L. Hwang, and B. Torrésani, “Identification of Chirps with Continuous Wavelet Transform,” in *Wavelets and Statistics* (A. Antoniadis and G. Oppenheim, eds.), pp. 95–108, New York, NY: Springer, 1995.

[41] S. Mallat, *A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way*. USA: Academic Press, Inc., 3rd ed., Nov. 2008.

[42] D. R. Brumley, K. Y. Wan, M. Polin, and R. E. Goldstein, “Flagellar synchronization through direct hydrodynamic interactions,” *eLife*, vol. 3, p. e02750, July 2014. Publisher: eLife Sciences Publications, Ltd.
# VI. SUPPLEMENTARY MATERIALS

## A. Python packages used for the simulations

The Python modules used to model the are `numpy`, `scipy`, `scipy.integrate-odeint`, `sdeint`, `pycwt`, `PyWavelets`, `mat-plotlib`, `pandas`.

## B. Morlet Wavelet transform

<figure><img src="image_7.png" /><figcaption>FIG. S1. Wavelet transform of a signal for different wave numbers: $n_0 = 6$ (left) and $n_0 = 12$ (right). (a) Re{$\Psi_{n_0}$}. (b) Signal: 1Hz sine function (0s-4s) combined to a 3Hz sine function (4s-8s). (c) $|W(t, f)|$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. S2. Real part of the Morlet wavelet in time (left) and spectral (right) domain. From top to bottom, the scale parameter varies from $a = f^*/f = 1$ to 1/4. The Morlet mother wavelet for CWT is a sine function convoluted with a Gaussian function, given by $\Psi_{n_0}(t) = \pi^{-1/4}e^{in_0t}e^{-t^2/2}$, where $n_0 > 5$ represents the wave number of the wavelet, proportional to the number of oscillations of the mother wavelet in time [14]. What distinguishes the wavelet transform from other time-frequency transforms is its constant quality factor $Q$ across frequencies (right panels). The $Q$-factor is defined by $Q = f^*/\Delta f$ where $f^*$ is the frequency of the maximum of the mother wavelet in Fourier space and $\Delta f$ its Full-Width-Half-Maximum. Increasing $n_0$ results in a better frequency accuracy (higher $Q$-factor) at the expense of time accuracy.</figcaption></figure>
C. Time frequency PSI in dynamical systems

Fig. S4(d) displays a time frequency PSI: $Υ$, where temporal averaging is performed over a constant time window $ΔT = 1$s. This basic PSI computation does not take into account the dynamics of the system and loses its precision over time when the overall frequencies of the two oscillators increase. This figure illustrates the relevance of computing the PSI over a real-time adaptive wavelet window $ΔT_Ψ$ (Fig. S4(e)) using Eq. 5, where $ΔT_Ψ = [ΔT_Ψ(f_1^*)+ΔT_Ψ(f_2^*)]/2$.

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

FIG. S3. System of two coupled oscillators with Gaussian noise, modeled by Eq. (2). Left column corresponds to strong negative coupling ($B/2\pi = -20$ Hz $< 0$, $|2B| > \Delta\omega$), with $\varphi_1(0) = 0$, $\varphi_2(0) = \pi$: leads to stable in-phase synchronization (same data as Fig. 2(a)). Central column corresponds to weak negative coupling ($B/2\pi = -3$ Hz $< 0$, $|2B| < \Delta\omega$), with same initial phases: leads to desynchronization with transient coordination (same data as Fig. 2(b)). Right column corresponds to strong positive coupling ($B/2\pi = 20$ Hz $> 0$, $|2B| > \Delta\omega$) with $\varphi_1(0) = \varphi_2(0) = 0$: leads to stable anti-phase synchronization. (same data as Fig. 2(c)). (a) Simulated signal following Eq. (2). Gray: $\cos(\varphi_1)$. Black: $\cos(\varphi_2)$. (b) $|W_1(t, f)|$, $n_0 = 6$. Gray line: local frequency maxima $f_1^*(t)$. (c) $|W_2(t, f)|$, $n_0 = 6$. Black line: local frequency maxima $f_2^*(t)$. (d) $Υ_Ψ(t, f)$ computed with a variable window $ΔT_Ψ = n_0/f^*$ (Eq. 5), $n_0 = 6$. Gray: $f_1^*(t)$. Black: $f_2^*(t)$. (e) Plain line: $Υ_Ψ(t, f_{12}^*(t))$ where $f_{12}^*(t) = (f_1^*(t) + f_2^*(t))/2$. Note that the influence cone is visible in panels (d-e), more details in [13].
<figure><img src="image_10.png" /><figcaption>FIG. S4. System of two coupled oscillators with a linear frequency evolution (chirp). (a) Simulated signal following Eq. (2). Gray: $cos(\varphi_1)$. Black: $cos(\varphi_2)$. Frequencies: $\omega_1/2\pi = 5$ Hz, $\omega_2/2\pi = 4$ Hz. Coupling strength: $B_1/2\pi = B_2/2\pi = 0.5$ Hz (b) $|W_1(t, f)|$, $n_0 = 6$. Gray line: local frequency maxima $f_1^*(t)$. (c) $|W_2(t, f)|$, $n_0 = 6$. Black line: local frequency maxima $f_2^*(t)$. (d) $\Upsilon(t, f)$ computed with a constant window $\Delta T = 1$ s. Gray: $f_1^*(t)$. Black: $f_2^*(t)$. (e) $\Upsilon_{\Psi}(t, f)$ computed with a variable window $\Delta T_{\Psi} = n_0/f^*$ (Eq. 5), $n_0 = 6$. Gray: $f_1^*(t)$. Black: $f_2^*(t)$. (f) Plain line: $\Upsilon_{\Psi}(t, f_{12}^*(t))$ where $f_{12}^*(t) = (f_1^*(t) + f_2^*(t))/2$. Dashed line: $\Upsilon(t, f_{12}^*(t))$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. S5. Temporal derivative of the δ computed using a Gaussian derivative wavelet [41]. (a) δ(t) using Fig. 3 signal. (b) Wavelet transform modulus of δ(t) shown in (a) with a Gaussian derivative ψ̂$_{G}$ mother wavelet. We have the relation:</figcaption></figure>

$$W_{\psi_{\hat{G}}}[\delta](t, f) \propto W_{\psi_{G}}[\delta](t, f).$$

The horizontal lines correspond to different frequencies (20, 60, and 100 Hz). (c) Profiles $\dot{\delta}_W(t) = W_{\psi_{\hat{G}}}[\delta](t, f)$ for the frequencies displayed in (b). (d) Classical derivative of the signal $\dot{\delta}(t) \sim (\delta(t + dt) - \delta(t))/dt$ with $dt = 10$ ms. We note that the scale of the wavelet changes the amplitude of the phase jumps (singularities of $\delta$), the smaller the scale ($af^*/f$), the sharper the local extrema of $\dot{\delta}$. However the temporal localization of these jumps can still be achieved for larger wavelet scales.
D. Breaststroke frequency of selected algae

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

FIG. S6. Distribution of the mean beating frequency observed in breaststroke swimming (1 Hz bins). The mean of the Wavelet transform Modulus of the $u_x$ signal (Y=0), is computed on breaststroke swimming regions and the frequency associated to the maxima is taken. Total of 44 frequencies, associated to 44 photoshock runs on 13 different cells. (a) Average frequency distribution before the photoshock (between $t = -1000$ ms and $t = -50$ ms [Fig. 5]). (b) Average frequency distribution after the photoshock (between $t = 2000$ ms and $t = 3000$ ms [Fig. 5]).