Many chemical and biological systems involve reacting species with vastly different numbers of molecules/agents. Hybrid simulations model such phenomena by combining discrete (e.g., agent-based) and continuous (e.g., partial differential equation- or PDE-based) descriptors of the dynamics of reactants with small and large numbers of molecules/agents, respectively. We present a stochastic hybrid algorithm to model a stage of the immune response to inflammation, during which leukocytes reach a pathogen via chemotaxis. While large numbers of chemoattractant molecules justify the use of a PDE-based model to describe the spatiotemporal evolution of its concentration, relatively small numbers of leukocytes and bacteria involved in the process undermine the veracity of their continuum treatment by masking the effects of stochasticity and have to be treated discretely. Motility and interactions between leukocytes and bacteria are modeled via random walk and a stochastic simulation algorithm, respectively. Since the latter assumes the reacting species to be well mixed, the discrete component of our hybrid algorithm deploys stochastic operator splitting, in which the sequence of the diffusion and reaction operations is determined autonomously during each simulation step. We conduct a series of numerical experiments to ascertain the accuracy and computational efficiency of our hybrid simulations and, then, to demonstrate the importance of randomness for predicting leukocyte migration and fate during the immune response to inflammation.