Understanding the propagation of earthquake-induced waves and hydromechanical response in the near-surface layers is crucial to predict eventual large deformations of the ground. Recent advances in the constitutive modelling of sands have successfully captured the element-level behavior under undrained cyclic loading. However, the system-level response of geosystems, such as foundations and slopes, has not been analysed in a large deformation framework. Besides, the non-linear solid-liquid interaction between solid grains and pore fluid is commonly neglected, preventing a reliable calculation of pore water pressure generation and depth of liquefaction hazard. The Material Point Method (MPM) is an advanced particle-based numerical method capable of modelling large-strains and multi-phase problems. The objective of this research is to establish an MPM framework capable of capturing earthquake site response using an advanced constitutive model in a system-level simulation. In this context, a non-zero kinematic and tied boundary condition with a particle domain relocation technique is used to reduce wave reflection on artificial boundaries. Three effective stress plasticity constitutive models are incorporated into the MPM framework to simulate the non-linear cyclic hysteresis of soil due to earthquake shaking. Namely, the sand-like behavior is simulated in a parametric manner using the Hypoplasticity model [1], Sanisand bounding surface plasticity [2], and Intergranular Strain Anisotropy [3]. The extensive constitutive model calibration performed by [4] is considered, and the corresponding parameters are used in the MPM model to simulate Karlsruhe fine sand. An undrained effective stress formulation is adopted whereby pore pressure generation is calculated based on the volumetric strain predicted by the constitutive model. First, the site response is assessed using the three constitutive models through a soil column simulation. Second, a large-scale model of a shallow foundation is adopted to investigate the system response in terms of deformation and pore pressure evolution. Finally, limitations and future advancements in the MPM numerical framework are highlighted to motivate the modelling of multi-phase interactions associated with earthquake-triggered liquefaction.