Do non-rearfoot runners experience greater second metatarsal stresses than rearfoot runners?

Stress fracture of the second metatarsal is a common and problematic injury for runners. The choice of foot strike pattern is known to affect external kinetics and kinematics but its effect on internal loading of the metatarsals is not well understood. Models of various complexities can be used to investigate the effects of running characteristics on metatarsal stresses. This study aimed to compare second metatarsal stress between habitual rearfoot and non-rearfoot strikers during barefoot running, using a novel participant-specific finite element model, including accurate metatarsal and soft tissue geometry. Synchronised force and kinematic data were collected during barefoot overground running from 20 participants (12 rearfoot strikers). Stresses were calculated using a previously evaluated and published 3D finite element model. Non-rearfoot strikers demonstrated greater external loading and joint contact forces than rearfoot runners, but there were no differences in stresses between groups. Additionally, the study allowed for a qualitative assessment of bone geometries and stresses. No correlation was found between bone volume and stresses, however, there was found to be a large variation in metatarsal shapes, possibly accounting for the lack of difference in stresses. This emphasises the importance of bone geometry when estimating bone stress and supports the suggestion that external forces should not be assumed to be representative of internal loading.


Introduction
Stress fracture to the second metatarsal is a common and problematic injury amongst runners (Bennell et al., 1996;Iwamoto and Takeda, 2003;Milgrom et al., 1985;Chuckpaiwong et al., 2007). A primary driver of stress fracture development is thought to be the magnitude of stress that the bone is subjected to, either directly via the accumulation of microdamage, or via an intermediate bone remodelling process (Milgrom et al., 2002;Schaffler and Jepsen, 2000;Carter et al., 1981). Measurement of stresses in the second metatarsal bone can be conducted directly, using bone staple strain gauges Milgrom et al., 2002), but the invasive nature of these studies limits their use when answering applied questions. Recent research has suggested that using external measurements, such as ground reaction forces, as a direct estimate of internal forces may not be valid (Matijevich et al., 2019;Ellison et al., 2020b). Mathematical modelling has been shown to be a valuable alternative to direct measurement of bone stresses on the metatarsal bones in the past (Gross and Bunch, 1989;Fung et al., 2017;Nunns et al., 2017) and these have ranged from beam theory models with limited participant-specific geometry (Gross and Bunch, 1989) to models incorporating complete participant-specific geometry and using finite element methods to estimate highly realistic bone deformations (Firminger et al., 2017;Li et al., 2017). It has been suggested that bone geometry is an important determinant of stress magnitude (Nunns et al., 2017). Evaluation of a three-dimensional finite element model of the second metatarsal, incorporating participant-specific bone geometry, was recently evaluated and deemed appropriate for the estimation of peak second metatarsal stresses on the dorsal and plantar surfaces of the bone during running (Ellison et al., 2020a).
There is some evidence that changing to run in so-called minimalist footwear may affect the risk of developing a stress fracture (Cauthon et al., 2013;Salzler et al., 2012), possibly due to the alteration of foot mechanics seen in runners who move from a cushioned shoe to a minimalist shoe (Firminger and Edwards, 2016). However, there is currently limited understanding of whether this is due to the change in cushioning, a change in landing mechanics or unaccustomed activity. If the change to minimalist footwear results in a more anterior foot-strike pattern, this could result in an increased magnitude of stress, leading to rapid accumulation of microdamage. Recently, a study investigating stresses in the second metatarsal found no statistically significant difference in peak stress between runners who landed on their heels and runners who did not (Ellison et al., 2020b) despite differences in ground reaction forces. This study utilised a beam theory model to determine stresses, which has certain limitations. For example, there is no ability to account for the cushioning soft tissue effects between the metatarsal and ground. Additionally, the beam theory model did not account for deformation of the metatarsal, which changes both the moment arm and inertial properties of the bone. A further study (Li et al., 2017) used a whole foot model to investigate stress differences between forefoot and rearfoot strike, finding increased stresses at mid stance with a forefoot strike. However, this was a theoretical study of one participant and was not able to establish whether such differences exist between habitual runners in practice.
This study aimed to use a participant-specific finite element model (Ellison et al., 2020a) to investigate second metatarsal stresses between different habitual foot-strike groups and to investigate whether there is any correlation between bone volume and peak stress. In addition, the data set used here has previously been presented using a participantspecific 2D beam theory model (Ellison et al., 2020b) and therefore this permits brief comparison of both modelling strategies.
It was hypothesised that there would be no differences in the magnitude of peak stresses on the second metatarsal between groups, but that the input forces would be higher in those who run with a more anterior foot strike compared with rearfoot strikers. Further it was expected that bone volume would have an inverse correlation to peak stress.

Methods
The data collection methods have been previously published (Ellison et al., 2020b;Ellison et al., 2020a), but are briefly described here for clarity.

Participants
20 injury free athletes participating in running activity (>3 times per week and >150 mins per week) were recruited. Sample size was estimated using G*Power (Faul et al., 2007) as detailed in Ellison et al. (2020b). No participants reported any current injuries affecting their running regimen and no participants had sustained any lower limb injuries that prevented their normal training within the last year. Eligible volunteer participants provided written informed consent. The study was approved by the Sport and Health Sciences Research Ethics Committee, University of Exeter (Ellison et al., 2020b). When analysing data, all 20 participants were included in the 2D beam theory model, however, data from one rearfoot participant was found to be unsuitable for use in the finite element analysis, therefore the participant numbers for finite element analysis are described in Table 1:

Running protocol
Synchronised kinematic, kinetic and plantar pressure data were collected during barefoot running at 3.6 ms − 1 using four CX1 units (200 Hz) (Codamotion, Charnwood Dynamics Ltd., U.K.) with an integrated force plate (1000 Hz, AMTI, MA, USA). Pressure data were collected at 200 Hz via Footscan software (RSscan Gait v7) and a plantar pressure plate (RSscan 0.5 m Hi-End Footscan, Belgium) placed within the boundaries of the force plate. An EVA runway was laid so that the contact surface was representative of a shoe midsole. 19 markers representing bony landmarks of the foot and shank were captured. Foot strike was assessed and participants categorised according to strike index (Cavanagh and Lafortune, 1980) as either rearfoot (RF) or nonrearfoot (NRF) runners. This assessment was conducted both shod (in participant's own running shoes) and barefoot for comparison (see Ellison et al. (2020b) for full details). Only two participants displayed different foot strikes between footwear conditions, and these were categorised as NRF runners based on their barefoot foot strike, as all further analyses were from barefoot trials.
Following a warm up and familiarisation, the experimental protocol involving ten barefoot running trials running at 3.6 m.s − 1 (±5%) with a right foot pressure plate contact was conducted, velocity was monitored using light gates (WITTY system, Microgate, Bolzano, Italy). The magnetic resonance (MR) imaging protocol is detailed in Ellison et al. (2020b).

Finite element analysis
The finite element model used in the present study was previously outlined in detail and evaluated (Ellison et al., 2020a). This model consisted of a second metatarsal mesh constructed from participantspecific MR images, with cortical and trabecular bone differentiated, and embedded within an encapsulating soft tissue (Fig. 1). The bone was fixed at the proximal end and positioned according to the sagittal and frontal plane metatarsal-ground vector angle obtained from kinematic data. The bone was loaded via a simulated floor using scaled vertical and anterior-posterior ground reaction forces and using joint contact forces at the distal end. The joint contact forces represented the forces from the floor transferred to the metatarsal via the toe, combined with the action of the plantar tendons (Ellison et al., 2020b). The model was simulated at three points of stance for each participant. Three time points of interest during the stance phase were chosen based on when maximal  horizontal or vertical forces occurred. These represent the time of maximum braking force (maximum negative horizontal ground reaction force); time of maximum vertical ground reaction force; and time of maximum propulsive force (maximum positive horizontal ground reaction force). These time points will henceforth be referred to as early, mid and late stance, respectively. Model simulation for each time point took approximately 35 min, after which the peak von Mises stress on the cortical bone was obtained for each participant at each of the three stance phases. Detailed model development, including segmentation protocol, selection of material properties, mesh sizing, loading and boundary conditions have been previously described and evaluated (Ellison et al., 2020a).

Analysis
Peak von Mises second metatarsal stresses during running were compared between foot strike groups, in addition to horizontal and vertical ground reaction forces, and horizontal and vertical joint contact forces. Statistical analyses were carried out in SPSS (Version 24.0, SPSS Inc., Chicago, IL, USA) with a significance level of P ≤ 0.05. Variables were examined using a Shapiro-Wilk test to confirm normality (P ≥ 0.05). Means were compared using an independent T-Test (normally distributed variables) or a Mann-Whitney-U test (non-normally distributed variables). Effect sizes were calculated using Cohen's d (Cohen, 1988) for all variables where P ≤ 0.05. Additionally, a qualitative analysis of bone shape and the location of peak stresses is given. The location information stated whether maximum stress occurred on the dorsal or plantar aspect of the bone and also the location along the shaft, where proximal shaft, mid-shaft and distal shaft refer to the approximated thirds of the length.
The volumes of the cortical, trabecular and total bone were correlated with the peak stress using a Pearson's correlation coefficient, to establish the importance of participant-specific geometry when estimating bone stress.

Discrete analysis
There was no difference in maximum stress values between groups in early stance, mid stance or late stance (Table 2; Fig. 2). Vertical ground reaction forces were significantly greater in the non-rearfoot runners than the rearfoot runners at all three time points. Horizontal ground reaction forces were significantly greater in the non-rearfoot group in early and late stance but were similar in mid stance. Both vertical and horizontal joint contact forces were greater in the non-rearfoot group in early and mid stance but were similar in late stance.

Qualitative assessment of metatarsal stresses
There was a large range of maximum second metatarsal stresses across participants during running (28.12 to 79.02 MPa). Group mean peak stresses were 19.81 MPa, 44.99 MPa and 33.54 MPa at early, mid and late stance, respectively. Maximum stress at mid stance was located within the diaphysis in all cases (e.g. Fig. 3). For each participant, at all three simulated time points, the maximum stress magnitude at the dorsal surface was similar to the maximum stress magnitudes at the plantar surface. Importantly, there were differences in both the shapes and proportions of the metatarsals in the present study, with some much more slender relative to their length than others (Fig. 4).
Maximum stresses were more likely to be located at approximately mid-shaft with 64%, 70% and 82% of participants experiencing the maximum stress in the mid-shaft during early, mid and late stance, respectively. No participants displayed maximum stresses in the distal shaft at mid stance, which was the time of overall maximum stress of the three time points analysed in all but one participant. Bold text/* Significant (P ≤ 0.05) between groups, negative stress values indicate compression, GRF = Ground Reaction Force, JCF = Joint Contact Force. Table 3 shows the average metatarsal bone volume of each group, with a further breakdown by cortical and trabecular bone types. Examination using a Pearson's correlation coefficient reveals that total bone stress was not associated with total bone volume (r = 0.101, P = 0.699) nor with trabecular or cortical bone volume (Cortical vs. Stress: r = − 0.139, P = 0.596; Trabecular vs. Stress: r = 0.263, P = 0.308).

Discussion
This study compared peak second metatarsal stresses between rearfoot and non-rearfoot runners at three points of stance using a novel, finite element model of the second metatarsal (Ellison et al., 2020a). It also provided a qualitative analysis of bone geometries and locations of maximum stresses.
Regarding the discrete analysis of stresses between groups, the results here support the results from the beam theory model, finding that whilst input forces were significantly greater throughout stance in nonrearfoot runners compared to rearfoot runners, the peak stresses were not significantly different between groups. However, at mid and late stance, there were non-significantly higher peak stresses in the nonrearfoot than rearfoot runners with P values < 0.06. In combination with the large effect size of 1.0, this suggests a trend of greater stresses for NRF runners than RF runners. Peak stresses show greater relative variability than peak external forces, which partly reflects the complex interplay between external forces, bone geometry and bone orientation. On balance, the current results suggest there may be a tendency towards higher stresses in NRF runners than RF runners, but not to the extent indicated by the external forces. Therefore, it should not be concluded that foot strike per se influences the magnitude of peak metatarsal stress during running, and therefore may not influence the risk of stress fracture via this mechanism. These results also support work by Matijevich et al. (2019) who found that greater ground reaction force metrics did not correlate to greater tibial forces, and suggests that externally measured forces may not be a valid method for inferring internal loading, and therefore injury risk.
Bone geometry has previously been noted as an important determinant of bone stress (Nunns et al., 2017) and this may partly explain why differences in external forces between groups did not equate to differences in internal forces. The fact that no association between total volume of the bone and maximum stress was observed in the present study suggests that some aspect of the shape rather than size per se may influence stress magnitude. However, there are numerous ways to quantify bone geometry and the geometry may influence peak stress differently according to the phase of stance. There is evidence that bone geometry may adapt over time through habitual loading, for example elite runners show thickened cortical bone and reduced trabecular bone compared to recreational runners (Hart et al., 2017). However, it is not known what magnitude of force or number of cycles induces these levels of adaptation. All the participants in the present study were recreational runners, reporting similar levels of activity, so it was assumed that the level of adaptation between participants would be equal unless influenced by their habitual foot strike. It is also possible that the geometry of the soft tissues may mitigate the forces more in some participants than others. For example if the distance between the most plantar aspects of  the metatarsal head and the sole of the foot is greater (i.e. soft tissue is thicker) it can be assumed there would be more cushioning and that forces may be damped more than with less soft tissue. The stress distribution on the metatarsal bones in the present study for all the participants in this study is visually similar to that seen in research by Firminger et al. (2017), with the stresses concentrated at the mid-shaft in the dorsal and plantar regions. The medial and lateral regions of each metatarsal displayed very low stresses, but it should be noted that this pattern of stress distribution is influenced by the modelling of the metatarsal as a cantilever with no mediolateral or torsional forces.
The shape of the metatarsals highlights the variability of individual geometry and is one potential explanation for the lack of statistical difference between groups, despite there being no correlation between bone volume and stress. For example, some bones are slender and long, whilst others are wider and shorter. This difference would result in potentially similar volumes, but two very different mechanical responses to the input forces.
The magnitudes of peak stresses on the plantar and dorsal surfaces for each simulation were similar meaning that identifying whether the maximum stress occurred on the plantar or dorsal surface is of relatively low importance, whereas it may be important to identify the proximaldistal location of peak stress. However, it should be noted that whilst the magnitudes were similar between the dorsal and plantar surfaces, the dorsal surface is predominantly under compression throughout stance, whilst the plantar surface is under tension. Bone is known to have different mechanical responses to compression vs. tension (Martin et al., 2015), being stronger in compression. Therefore, a given magnitude of compression on the dorsal surface may not induce either the same level of microdamage, or the same level of remodelling response as the same magnitude of tension on the plantar surface, and stress fracture development may be more likely on the plantar surface.

Comparison of modelling methods
The data presented in this study have previously been used in a 2D beam theory model (Ellison et al., 2020b) to compare metatarsal stresses between RF and NRF runners. Interestingly, the non-significant differences in peak stresses between groups were more marked using the finite element model in the present study than when comparing these two groups using the 2D beam theory approach. This suggests the present model is more sensitive which may be important for clinical research.
The finite element method has distinct advantages when considering the quality and nature of the output data. The model considers all points along the metatarsal and stress is computed at every node. This allows the peak stress and its location to be identified and permits qualitative assessment of the data. The beam theory method uses only individual slices of MR data to determine participant-specific parameters, which, whilst less time consuming, provides much less detailed information. The results from the finite element model show that the peak stress can be located differently for any given time point of stance, with only 70% of participants showing a peak stress in the mid-shaft at mid stance. The 2D beam theory model cannot provide this information and only the peak dorsal and plantar stresses at the selected cross-section are obtained. Conversely the beam theory model allows computation of the stress during the entirety of stance, allowing 1D timeseries analyses such as Statistical Parametric Mapping to be used. This allows for detailed information about differences in stress magnitudes at different stance phases. However, in the context of understanding injury risk, it is believed that the magnitude of peak stress is important, and therefore stresses throughout stance that are submaximal may not be highly important.
Whilst it should be noted that the overestimation of the stresses by the beam theory model is a limitation to the interpretations of its results, there is currently no well-defined stress magnitude that will induce the formation of a stress fracture over a fixed number of cycles. Therefore, either model could be used to investigate differences between groupings of participants, but neither can discern if the stress estimated might lead to a stress fracture.
The limitations of the model used in this study have been outlined in Ellison et al. (2020a).This includes the limitation of reducing the foot to a single bone modelled as a fixed cantilever without a deformable arch and other plantar tissues, such as the long and short plantar ligaments. In terms of the data collection protocol, there are limitations introduced by asking participants to run barefoot, where they are habitually shod. Furthermore, NRF runners included those who run with a midfoot strike, forefoot strike and toe runners (Nunns et al., 2013), and these groups have some distinct characteristics.

Conclusion
Habitual rearfoot and non-rearfoot runners experience similar magnitudes of second metatarsal stresses during running, despite nonrearfoot runners experiencing greater external loading under the metatarsal head and greater joint contact forces. This further supports recent work suggesting that external loading measures should not be used as a proxy for internal bone loading.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.