In this study, we investigate the application of fractional calculus to the mathematical modeling of biological systems, focusing on fractional-order-in-time partial differential equations (FTPDEs). Fractional derivatives, especially those defined in the Caputo sense, provide a useful tool for modeling memory and hereditary characteristics, which are problems that are frequently faced with integer-order models. We use the Chebyshev spectral approach for spatial derivatives, which is known for its faster convergence rate, in conjunction with the L1 scheme for time-fractional derivatives because of its high accuracy and robustness in handling nonlocal effects. A detailed theoretical analysis, followed by a number of numerical experiments, is performed to confirmed the theoretical justification. Our simulation results show that our numerical technique significantly improves the convergence rates, effectively tackles computing difficulties, and provides a realistic simulation of biological population dynamics.