## 2 Materials and Methods
### Processing
The model with 12 segments, 19 degrees of freedom, and 92 musculotendon actuators [44] was used to estimate lower extremity muscle and joint forces. A residual reduction algorithm [22] was first used to resolve dynamic inconsistency between the kinematics of the measured GRF and the model, and muscle kinetics were then quantified using static optimization, as described by [PERSON] et al. [45]. From the above static optimization processes, the peak model muscle forces during the stance phase were extracted after being normalized to bodyweight (BW). Furthermore, the impulses (BW-ms) during the stance phases for the aforementioned muscles were also extracted using a trapezoidal function. In addition, the biceps femoris long-head, biceps femoris short-head, semitendinosus, and semimembranosus muscle forces were summed to create the total hamstring muscle force, while the rectus femoris, vastus lateralis, vastus medialis, and vastus intermedius forces were also summed to create the total quadriceps muscle force. The peak total hamstring and quadriceps forces, as well as their impulses during the stance phase, were also extracted. In accordance with [PERSON] et al. [23], the mean (%) contributions of each of the individual hamstring and quadriceps muscles to the total hamstring and total quadriceps forces during the stance phase were quantified. Finally, the mean co-contraction ratio was calculated by dividing the total hamstring force impulse by the total quadriceps impulse and peak co-contraction ratio, by dividing the total hamstring force at the instance of peak patellofemoral force by the total quadriceps force at the same time point. As muscle forces are the main determinant of joint compressive forces [46], following the static optimization process the peak compressive patellofemoral joint force (BW) was calculated via the joint reaction analysis function, using the muscle forces generated from the static optimization process as inputs [22]. Furthermore, peak patellofemoral stress, expressed as kilopascals per bodyweight (KPa/BW), was quantified by dividing the patellofemoral force by the joint contact area. Patellofemoral contact areas were obtained by fitting a polynomial curve to the sex-specific data of [PERSON] et al. [47], who estimated patellofemoral contact areas as a function of the knee flexion angle using magnetic resonance imaging. In addition, peak patellofemoral force (BW-ms) and stress (expressed as kilopascals per bodyweight per millisecond--KPa/BW-ms) impulses during the stance phase were extracted, also using a trapezoidal function. Finally, patellofemoral force/stress load rates (BW/s and KPa/BW/s, i.e., kilopascals per bodyweight per second) were extracted by obtaining the maximum increase in force/stress between adjacent data points using a first derivative function. In addition to the above, the patellofemoral force (BW) and patellofemoral stress (KPa/BW) experienced per mile were also quantified in order to explore the effect of step length on loading at this joint [48]. The number of steps needed to complete 1 mile was  
Figure 2: Definitions of footstrike modalities according to the strike index [39] (**a**) Rearfoot strike; (**b**) midfoot strike; and (**c**) forefoot strike. quantified using the step length (m), which was determined by taking the difference in the horizontal position of the foot's centre of mass between the right and left limbs at footstrike from each running trial [49]. To obtain the force/stress experienced per mile, impulse values for these measures outlined above were multiplied by the number of steps required to run a mile.