Three- and four-jet final states have been measured in photoproduction at HERA, using the ZEUS detector and an integrated luminosity of 121 pb$-1$.  
Both the three- and four-jet events have been studied in a semi-inclusive,
$M_{nj}>25$ GeV, and high-mass region, $M_{nj}>50$ GeV.  The three-jet cross
sections have been compared to an $O(\alpha\alpha_s^2)$ perturbative 
QCD calculation and are shown to be sensitive to higher-order processes at
low $M_{3j}$. The calculation describes the higher $M_{3j}$ region reasonably well.
In addition, the three- and four-jet cross sections have been compared with
leading-logarithmic Monte Carlo models that rely on parton showers to generate
events with high jet-multiplicities. The Monte Carlo models describe the data
well once multi-parton interactions are included. Results on underlying events
  are also presented in terms of the energy flow between jets and energy correlations 
  in an inclusive photoproduction sample. Finally, leading neutron data measured at
HERA in photoproduction and deep inlelastic scattering results are presented and
compared to different models on pion-exchange followed by secondary interactions.