At least in the Kuramoto model, it's a little more complicated than the OP makes it sound (but not much). Step 1 is more like "if you see a flash, bump your clock closer to midnight."
Specifically, an ODE is given for firefly $i$'s phase angle $
\theta_i$
$$\dot\theta_i = \omega_i + K / N \sum_{j=1}^N \sin(\theta_j - \theta_i)$$
This gets interesting because the flies natural frequencies $\omega_i$ are also assumed to be randomly distributed. So you don't get perfect phase-synchronization--flies with fast natural frequencies lead the pack as it loops around the phase ring, and flies with slow natural frequencies are dragged along at the back. Relative to the mean phase, your excess phase approaches a smooth increasing function of your natural frequency. But, for high enough values of $K$, you do get frequency-synchronization--everyone oscillates at the average frequency.
(For low values of $K$, or too-large natural frequencies, "rogue oscillators" emerge in a SNIPER bifurcation. They zoom around the phase ring at a different frequency, briefly slowing as they pass through the cloud of synchronized oscillators. Also applies for too-slow rogues.)
In the video and the OP's simulation, it looked the natural frequencies were all pretty similar, if not the same. However, there was a second addition that is not in the original Kuramoto model (but is in most subsequent models): rather than observing all other flies, only observe nearest neighbors. This can be added by putting a symmetric boolean adjacency matrix $A_{i,j}$ right before the $\sin$ in the previous equation.
This has the effect of making excess phase a smooth function of not only the natural frequency, but also some structural feature imposed by the network. In random networks like Erdos-Renyi, this feature is the node degree, but in the video it looks like it might be the long axis of the bush (so, like, one of the eigenvectors of the graph Laplacian).
(In general, coupled oscillators, such as circadian gene clocks, show this smooth dependence of excess phase on per-unit heterogeneities, which is the topic for the first half of my PhD dissertation. The second half is figuring out what the heterogeneities are when you only have recordings of the dynamics to go by.)