# The two-community Kuramoto model - Mathematical approach

As written in the main article we can represent the interactions between the neurons in the body-clock as a network. In our case the network has two communities, where each community has many edges between the nodes, while there are only few edges between the nodes in different communities. The network describes which neurons interact with each other, but it does not yet tell us how they respond to one another.

This is where the Kuramoto model comes in. The state of each neuron at each point in time $t$ can be thought of as the position of a dot on a circle, or equivalently as an angle $\theta_{i, j}(t)$ between $0$ and $2\pi$. Here the subscript $i$ tells us which community the neuron is in (so $i=1$ for community $1$ or $i=2$ for community $2$), and the subscript $j$ tells us which neuron in that community it is, where we assume that there are $N$ neurons in each community, (so $j=1, \ldots, N$). For the interaction we want that a neuron adjusts the time at which it fires as a result of some other neuron firing before it, either by firing a little earlier, or by firing a little later. The easiest way to do this is to adjust the neuron represented by angle $\theta_{i, j}(t)$ as a response to the neuron represented by the angle $\theta_{k, l}(t)$ via the sine of the difference between the angles, i.e., $\sin(\theta_{i, j}(t)-\theta_{k, l}(t))$. Since a neuron is responding to all the neurons it is connected to simultaneously, it will adjust its angle by an amount that equals the sum of all the adjustments it would make for each neuron separately. In order to model the fact that the neurons are positioned in a fluid, we add a small random term to the total adjustment.

### Kuramoto's model of synchrony

The Kuramoto model has been widely studied since it was conceived in 1975. Most studies have focused on the case where all angles interact with all other angles (i.e., the network is on the complete graph). Since this is clearly not reflecting the structure of the body-clock, we analyze what happens if we consider the Kuramoto model on a two-community network. For this purpose, we take the neurons in the same community to interact with a strength given by a parameter $K>0$ and the neurons in different communities to interact with a strength given by parameter $L$, which can be positive, zero or negative. When $L$ is positive the two communities try to align with each other, when $L$ is zero there is no interaction and when $L$ is negative the two communities repel each other. We are interested in studying which states are possible for the system as a whole when we fix $K$ and $L$ and consider the number of neurons in each community to be very large.

What may we expect? In the Kuramoto model with the all-to-all interaction i.e., all neurons interact with all other neurons with strength $K$, synchronization appears as a threshold phenomenon. This means that for values of $K$ below a critical value $K_{c}$ the neurons fire completely randomly, while for values of $K$ above $K_{c}$ the neurons start to fire in a synchronized way. When $K$ exceeds $K_{c}$, synchronization is measured by what is called an order parameter. The order parameter for the synchronization is $r(t)$, which takes values between $0$ and $1$, where $0$ represents the case where the neurons are completely unsynchronized and $1$ represents the case where the neurons are completely synchronized. We may therefore expect that synchronization in the two-community Kuramoto model will also appear as a threshold phenomenon.

If the environment in which the body-clock resides remains constant, then we may expect that the system will reach a stable and stationary synchronization level after some time. This means that the synchronization level will remain the same as long as the environment does not change.

Of course, in the two-community case we do not only want to know the level of synchronization in the body-clock, but also the synchronization in each community. In order to characterize all the possible stationary states, we also want to know the difference between the average angles in the two communities, which we call the phase difference and denote by $\psi$. We found that the phase difference $\psi$ can only be $0$ or $\pi$ in the stationary state. In both of these cases, we can split the possible stationary solutions into three categories. In category (I) both communities are unsynchronized; in category (II) both communities are synchronized to the same degree (symmetrically synchronized); in category (III) both communities are synchronized, but one community is more synchronized than the other (nonsymmetrically synchronized).

The two-community version of the model has two parameters, $K$ and $L$, which means that we do not have a single critical value, like in the one-community version, but have a critical line. The critical line, splitting the unsynchronized state (category (I)) from the symmetrically synchronized state (category (II)), is $K+L=2$. The nonsymmetrically synchronized state (category (III)) is only possible in a part of the symmetrically synchronized section of the plane and requires $L<0$ when $\psi=0$ and $L>0$ when $\psi=\pi$. The critical line dividing the symmetrically synchronized state and the nonsymmetrically synchronized state is identified as the solution to an equation involving special functions, which means that we cannot write down the solution explicitly. In order to visualize the situation, we plot the critical lines in Figure 1 below and indicate in which regions which solutions are possible. Figure 1: In the light red region there is one pair of solutions: unsynchronized (U). In the light green region there are two pairs of solutions: unsynchronized and symmetric synchronized (S). In the light blue region there are three pairs of solutions: unsynchronized, symmetric synchronized and non-symmetric synchronized (NS). The left image is for the case $\psi=0$ and the right image is for the case $\psi=\pi$.

### Simulations

The previous section dealt with the possible stationary states. In the experiments with rats and mice one of the most interesting features that was observed was the possibility to see two different types of transitions to the phase-split state, which is certainly not something that can be understood by studying stationary states. In our model we did, however, see that stationary states exist in which the two communities have a phase difference of $\pi$, which corresponds to the phase-split state observed in experiments. To study the transitions between these states (the one with $\psi=0$ and the one with $\psi=\pi$) is a mathematically intractable problem at the moment, but it is possible to see what sort of transitions are possible. The possible transitions that we observed are shown in Figures 2 and 3 below. Here, the average angle for community 1 and community 2 are given by $\psi_{1}$ and $\psi_{2}$, and the synchronization levels are given by $r_{1}$ and $r_{2}$. Figure 2: Simulation of 10000 oscillators per community with $L=-2$ (A) or $L=2$ (B) and $K=5$. The time step is set at ${\rm d}t=0.01$. The top images show the synchronization levels ($r_{1}$ and $r_{2}$), the bottom images show the average phases ($\psi_{1}$ and $\psi_{2}$). The phases are plotted in such a way that the different communities can be seen as the two components of activity in the suprachiasmatic nucleus. The time on the vertical axis can be regarded as consecutive days. (A) and (B) show the case where $r_{1}$ and $r_{2}$ are the same. Figure 3: Simulation of 10000 oscillators per community with $L=-2$ (C) and $L=2$ (D) and $K=7$. The time step is set at ${\rm d}t=0.01$. The top images show the synchronization levels ($r_{1}$ and $r_{2}$), the bottom images show the average phases ($\psi_{1}$ and $\psi_{2}$). The phases are plotted in such a way that the different communities can be seen as the two components of activity in the suprachiasmatic nucleus. The time on the vertical axis can be regarded as consecutive days. (C) and (D) show the case where $r_{1}$ and $r_{2}$ differ at first, but both approach 1 when the stable state is reached.

### Interpretation

The simulations also show two different transitions between the symmetrically synchronized state with $\psi=0$ and the symmetrically synchronized state with $\psi=\pi$, and back. In the first type of transition, both communities move to a higher level of synchrony together while their average phases either move further apart or move closer together. In the second type of transition, one community is forced to a much lower level of synchrony before it moves to the higher level of synchrony, while the other community moves directly to the higher level of synchrony. The first situation could possibly correspond to the smooth transition observed in the experiments, while the second transition could correspond to the chaotic transition.

There are two factors that seem to play a role in determining which transition occurs. The first is the initialization of the system and the second is the existence of the nonsymmetrically stationary state. When we start the simulation, we have to chose an initial state for the angles of the neurons. In experiments this state is difficult to interpret, but we can think of it as the situation that prevails immediately after a change of the environment in which the animals live has occurred (e.g. switch from light-dark cycle to constant light). The second factor is more easy to interpret: the strength of the interactions in the network is given by $K$ and $L$, and if these are such that the nonsymmetrically synchronized state is possible, then the system might be forced to the nonsymmetrically synchronized state before it can move to the phase-split state.

Proceed to the Animation of the two-community Kuramoto model.