Mastication is a vital human function and uses an intricate coordination of muscle activation to break down food. Collection of detailed muscle activation patterns is complex and commonly only masseter and anterior temporalis muscle activation are recorded. Chewing is the orofacial task with the highest muscle forces, potentially leading to high temporomandibular joint (TMJ) loading. Increased TMJ loading is often associated with the onset and progression of temporomandibular disorders (TMD). Hence, studying TMJ mechanical stress during mastication is a central task. Current TMD self-management guidelines suggest eating small and soft pieces of food, but patient safety concerns inhibit in vivo investigations of TMJ biomechanics and currently no in silico model of muscle recruitment and TMJ biomechanics during chewing exists. For this purpose, we have developed a state-of-the-art in silico model, combining rigid body bones, finite element TMJ discs and line actuator muscles. To solve the problems regarding muscle activation measurement, we used a forward dynamics tracking approach, optimizing muscle activations driven by mandibular motion. We include a total of 256 different combinations of food bolus size, stiffness and position in our study and report kinematics, muscle activation patterns and TMJ disc von Mises stress. Computed mandibular kinematics agree well with previous measurements. The computed muscle activation pattern stayed stable over all simulations, with changes to the magnitude relative to stiffness and size of the bolus. Our biomedical simulation results agree with the clinical guidelines regarding bolus modifications as smaller and softer food boluses lead to less TMJ loading. The computed mechanical stress results help to strengthen the confidence in TMD self-management recommendations of eating soft and small pieces of food to reduce TMJ pain.