This is actually an interesting problem in classical mechanics, dating back to Huygens. We'll work with the three variables you define in the question, namely $(x,\theta_1, \theta_2)$. Also we'll set your $m = l = g = 1$ for simplicity.
Kinetic Energy
The position vector of the first pendulum bob is
$$\mathbb{r}_1 = (x+\sin\theta_1, \cos\theta_1)$$
whence we deduce its kinetic energy to be
$$2T_1 = \dot{x}^2 + 2\dot{x}\dot{\theta_1}\cos\theta_1 + \dot{\theta_1}^2$$
We can similarly find the kinetic energy of the second bob.
Finally we must take into account the kinetic energy of the support with mass $M$.
$$2T_3 = M\dot{x}^2$$
It's interesting to keep $M\neq 0$ since different values of $M$ give different behaviour.
Potential Energy
We assume gravity acts on the bobs as usual, producing potential energy terms of form $\cos\theta$. We also assume an elastic potential $kx^2$ pulling the table back to equilibrium. Overall we have
$$V = kx^2 - \cos\theta_1-\cos\theta_2$$
Equations of Motion
One can easily write down the Lagrangian $L=T-V$ and from it deduce the equations of motion
$$\ddot{\theta}+\ddot{x}\cos\theta+\sin\theta-\dot{x}\dot{\theta}\sin{\theta} = 0$$
$$(M+2)\ddot{x}+\ddot{\theta_1}\cos\theta_1+\ddot{\theta_2}\cos\theta_2 - \dot{\theta_1}^2\sin\theta_1-\dot{\theta_2}^2\sin{\theta_2} + kx = 0$$
Interestingly when I put these into Mathematica, there was no synchronization! It turns out the missing ingredient is damping.
Damping
Intuitively the phase difference between the pendulums must drift in a periodic way in the absence of any dissipative effects. Indeed that's what you see when numerically solving the above equations with Mathematica.

Recall that damping is usually modelled as an additive term proportional to the velocity. Adding in such terms for $\theta_1$, $\theta_2$ and $x$ now does produce the desired synchronization behaviour. For my initial conditions and choice of constants we get antiphase locking.

Summary of the Physics
Momentum transfer through a connecting medium coupled with dissipative effects leads to synchronization.
Better Models
To fully model the video mentioned you'd need a forcing term from the escapement mechanism of the metronomes. You can read about such an approach here. See also this Wolfram demonstration and the papers it references.
Towards Chaos
Evidently this setup is nonlinear and so generically displays chaotic behaviour. The study of such systems is particularly important in chemistry and biology. Here is a good introduction.
If you want to play around with this behaviour yourself, here's my rudimentary Mathematica code. Try playing with the constants and initial conditions.
sol = NDSolve[{30 x''[t] + y''[t] Cos[y[t]] + z''[t] Cos[z[t]] -
y'[t]^2 Sin[y[t]] - z'[t]^2 Sin[z[t]] + 30 x[t] + 2 x'[t] == 0,
y''[t] + x''[t] Cos[y[t]] + Sin[y[t]] - x'[t] y'[t] Sin[y[t]] +
0.02 y'[t] == 0,
z''[t] + x''[t] Cos[z[t]] + Sin[z[t]] - x'[t] z'[t] Sin[z[t]] +
0.02 z'[t] == 0, x[0] == 0, x'[0] == 0, y[0] == Pi/10,
y'[0] == 0, z[0.5] == 1, z'[2] == 0}, {x, y, z}, {t, 0, 1000}]
Plot[{Evaluate[y[t] /. sol], Evaluate[z[t] /. sol]}, {t, 0, 250},
PlotRange -> All]