The work presents the qualitative analysis of the free boundary value problem related to the biosorption process in multispecies biofilms. In the framework of continuum biofilm modeling, the mathematical problem consists of a system of nonlinear hyperbolic partial differential equations for microbial species growth and spreading, a system of semilinear parabolic partial differential equations describing the substrate trends and a system of semilinear parabolic partial differential equations accounting for the diffusion, reaction and biosorption of different agents on the various biofilm constituents. Two systems of nonlinear hyperbolic partial differential equations have been considered as well for modeling the dynamics of the free and bounded sorption sites. The free boundary evolution is regulated by a nonlinear ordinary differential equation. Overall, this leads to a free boundary value problem essentially hyperbolic. The main result is the existence and uniqueness of the solutions to the stated free boundary value problem, which have been derived by converting the partial differential equations to Volterra integral equations and then using the fixed point theorem. Moreover, the work is completed with numerical simulations for a real case examining the growth of a heterotrophic–autotrophic biofilm devoted to wastewater treatment and acting as a sorbing material for heavy metal biosorption.