In MMAS-based biological system simulation, it is a challenging task to deal with numerous interactions among a vast number of autonomous agents. In our work, a hybrid massively multi-agent systems (MMAS) model is developed, and it incorporates the characteristics of cellular automaton (CA) and system-level mathematical equation modeling to simulate HIV-immune interaction dynamics. The mathematical equations are adopted within the site of a two-dimensional lattice. As the average high density, agent interactions can be calculated according to the equations without significantly affecting the performance of the systems studied. In the mean time, the CA model keeps the spatial characteristics of HIV evolution among the sites. The simulation based on the implemented MMAS discovers the dynamics of HIV evolution over different temporal and spatial scales, and reproduces the typical three-stage dynamics of HIV infection.