Nonreciprocal synchronization in embryonic oscillator ensembles

Significance Synchronization within networks of coupled oscillators is a universal phenomenon found across widely different contexts, such as the unison applause we witness in concert halls and fireflies that flash in synchrony. Despite apparent differences in scale and context, synchronization depends on a set of general coupling rules that are universal. Different rules for oscillating couplings lead to different emerging collective behaviors. Here, we combine experiments and theory to identify the coupling rules for an embryonic oscillator cell ensemble. Our results reveal “winner-takes-it-all” synchronization dynamics and more generally, the importance of nonreciprocal cell interactions as a general theme in biology.

From the generic Winfree model of phase oscillator coupling to Kuramoto's phase-difference dependent coupling function.We start with the following differential equations to describe the phase dynamics of coupled oscillators A and B (1): φA = ωA + c Z(ϕA) Γ(ϕA, ϕB) φB = ωB + c Z(ϕB) Γ(ϕB, ϕA) [1] where c represents the coupling strength, and Z describes the response of each oscillator to the external perturbation Γ imposed by the other oscillator.We assume that the coupled oscillators reach a common frequency Ω, and perform the following change of variables: Injecting this Ansatz into Eq. 1 gives the differential equations for ΨA and ΨB: ΨA = ωA − Ω + c Z(ΨA + Ω t) Γ(ΨA + Ω t, ΨB + Ω t) ΨB = ωB − Ω + c Z(ΨB + Ω t) Γ(ΨB + Ω t, ΨA + Ω t) [3] This change of variable allows us to simplify the coupling terms.Indeed, if ΨA and ΨB vary slowly, we can safely average Eq.
3 over one period T = 2π/Ω to get the following equations: where Z(u) Γ(u, ΨB − ΨA + u) du [5] where we performed a change of variable u = Ωt + ΨA to eliminate the Ω dependency, using the fact that Z and Γ are 2π-periodic in their phase arguments to perform the period averaging.This procedure renders explicit the dependency on the phase difference ΨB − ΨA.Since ΨB − ΨA = ϕB − ϕA, we can also go back to the initial variables: φA = ωA + c H(ϕB − ϕA) φB = ωB + c H(ϕA − ϕB) [6] which are Eq. 3 & 4 given in the main text.It is of course always possible to study directly such system of equations, but notice it is only equivalent to the full version when both oscillators are assumed to lock to a common frequency.In his 1984 paper, Kuramoto studied the synchronization dynamics of Eq. 6 with H(∆ϕ) = sin(∆ϕ) (2).A follow-up paper published 2 years later with Sakaguchi extended the analysis to sine functions with any offset α: H(∆ϕ) = sin(∆ϕ + α) (3).

Construction of the rectified Kuramoto (ReKu) coupling function.
To determine which properties the H function must have to produce winner-takes-it-all synchronization, we first note that the H function directly governs the temporal evolution of a given coupled oscillator relative to its reference oscillator.Indeed, since the temporal evolution of the reference oscillators is respectively given by φr A = ωA and φr B = ωB, we have: The winner-takes-it-all dynamics means that at least one of those equations is always 0, irrespective of the (changing) values of ϕA and ϕB.Formally, setting ∆ϕ = ϕA − ϕB, this means that for all values of ∆ϕ taken experimentally.Eq. 9 implies that H(∆ϕ = 0) = 0, which is consistent with the fact that oscillators eventually synchronize with no change of intrinsic frequencies (Fig. S4A).Furthermore, Eq. 9 suggests that H has a singular, asymmetric behaviour close to 0. To see this, let us assume that H is differentiable in 0, meaning that sufficiently close to |∆ϕ| = 0, H(∆ϕ) ∼ A∆ϕ n where n is the order of the first non-zero derivative close to 0. Then, Eq. 9 for non-zero ∆ϕ imposes that A is 0, which implies that H = 0 everywhere.To solve this conundrum, the simplest solution is to assume that H is not differentiable in 0, and with Eq. 9 the simplest hypothesis is to assume that H is linear on one side and identically 0 on the other side.There are thus four possible qualitative behaviours depending on the sign of the derivative in 0 and the choice of the non-zero side of H.Only two of these four choices give stable synchronized coupling: where a > 0 ensures that ∆ϕ = 0 is stable, which means that synchronized oscillators stay synchronized.Periodic generalization of this rectified linear behaviour close to 0 is achieved by considering only the rectified part of the Kuramoto coupling function i.e., taking: H + (∆ϕ) = a max(0, sin(∆ϕ)) or H − (∆ϕ) = a min(0, sin(∆ϕ)) [11] as indicated in Eq. 6 and 7 of the main text.
Delayed coupling and the Kuramoto-Sakaguchi model.The goal of this section is to show that, in our theoretical framework, models based on delayed Kuramoto coupling are equivalent to the Kuramoto-Sakaguchi model.Following (4), we consider phase oscillators coupled via the Kuramoto coupling function with a delay: where τ accounts for various biochemical delays.Considering two coupled oscillators with a common frequency ω, we get: To study the coupling dynamics, we assume that both oscillators end up oscillating with a common frequency Ω, and we define: where ΨA and ΨB respectively define the relative phases of oscillators A and B close to stationarity (i.e. they vary very slowly compared to Ω, similar to what is done in the first section of this Supporting Information Text).The differential equations for ΨA and ΨB are: Taking the difference, at stationarity we get: In general, Eq. 16 implies that for positive c, the only stable phase relation is ΨB − ΨA = 0 mod 2π.Injecting this phase relation into equation 15, we get (4): This is a transcendental equation and there is no analytical solution, but, assuming that τ is small, we get the following approximate expression: This shows that if the delay is high enough, the coupled oscillators go slower than their "intrinsic" frequency ω.Now, we inject Eq. 17 into Eq.13: Finally, we use Eq. 14 close to stationarity, i.e. assuming that ΨA and ΨB are almost constant (ΨA(t) ∼ ΨA(t − τ ) ∼ ΨA), to write : which allows us to get our final result: Therefore, close to stationarity, introducing a delay τ in the Kuramoto coupling function corresponds to using the Kuramoto-Sakaguchi coupling function, H(∆ϕ) = sin(∆ϕ + α) − sin(α), with a value of α = −Ωτ (see Eq. 8 of the main text).Using the parameters from (4), i.e a delay of τ = 7 min, and a collective period of T = 39.1 min, corresponding to a collective frequency of Ω = 2π/T = 0.161 min −1 , we get a value of α = −1.13.As explained in the main text, the Kuramoto-Sakaguchi model with values of α close to −1 display winner-takes-it-all dynamics for a large interval of initial conditions for the two coupled oscillators' phases (Fig. 5B), but fails to generate winner-takes-it-all dynamics for initial phase differences close to π (Fig. S8A), in contrast with the RAFL data (Fig. 3B, second row, and Fig. S5B).Pulsed-coupling models and piece-wise linear versions.To build the pulsed-coupling model, we use Winfree's theoretical framework (see (1), and the first section of this Supporting Information Text) with two coupled oscillators ϕA and ϕB:

Christine
with a Gaussian signal function S: and a sine response function R: In the rectified pulsed-coupling model, we replace the response function R by an expression similar to Eq. 9 of the main text: In the piece-wise linear rectified pulsed-coupling model, we replace the response function R by the following function: where parameter γ plays the role of parameter β (Eq.26).Similarly, in the piece-wise linear rectified Kuramoto model, the H function of ODEs: becomes: where again, parameter γ plays the role of parameter β (Eq. 9 of the main text).

Methods
Mouse lines.The majority of the experiments were performed using a transgenic mouse line expressing a dynamic Notch signaling reporter controlled by the Lfng promoter and also have a constitutive expression of H2B-mCherry.This mouse line is named: H2BmCherry/LuVelu.The generation of the LuVelu mouse line was previously described in (5).
Medium preparation.Dissection and culture media were freshly prepared as indicated in Table S1: culture medium was filtered using a PVDF filter, pore size 0.22 µm (Merck) and kept in the incubator at 37 °C for at least 15 minutes to equilibrate it.
Randomization assay.To perform randomization experiments, we used the H2BmCherry/LuVelu mouse line.Mouse embryos were collected at 10.5 dpc (days post coitum) in a dissection medium containing DMEM/F-12 (Cell Culture Technologies) supplemented with 2 µM of Glucose, 2 µM L-Glutamine (Thermofisher 25030081), 0,5 % of BSA, HEPES and penicillin/Streptomycin (Thermofisher 15140122).Using a scalpel, PSM tissues were cut at the very posterior tip of the mouse tails to ensure a narrow phase and frequency distribution between cells.Randomization process included mechanical dissociation, cell filtration through a 30 µm filter (ParTec) and reaggregation of the PSM tissue in the micro insert 4 well Ful trac (Ibidi 80486) coated with fibronectin (Sigma-Aldrich F1141).Finally, PSM randomized cells were cultured overnight at 37 °C and 5 % CO2, in DMEM/F-12 medium supplemented with 2 µM of Glucose, 2 µM L-Glutamine and 0.5 % of BSA.
Phase synchronization assay.To study the phase synchronization dynamics, we used two very posterior PSM tissues from two different mouse embryos.Each of them was randomized as described above.One part of each population was kept as a population reference while the other part was mixed in the ratio of 1:1 with the other cell population.To distinguish the origin of two different cell populations, one population carried mCherry as reporter associated with H2B, while the other one didn't.
Video-microscopy.Imaging was performed using a Zeiss LSM780 laser-scanning microscope featuring an incubator chamber (EMBL Mechanical Workshop) for CO2 and temperature control.Samples were excited with an Argon laser at 514 nm for Venus and for mCherry.An 20X apo objective was used.Every 5 or 10 minutes a stack of 4 planes at 6 µm was scanned.
Multiple samples were recorded using a motorized stage during each experiment.Movies were recorded in 512 × 512 pixels, 1.38 µm per pixel. .Tail experiments and phase resetting.(A) The full PSM is placed in a dish and its oscillations are monitored (dark red), before a RAFL experiment is performed with cells from the PSM's tailbud (orange).The tail and RAFL oscillation phases, resp. in solid black and grey, are extracted using the pyBOAT software.To estimate the resetting due to the RAFL procedure, we compare the phase at the moment when the cells are placed in the incubator (black dot, ϕ tail ) to the RAFL phase two periods later (grey dot, ϕ RAF L ).Since we cannot image the cells as they are placed in the incubator, we estimate ϕ tail by computing the tail phase one period before incubation.A Ho, Laurent Jutras-Dubé, Michael Zhao, Gregor Mönke, Istvan Z. Kiss, Paul François, and Alexander Aulehla 3 of 19

Fig. S2 .SUPP
Fig. S2.Twin experiments.(A) Schematic of the twin experiments: two RAFLs are independently performed with cells coming from the same embryo.(B) Time series and polar plots for all twin experiments (n=8).The thick black line indicates the time of the associated polar plot.

Fig. S4 .Fig. S6 .SUPP. FIG. 7 CFig. S7 .Fig. S8 .CFig. S9 .
Fig. S4.All synchronization experiments (n=32).(A) Experiments with winner-takes-it-all dynamics (n=20).The reference population from embryo A is shown in orange, the reference population from embryo B is shown in blue, and the mixed population AB is shown in purple.The thick black line indicates the time of the associated polar plot.The proportion of cells from embryo A (resp.B) in the mix is indicated in orange (resp.blue) next to the polar plot.If no proportion is indicated, then the proportions are 50/50.(B) Time series and polar plots of experiments in which the reference populations are synchronized (n=11).(C) Only one experiment underwent phase averaging.Christine Ho, Laurent Jutras-Dubé, Michael Zhao, Gregor Mönke, Istvan Z. Kiss, Paul François, and Alexander Aulehla 9 of 19

Fig. S12 .
Fig. S12.Testing the impact of varied proportions of cells from each input population in the mixed population.(A) Simulation results of the ReKu H + model with varied proportions of cells from input population A. The blue lines show transient behavior that eventually become fully winner-takes-it-all, "ahead wins" synchronization at later times.(B) Representative time series and polar plots of the simulations shown in A. (C) Simulation results of the Kuramoto model with varied proportions of cells from input population A. (D) Representative time series and polar plots of the simulations shown in C.