Calculate the Laplacian of Kuramoto model

Question

I am trying to calculate the Laplacian of the Kuramoto model which describes the phase dynamics of a set of $N$ phase oscillators $i$ with natural frequencies $\omega_i$ $$ \dot{\theta}_i = \omega_i + \frac{K}{N}\sum_{j=1}^N \sin(\theta_j - \theta_i). $$ The idea of the Laplacian as I understand it is that it would satisfy $\dot{\theta} = -L \theta$ where $\theta$ is a vector. I found a solution in this paper where the equation is first rewritten using order parameters of $Re^{i\psi} = \frac{1}{N}\sum_{j=1}^{N}e^{i \theta_j}$ in the usual way $$ \dot{\theta}_i = \omega_i - \frac{K}{N} R \sin(\theta_i - \psi). $$ (Actually it seems to me that they should replace $\frac{K}{N}$ with just $K$). From here they say that the Laplacian is given by $L = B B^T$ where $B$ is the incidence matrix, and they then claim that the incidence matrix satisfies $$ \dot{\theta} = \omega - \frac{K}{N} B \sin(B^T \theta) $$ where $\theta$ and $\omega$ are vectors. I see the analogy between these last two equations, but I don't understand at all how they have derived the incidence matrix or even what the incidence matrix would look like if it were written out in full. Does anyone have any ideas? Are there maybe also other ways to derive the Laplacian?