We extended the nutrient–phytoplankton–zooplankton model involving variable-order fractional differential operators of Liouville–Caputo, Caputo–Fabrizio and Atangana–Baleanu. Variable-order fractional operators permits model and describe accurately real world problems, for example, diffusion or spread of nutrients or species in different states. Particularly, we model the interaction of nutrient phytoplankton and its predator zooplankton. The variable-order fractional numerical scheme based on the fundamental theorem of fractional calculus and the Lagrange polynomial interpolation was consider. Numerical simulation results are provided for illustrating the effectiveness and applicability of the algorithm to solve variable-order fractional differential equations.