In this talk, I will present a novel numerical method for a class of continuity equations with concentration-dependent mobilities, which arises widely in applications in biology and material sciences such as tumor growth, swarming dynamics, solid-state wetting/dewetting and thin film surfactant dynamics. This kind of nonlinear continuity models can be regarded as the gradient flows with respect to some Wasserstein-like transport distances. By leveraging the variational structure, along with the dynamical characterization of the Wasserstein-like transport distance, we construct a fully discrete scheme than ends up with a minimization problem with strictly convex objective function and linear constraint, which can be solved by primal dual operator splitting schemes and its accelerated versions. Our method has built-in positivity or bounds preserving, mass conservation, and entropy decreasing properties, and overcomes stability issue due to the strong nonlinearity and degeneracy. I will show a suite of simulation examples to demonstrate the effectiveness of our algorithm.