In this paper, we present the first implementation of a deep neural network-based variational Monte Carlo method for multi-quark bound states. To address the unique challenges of multi-quark systems (strong correlations, additional discrete quantum numbers, and intractable isolation interactions) that are more complex than electron or nucleon systems due to strong SU(3) color interactions, we design a novel, highly efficient architecture, DeepQuark. DeepQuark demonstrates competitive performance with state-of-the-art approaches, including diffusion Monte Carlo and Gaussian expansion methods, in the nucleon, doubly baryon quadruple quark, and fully baryon quadruple quark systems, and outperforms existing calculations in particular for triple-baryon quintuplet quarks. For nucleons, we successfully incorporate the three-body flux-tube isolation interactions without additional computational overhead. In the quadruple quark system, we consistently describe the hadron molecule $T_{cc}$ and the compact quadruple quark $T_{bb}$ in the unbiased form of the wavefunction ansatz. In the field of quintuple quarks, we have obtained the weakly bound $\bar D^*\Xi_{cc}^*$ molecule $P_{cc\bar c}(5715)$ and its ground partner $P_{bb\bar b}(15569)$ with $S=\frac{5}{2}$. They can be regarded as analogues of the molecule $T_{cc}$. We encourage experimental exploration of $P_{cc\bar c}(5715)$ in the D-wave $J/\psi \Lambda_c$ channel. DeepQuark holds great promise to overcome the computational barriers of existing methods and to extend the reach to larger multi-quark systems, serving as a powerful framework to explore isolation mechanisms beyond two-body interactions in multi-quark states, which may provide valuable insights into nonperturbative QCD and general many-body physics.
