MOTIVATION: Single-cell time series data often exhibit significant variability within an isogenic cell population. When modeling intracellular processes, it is therefore more appropriate to infer parameter distributions that reflect this variability, rather than fitting the population average to obtain a single point estimate. The Global Two-Stage (GTS) approach for nonlinear mixed-effects (NLME) models is a simple and modular method commonly used for this purpose. However, this method is computationally intensive due to its repeated use of nonconvex optimization and numerical integration of the underlying system. RESULTS: We propose the Gradient Matching GTS (GMGTS) method as an efficient alternative to GTS. Gradient matching offers an integration-free approach to parameter estimation that is particularly powerful for systems that are linear in the unknown parameters, such as biochemical networks modeled by mass action kinetics. By incorporating gradient matching into the GTS framework, we expand its capabilities through uncertainty propagation calculations and an iterative estimation scheme for partially observed systems. Comparisons between GMGTS and GTS across various inference setups show that our method significantly reduces computational demands, facilitating the application of complex NLME models in systems biology. AVAILABILITY AND IMPLEMENTATION: A Matlab implementation of GMGTS is provided at https://github.com/yulanvanoppen/GMGTS (DOI: http://doi.org/10.5281/zenodo.14884457).