Chemotaxis is a biological phenomenon whereby unicellular organisms direct their movements in response to certain chemicals in their habitat. This study presents a numerical investigation of the fractional Keller–Segel model describing the aggregation of cellular slime molds and bacterial chemotaxis. Two numerical schemes are provided to solve this model; primarily, a meshfree numerical scheme based on the local radial basis function partition of unity method is presented. In this approach, the domain is split up into a number of smaller, overlapping subdomains, and the radial basis function interpolation is performed separately on each of these. On the other hand, a numerical method employing the L1 scheme for temporal discretization and centered difference for spatial discretization is introduced to compare the primary proposed method solutions with the simulations acquired by this method. Stability and convergence of the time-discrete algorithm are rigorously established. The strengths of the carried work is that the proposed approach is meshfree, where as the classical methods like finite difference/element approaches depends on mesh. Also, as per the best of authors knowledge, the analytical solutions of the considered fractional model are not known in literature, which makes the carried numerical investigation innovative. Computational experiments are carried out, and simulation results of both schemes are compared. Also, the density plots of cellular slime mold and the chemical attractant for specific biological parameters are illustrated to observe their biological behavior.