  The possibility of the appearance of $\Delta$ isobars in neutron star matter
and the so called {\it $\Delta$ puzzle} is studied in a modified quark meson
coupling model where the confining interaction for quarks inside a baryon is
represented by a phenomenological average potential in an equally mixed
scalar-vector harmonic form. The hadron-hadron interaction in nuclear matter is
then realized by introducing additional quark couplings to $\sigma$, $\omega$,
and $\rho$ mesons through mean-field approximations. The couplings of the
$\Delta$ to the meson fields are fixed from available constraints while the
hyperon couplings are fixed from the optical potential values. It is observed
that within the constraints of the mass of the precisely measured massive
pulsars, PSR J0348+0432 and PSR J1614-2230, neutron stars with a composition of
both $\Delta$ isobars and hyperons is possible. It is also observed that with
an increase in the vector coupling strength of the $\Delta$ isobars there is a
decrease in the radius of the neutron stars.
