The mechanical properties of vertebrate bone are largely determined by a process which involves the complex interplay of three different cell types. This process is called bone remodeling and occurs asynchronously at multiple sites in the mature skeleton. The cells involved are bone resorbing osteoclasts, bone matrix producing osteoblasts, and mechanosensing osteocytes. These cells communicate with each other by means of autocrine and paracrine signaling factors and operate in complex entities, the so-called bone multicellular units (BMUs). To investigate the BMU dynamics in silico, we develop a novel mathematical model resulting in a system of nonlinear partial differential equations (PDEs) with time delays. The model describes the osteoblast and osteoclast populations together with the dynamics of the key messenger molecule RANKL and its decoy receptor OPG. Scaling theory is used to address parameter sensitivity and predict the emergence of pathological remodeling regimes. The model is studied numerically in one and two space dimensions using finite difference schemes in space and explicit delay equation solvers in time. The computational results are in agreement with in vivo observations and provide new insights into the role of the RANKL/OPG pathway in the spatial regulation of bone remodeling. © 2010 Society for Industrial and Applied Mathematics.