Functional heterogeneity of a skeletal muscle refers to its ability to produce a variety of force vectors through the localised, regional recruitment of motor units. Existing three-dimensional, macroscopic, continuum-mechanical muscle models (3D models) neglect motor unit anatomy and are thus ill-suited to investigate a muscle's functional heterogeneity. Here, we develop a method to generate and enrich 3D models with motor unit anatomy and activity information: First, a virtual innervation process is applied to group idealised computational fibres into motor units. Second, the discrete microstructure is homogenised to form continuous volumetric motor unit distributions. Third, motor unit activity is computed as a function of its firing times and twitch properties. As a proof of concept, we generated multiple masticatory models, each with different (masseter) motor unit distributions. The masticatory models were supplied with identical motor unit activity to simulate a sub-maximal static bite. The virtual innervation method captured physiological hallmark features of motor unit anatomy such as size, shape, position, and overlap. We found that motor unit anatomy significantly influenced bite force direction while having minimal impact on magnitude. Motor unit position played a more significant role than motor unit overlap. These results highlight the importance of modelling motor unit anatomy in understanding movement, particularly for muscles with complex internal architecture. This method provides an environment to investigate the interplay between motor control, skeletal muscle mechanics, and joint kinematics. Such investigations are relevant for the functional heterogeneity of healthy and pathological muscles.