Plantar fascia loading at different running speed: a dynamic finite element model prediction

ds on the plantar fascia could be influenced by running speed and relate to its pathology. This study calculated and compared plantar fascia strains under different running speed conditions using a dynamic finite element foot model and computational simulations. The model was previously validated featuring twenty bones, bulk soft tissue, muscles/ligaments, and a solid part of plantar fascia. A runner performed running trials under one preferred speed (PS), two lower (PS - 10% and PS - 20%) and two higher (PS + 10% and PS + 20%) speed conditions. The movement data were processed to drive musculoskeletal modelling and calculated boundary/loading conditions for the subsequent finite element analyses. The results show that peak strains of the plantar fascia increased with increasing running speed. From PS – 20% to PS + 20%, peak strain in the proximal and distal fascia regions increased by 96.78% and 58.89% respectively. Running speed could directly affect plantar fascia loading, which should be considered in running regimens and rehabilitation programmes. However, prescribing speed control for runners is worth pondering as it influences the trade-off between maximum single- step loads and loading frequency, which in coalescence determine the risk of plantar fascia injury and warranted further investigations. KEYWORDS Plantar fascia; running; tensile strain; finite element modelling; computational simulation CONTACT Ming Zhang Received 26 March 2020 ming.zhang@polyu.edu.hk et al., 2000; Hreljac, 2005). A positive relationship between ground reaction force, foot joint motion and running speed has been well established (Edwards et al., 2010; de David et al., 2015). Anecdotal information suggests that a similar relationship existed for loads on plantar fascia (Welk et al., 2015). A reduction in running speed is expected to alleviate tension on the plantar fascia and facilitate injury prevention, although evidence in the mechanical context is currently scarce. Identifying the loading patterns of the fascial band may therefore supplement the insight into the scenario. While in-vivo measurements of the plantar fascia force were usually difficult for human locomotion, computational simulation has emerged as a substitute method

is expected to alleviate tension on the plantar fascia and facilitate injury prevention, although evidence in the mechanical context is currently scarce. Identifying the loading patterns of the fascial band may therefore supplement the insight into the scenario. While in-vivo measurements of the plantar fascia force were usually difficult for human locomotion, computational simulation has emerged as a substitute method for solving the problem, which allows direct estimation of internal strains on foot soft tissues during a dynamic task. The purpose of this study is to determine the influence of changing running speed on the plantar fascia loading in running. Increased running speed was hypothesised to produce increments in tensile strains on the fascia band. To answer the research question, the stance phase of running trials was simulated under five speed conditions using an established finite element foot model. 2. Methods 2.1. Participant A healthy male (age: 30, height: 170.5 cm, weight: 1. Introduction Due to its known healthy benefits, running is a popular activity for many people (van Gent et al., 2007). Unfortunately, running-related injuries have remained highly prevalent as 26–75% of runners—both recreational and competitive ones – are estimated to sustain some forms of overuse injuries in a given year (Hryvniak et al., 2014). Plantar fasciitis was one of the most common foot problems occurring to runners (Taunton et al., 2002). More than two million individuals were diagnosed with plantar fasciitis in the United States annually (Roxas, 2005) and about 10% of the population were predicted to develop plantar fasciitis in their lifetime (Riddle et al., 2003). Plantar fascia injury was primarily caused by micro- tears on the fascial band secondary to repeated mechanical overload and excessive strain (Cutts et al., 2012). Plantar fascia attaches to the calcaneus and the five digits on the plantar facets and regulates motions of the foot arch through passive stretching (Wearing et al., 2006). External loading on the lower limb causes structural deformation of the foot arch and strains plantar fascia. Biomechanical factors in movement strategies, e.g. running speed, have been considered one of the factors that influences the force produced through

and the five digits on the plantar facets and regulates motions of the foot arch through passive stretching (Wearing et al., 2006). External loading on the lower limb causes structural deformation of the foot arch and strains plantar fascia. Biomechanical factors in movement strategies, e.g. running speed, have been considered one of the factors that influences the force produced through foot osteofascial structures and subsequent risk of injury (Wilk et al., 2000; Fredericson and Misra, 2007). Running speed has been frequently modulated for load management and injuries control for runners (Hreljac Plantar fascia loading at different running speed: a dynamic finite element model prediction Tony Lin-Wei Chen 1 , Duo Wai-Chi Wong 1,2 , Yan Wang 1,2 , Ming Zhang 1,2 1 Department of Biomedical Engineering, Faculty of Engineering, The Hong Kong Polytechnic University, Hong Kong, People’s Republic of China 2 The Hong Kong Polytechnic University Shenzhen Research Institute, Shenzhen, People’s Republic of China

HKIE Transactions | Volume 28, Number 1, pp.14-21 15 68.1 kg) volunteered for the study. He was an experienced runner with a weekly training volume of 15 kilometres and a usual pacing of 6 minutes per kilometre. He was recruited for both computer model construction and gait analysis. The participant did not report any musculoskeletal diseases/ injuries that could influence his running performance at the time of experiment. He was fully informed of the research procedure and signed the consent form. The experimental protocol was approved by the institution (NO. HSEARS20170626003). 2.2. Experimental procedure The experiment comprised several running trials under five running speed conditions: the runner’s preferred speed (PS, 10 km/h as reported by the runner), two decreased speed conditions (PS - 10% and PS - 20%), and two increased speed conditions (PS + 10% and PS + 20%). The participant was given ample time to acclimatise to the running condition before the experiments. The five speed conditions were tested in a random sequence with a five-minute rest interval. Footwear was excluded from the tests so as to eliminate its effects on the running gait pattern. Data collection began with a static calibration trial that would later be used to scale an established musculoskeletal model. Afterwards, the participant ran through the motion capture region where two pairs of photoelectric cells were placed 2.6 m apart along the runway to monitor running speed (Hamill et al., 2014). Data generated from the running trials were considered valid when the measured running speed fell with 5% variance of the target value and the runner’s footsteps landed completely within the force platforms. A motion capture system with eight cameras (Vicon, Oxford Metrics Ltd., Oxford, UK) and two force platforms (OR6, AMTI, Watertown, USA) were used to record marker trajectories and ground reaction force at the sampling rate of 250 Hz and 1000 Hz respectively. Thirty- six retroreflective markers were affixed to the following anatomic landmarks: acromioclavicular joints, posterior/ anterior iliac spines, greater trochanters, lateral/medial femoral epicondyles, lateral/medial malleoli, calcaneal tuberosity, the base/head of the first and fifth metatarsals, and the distal phalanx of

were used to record marker trajectories and ground reaction force at the sampling rate of 250 Hz and 1000 Hz respectively. Thirty- six retroreflective markers were affixed to the following anatomic landmarks: acromioclavicular joints, posterior/ anterior iliac spines, greater trochanters, lateral/medial femoral epicondyles, lateral/medial malleoli, calcaneal tuberosity, the base/head of the first and fifth metatarsals, and the distal phalanx of the hallux. 2.3. Musculoskeletal model A representative trial of each speed condition was selected for computational simulations. The kinematic and kinetic data were first input to the OpenSim platform (version 4.0, National Centre for Simulation in Rehabilitation Research, Stanford, USA). A generic model (John et al., 2013) with 12 rigid-body segments, 23 degrees of freedom, and 92 musculotendinous units was scaled to accommodate the mass and anthropometry of the participant. Inverse kinematics was then solved and the dynamic inconsistency was reduced by small adjustments to model mass properties (Arnold et al., 2010). The built- in tool Computed Muscle Control was implemented to compute muscle excitations that drove the desired movement patterns. Upon completing the analysis steps, forces of the external foot muscles, joint reaction force, and segmental kinematics were exported as the boundary and loading conditions for the subsequent finite element analysis. 2.4. Finite element model 2.4.1 Geometry reconstruction A finite element model was built on the MRIs of the participant’s left foot. The MRIs were acquired by a 3.0 scanner (GoldSeal Certified Signa HDxt, General Electric Company, Boston, USA) at T1 sequence, 1-mm slice interval, and 0.625 mm pixel size with the participant’s leg fixed at neutral position by a customised ankle-foot orthosis. The image was processed to obtain the 3D foot geometries with an aid of segmentation software (Mimics and 3-matics, Materialise, Leuven, Belgium). The foot model was validated against a cadaveric experiment in our previous study (Chen et al., 2020). Simply speaking, the foot model was proximally fixed and sustained a pulling force of 98.1 N on the plantar forefoot to mimic a loading condition of plantar fascia-stretching exercise. Strains on selected sites of the fascia band were compared between our simulations and the cadaveric study. A

validated against a cadaveric experiment in our previous study (Chen et al., 2020). Simply speaking, the foot model was proximally fixed and sustained a pulling force of 98.1 N on the plantar forefoot to mimic a loading condition of plantar fascia-stretching exercise. Strains on selected sites of the fascia band were compared between our simulations and the cadaveric study. A good agreement was established between the predicted and measured strain values through a regression analysis (Chen et al., 2020). The model included 20 bony segments, 26 ligaments, 9 intrinsic foot muscles, and 11 extrinsic foot muscles, as shown in Figure 1 (Chen et al., 2020). The bulk soft tissue was modelled by a cluster of smooth particle hydrodynamics (SPH) particles encapsulated in a shell unit using the cohesive property. The SPH particles were shown to validly incorporate themselves into regular deformable elements in the finite element domain (Chen, Wong, et al., 2019) and had the merits of solving high-impact problems (Zhang and Batra, 2009). The shell unit consisted of interior and exterior surfaces representing the periosteum layer (tied to the bony structures) and skin layer (in contact to the ground) respectively. Ligaments and intrinsic foot muscles were modelled as truss units while the extrinsic foot muscles were embodied by mechanical connectors. Both bones and plantar fascia were reconstructed as 3D solids. The geometry and placement of the plantar fascia were determined by its manifests in MRIs and the human anatomy atlas (Standring, 2020) using computer-aided design software (Solidworks v2014, Dassault Systèmes, MA, USA). From the origin to distal end, the thickness of the fascial bands gradually reduced from 3.0 mm to 1.3 mm. The plantar fascia was tied to the inferior calcaneus proximally and to the heads of five phalanges distally.

T L W CHEN ET AL. 16HKIE Transactions | Volume 28, Number 1, pp.14-21 Figure 1. Outline of the finite element modelling and the boundary/loading conditions. Solid black arrows denote the names of the parts. Solid maize arrows represent boundary/loading conditions. Dashed frames denote the interactions among parts. The bulk soft tissue was Figure 1. Outline of the finite element modelling and the boundary/loading conditions. Solid black arrows denote the names of the parts. Solid maize arrows represent boundary/loading conditions. Dashed frames denote the interactions among parts. The bulk soft tissue was modelled as SPH particle elements and encapsulated in a shell unit that possessed an interior profundal fascia layer (tied to the skeletal structures) and exterior skin layer (in contact with the fixed ground). The foot segments were connected by ligaments (truss unit), intrinsic foot muscles (truss unit), and the plantar fascia (three-dimensional solid). Three-dimensional ankle joint reaction force, extrinsic foot muscle force, and initial transitional velocity were applied to the model to drive the simulation. DOFs = degrees of freedom; SPH = smoothed-particle hydrodynamics.

HKIE Transactions | Volume 28, Number 1, pp.14-21 17 2.4.2 Material properties and mesh creation All model components were assumed linearly elastic and isotropic except that plantar fascia was assumed to be hyperelastic and isotropic. Parameters for the material equations were acquired from the literature (Table 1). The model mesh was created by Abaqus v6.14 (Simulia, Dassault Systèmes, RI, USA). Previous mesh convergence tests had determined an optimal global mesh size of 3.5 mm for the osseous and soft tissue components (Chen, Wong, et al., 2019). A frictionless contact algorithm was used for bone-to-bone interface to mimic the function of cartilages (Athanasiou et al., 1998). The interaction between the skin and the ground plate was “hard” contact with a friction coefficient of 0.6 (Zhang and Mak, 1999). 2.4.3 Boundary and loading conditions Boundary and loading conditions for the finite element analyses were input as a tabulated time-series matrix from preceding musculoskeletal modelling. Concentric connector force was applied to the slip-ring connectors to simulate extrinsic foot muscle force. Three-dimensional ankle joint reaction force was loaded on the ankle surface of the talus (Figure 1). The whole foot model was initially positioned and oriented corresponding to the instant before foot strike. A pre-defined velocity was also assigned to the whole foot model. Gravity was enabled throughout the simulation steps using the force/mass ratio of 9.8. 2.5 Simulation solver and data output The Abaqus dynamic explicit solver (version 6.14, Simulia, Dassault Systèmes, RI, USA) was used to perform running simulations for each speed condition. To avoid artefacts of extreme deformation, principal tensile strain was reported for the central fascial band excluding regions that were attached to the bones at the proximal and distal ends. The central plantar fascia was equally divided into three portions from proximal to distal. Principal tensile strain was averaged across all elements of each portion and the mean values were compared across the conditions. Strain contours of the plantar fascia were also plotted. 3. Results Figure 2 displays the principal tensile strains on the plantar fascia (from the proximal to the distal regions) as a function of percentile stance

portions from proximal to distal. Principal tensile strain was averaged across all elements of each portion and the mean values were compared across the conditions. Strain contours of the plantar fascia were also plotted. 3. Results Figure 2 displays the principal tensile strains on the plantar fascia (from the proximal to the distal regions) as a function of percentile stance phase. The fascial band was increasingly loaded after initial contact and maximally loaded during the mid-to-terminal stance phase (50–90% stance). The results show that running speed posed an effect on fascia strain level. Compared with the PS condition, peak strain in the proximal fascia region decreased by 7.89% and 4.32% in the PS - 20% and PS - 10% conditions, and increased by 13.74% and 17.08% in the PS + 10% and PS + 20% conditions. Likewise, peak strain in the distal fascia region decreased by 3.95% and 5.52% in the PS - 20% and PS - 10% conditions, and increased by 15.65% and 23.17% in the PS + 10% and PS + 20% conditions. Peak strains in the mid fascia region were most reduced in the PS - 20% condition (by 20.45%) but roughly comparable between the PS and two acceleration conditions. Figure 3 shows the strain maps of the plantar fascia when it was most loaded at late stance for the five speed conditions. As running speed increased, there was clearly a wider distribution of warn- colour coded regions in the fascia band, which indicated a higher strain level according to the colour atlas. Figure 2. Changes in tensile strains on the three portions of the plantar fascia as a function of percentile running stance. A: the proximal portion; B: the middle portion; C: the distal portion.

T L W CHEN ET AL. 18HKIE Transactions | Volume 28, Number 1, pp.14-21 Figure 2. Changes in tensile strains on the three portions of the plantar fascia as a function of percentile running stance. A: the proximal portion; B: the middle portion; C: the distal portion. Figure 3. Strain maps of the plantar fascia for the five running speed conditions. The contour plots were extracted at the instant of maximal strain for each condition. The images are colour-coded based on the distribution of principal tensile strains. The highest strain is marked in red, and the lowest strain is marked in blue. Regions of the plantar fascia tied to the bony segments are removed from display. 4. Discussion Running towards higher speeds increases forces that apply to the lower limb (Schache et al., 2011; de David et al., 2015). Increased mechanical burden to the foot was believed to tighten the plantar foot soft tissues and exacerbate foot problems (Chen, Agresta, et al., 2019). However, there is limited quantitative evidence reflecting tension in the plantar fascia in Figure 3. Strain maps of the plantar fascia for the five running speed conditions. The contour plots were extracted at the instant of maximal strain for each condition. The images are colour-coded based on the distribution of principal tensile strains. The highest strain is marked in red, and the lowest strain is marked in blue. Regions of the plantar fascia tied to the bony segments are removed from display. 4. Discussion Running towards higher speeds increases forces that apply to the lower limb (Schache et al., 2011; de David et al., 2015). Increased mechanical burden to the foot was believed to tighten the plantar foot soft tissues and exacerbate foot problems (Chen, Agresta, et al., 2019). However, there is limited quantitative evidence reflecting tension in the plantar fascia in running tasks, which could be associated with injury risk, fitness, and rehabilitation program. For this reason, we calculated the plantar fascia strains during the stance of running gait and examined how it could be influenced by different running speeds. In accordance with our hypothesis, tensile

al., 2019). However, there is limited quantitative evidence reflecting tension in the plantar fascia in running tasks, which could be associated with injury risk, fitness, and rehabilitation program. For this reason, we calculated the plantar fascia strains during the stance of running gait and examined how it could be influenced by different running speeds. In accordance with our hypothesis, tensile strains on the fascia band (from proximal to distal portions) increased as the runner accelerated during the running trials. The increments in fascial strains were most distinct in the transition from mid- to terminal-stance. Similar to our previous work, the magnitude of peak fascial strain (1.22–3.22%) calculated in the current study agreed with a normal strain range (1.51–2.73%) for common activities, such as walking (Erdemir et al., 2004; McDonald et al., 2016), but apparently fell short to those (14.52%) reported by Chen (Chen et al., 2014). In Chen’s study, the effects of foot muscles and other planter soft tissues were omitted, which resulted in a large portion (20.1%) of the total foot loading solely sustained by the plantar fascia. This percentage of load sharing was considered to be overestimated because the yield strain and force for the dissected plantar fascia specimen were around 15% (Butler et al., 1986) and 1189 N (Kitaoka et al., 1994) respectively. Reaching over 90% of strength limit of the fascia material during a single walking step was physically improbable. The discrepancies in reported strain level of the plantar fascia among the studies reflect the existence of the synergic effects of other plantar foot soft tissues for foot supports (Kirby, 2017). Our predicted timing of peak fascial strain (80–90% stance phase) was in line with those (80–83.3% stance phase) of both Chen’s (Chen et al., 2014) and Erdemir’s (Erdemir et al., 2004). During the period, the plantar fascia was considered most tensioned by the activated windlass and turned the foot arch into a rigid lever for propelling the body forward (Bolgla and Malone, 2004). In line with previous studies (Hamill et al., 1983; Edwards et al., 2010; Schache et al., 2011; de David et al., 2015;

and Erdemir’s (Erdemir et al., 2004). During the period, the plantar fascia was considered most tensioned by the activated windlass and turned the foot arch into a rigid lever for propelling the body forward (Bolgla and Malone, 2004). In line with previous studies (Hamill et al., 1983; Edwards et al., 2010; Schache et al., 2011; de David et al., 2015; Boey et al., 2017), the increased fascia strain from slow to fast running conditions was probably attributed to change in force transmission along the lower leg. Without adjusting the landing strategy such as gait retraining, runners usually experienced increased impacts on the lower leg to cope with shock absorption in fast running. Increased running speed was associated with larger ground reaction force (Hamill et al., 1983; Nilsson and Thorstensson, 1989) and tibial deceleration (Boey et al., 2017). During the stance phase, the two forces apply to the base and the apex of the foot arch and form a three-point bending motion to the foot (Lieberman, 2012; Perl et al., 2012). Increments in both forces further deform the foot arch construct and bring in tension to the plantar connective tissues (Chen, Agresta, et al., 2019). In this study, the fascia strain increased by around 4% when running speed increased by 20%. These changes, though seemingly small in magnitude, were reportedly equivalent to windlassing the fascia band by dorsiflexing the first metatarsophalangeal joint by around 15° (Carlson et al., 2000). A slight reduction in fascia strains appeared to improve foot joint motions and benefit runners with constrained foot mobility because of a tightened fascia band (Irving et al., 2006). To avoid possible injury risk, research suggest decreasing gait speed as a means of controlling loading to the foot soft tissues and benefiting symptomatic runners (Robbins and Maly, 2009). However, pathomechanics of plantar fasciitis resembles the procedure for material fatigue (Chen, Wong, et al., 2019). It means that both loading magnitude and accumulated loading cycles determine the onset of injuries. Concerns were raised: reducing running speed but maintaining running volume could expose runners to prolonged loads and increased loading frequency on vulnerable anatomic

This study calculated and compared plantar fascia strains under different running speed conditions.