Some preliminary consideration:
In a uniform expanding universe the gravitational potential fulfills in atomic units the equation
$$\frac{GM}{R} \propto 1$$
with $c=1$ and $R$ is the reciprocal Hubble parameter and $M$ the total mass of the visible universe. Now the assumption of Brans was that not only the geometry defines gravitation, but also the mass. So the gravitational constant should be depending on the mass. When you write this equation like
$$ \phi= \frac{1}{G} \propto \frac{M}{R} \propto \sum_i \frac{m_i}{r_i}$$
you can interpret this as $G$ depending on the contribution of all matter to the inertial reaction. Here $G$ is assumed to be a function of $\phi$ so that the variation of both quantities is assumed to be equal. Then you get a roughly equal equation to that one above when assuming a inhomogeneous wave equation
$$ \square \phi \propto \rho$$
with a scalar mass density $\rho$. The concept of this kind of wave equation is that the quantity on the right hand side describes a so called source function. These source functions describes the effect of the sources, here the mass, on the medium carrying the waves. This approach arises from the implementation of Machs principle which leads to the assumption that not only the geometry of the Rieman manifold, but also the mass have an effect on the gravitation.
Outgoing from the Einstein Lagrangian
$$\mathcal{L} = R + \frac{16\pi G}{c^4} \mathcal{L}_{\mathcal{M}},$$
the approach is to replace the gravitational constant by the scalar field $\phi$. For that the Einstein Lagrangian is multiplied by $G^{-1}=\phi$:
\begin{align}
\mathcal{L}_{\text{JBD}} :&= \frac{\mathcal{L}}{G} = \frac{R}{G} + \frac{16\pi}{c^4} \mathcal{L}_{\mathcal{M}} \\
&= R \phi + \frac{16\pi}{c^4} \mathcal{L}_{\mathcal{M}}.
\end{align}
Because of the previous considerations we want a wave equation for $\phi$. So we add a term $\mathcal{L}_\phi \propto \phi_\mu \phi^\mu$:
$$\mathcal{L}_{\text{JBD}}=R \phi + \frac{16\pi}{c^4} \mathcal{L}_{\mathcal{M}} + \mathcal{L}_{\phi}(\phi, \phi_\mu).$$
Because of unit consistency you can show that $\mathcal{L}_\phi$ has to be second order in space-time and first order in $\phi$. The $\phi$ in the denominator of $\mathcal{L}_\phi$ arises, because in physics the coupling constants describing the strength of fundamental forces are chosen dimensionless. Because $\omega$ describes a coupling between the mass and the scalar field this quantity should also be dimensionless. With these assumptions you should be able to verify the desired equation. Keep in mind that $\omega$ is only a dimensionless constant without further meaning in this derivation.
