Domain decomposition methods (DDMs) enable the simulation of large-scale multiphysics systems by partitioning a computational domain into smaller, parallelizable subdomains. Classical formulations such as Schwarz, FETI, and mortar methods assume that all subdomains are governed by well-characterized physics and can be uniformly treated using traditional numerical schemes. However, many systems include regions where the underlying physics are computationally prohibitive to solve. We present a hybrid DDM framework that integrates physics-based solvers with machine learning (ML) surrogates to efficiently handle such regions. In the proposed approach, selected subdomains are replaced with neural surrogates trained to predict interior solutions from boundary or interface information. In contrast, the remaining subdomains are solved using standard numerical methods such as finite differences. To do the surrogate training in the absence of direct boundary data, we introduce a Gaussian process (GP)-based boundary synthesis technique. Smooth boundary functions sampled from GP priors serve as synthetic interface conditions that generate a physically plausible boundary. The coupling between physics-based and data-driven subdomains is achieved through an iterative interface exchange scheme that enforces Dirichlet–Neumann consistency, ensuring both accuracy and stability across interfaces. The framework is demonstrated on elliptic (Poisson), parabolic (heat), hyperbolic (wave), and mixed problems, spanning both steady and transient regimes. For time-dependent systems, interface data are exchanged at each time step to maintain continuity. Across different forcing conditions, the hybrid solver achieves accuracy comparable to that of full-physics DDMs while reducing computational cost. A detailed parametric analysis is conducted to investigate the impact of surrogate complexity and GP hyperparameters. These results demonstrate that hybrid domain decomposition informed by GP-augmented training provides an efficient approach for complex PDE systems where complete physical modeling is expensive.