Viscoelasticity in Foot-Ground Interaction

Mechanical properties of the plantar soft tissue, which acts as the interface between the skeleton and the ground, play an important role in distributing the force underneath the foot and in influencing the load transfer to the entire body during weight-bearing activities. Hence, understanding the mechanical behaviour of the plantar soft tissue and the mathematical equations that govern such behaviour can have important applica‐ tions in investigating the effect of disease and injuries on soft tissue function. The plantar soft tissue of the foot shows a viscoelastic behaviour, where the reaction force is not only dependent on the amount of deformation but also influenced by the deformation rate. This chapter provides an insight into the mechanical behaviour of plantar soft tissue during loading with specific emphasis on heel pad, which is the first point of contact during normal gait. Furthermore, the methods of assessing the mechanical behaviour including the in vitro/in situ and in vivo are discussed, and examples of creep, stress relaxation, rate dependency and hysteresis behaviour of the heel pad are shown. In addition, the viscoelastic models that represent the mechanical behaviour of the plantar soft tissue under load along with the equations that govern this behaviour are elaborated and discussed.


Introduction
The human heel pad is usually the first part of foot that contacts the ground during normal gait. The soft tissue structure, which is located underneath the calcaneus (heel bone) consists of fat pad and skin. While this fat pad, also known as corpus adiposum, works as shock absorber and dampens the ground reaction forces during weight-bearing activities like standing and walking, the skin has another important role to prevent tear and to work as an impermeable barrier to protect the underlying soft tissue [1]. Reported heel pad thickness varies from 12.5 to 24.5mm in different studies [2][3][4][5][6][7][8] using different imaging modalities including ultrasound [2][3][4], magnetic resonance imaging (MRI) [5] and radiography [6,7].
The plantar soft tissue structure is designed to bear large loads. Similar type of adipose tissue is found in other parts of body that normally need to bear compressive and shear stress, such as fingertips. The structure of the plantar adipocytes consists of a dense network of septa, which prevents free movement of fat cells while allows the lateral movement [9]. The unique structure of the plantar soft tissue enables it to bear large strain in reaction to the ground force. In each heel strike, vertical loading roughly equal to 110% body weight is applied on the heel, whereas 25% of the body weight is applied in anteroposterior and 10% in mediolateral directions [10]. Normalising the loading over an area results in an average pressure around 100-400 kPa for healthy individuals depending on the site of the foot [11]. These plantar pressure values can increase as a result of increase in the stiffness of plantar soft tissue, i.e., due to diabetes that results in a decrease in the contact area [11].
Furthermore, understanding the effect of different pathologies such as diabetes on the mechanical properties of human plantar soft tissue is paramount. While these pathological conditions may not affect the structure (geometry) of the plantar soft tissue they would affect the mechanical properties of the plantar soft tissue [12]. The knowledge of these mechanical properties that determine the behaviour of tissue under load can be utilised for diagnosis of foot pathologies as well as for treatment interventions such as foot orthoses and footwear.
In order to understand the mechanical behaviour of the plantar soft tissue, it is necessary to have an overview of the basic mechanical definitions.

Principles of mechanics of materials
Stress is a quantity which expresses the force that neighbouring particles apply on each other. Stress can also be defined as the amount of force, which is required for a unit of cross-section area to compress or extend in the normal direction. F A s = (1) where σ represents stress, F is the load and A is the cross-section area. The unit of the stress is Nm −2 .
Strain is the change in the length of the object in the axial direction which is normal to the surface of applied load. Strain can be defined as the amount of change in the length of the object over the original length as a result of applying load.
where ε is strain, dL is the change in length and L is the original length ( Figure 1); hence, strain does not have a unit. Shear stress strain is the response of a material to the force, which is applied parallel to the specific surface. This force makes the geometry of the structure to deform but not to stretch/ compress.

Parallel surface to force
where τ is the shear stress, F is the applied force parallel to the surface and A is the cross-section area (Figure 2).
The gradient of the force-deformation graph describes the stiffness of the material, and it is the quantification of the rigidity of the material and is expressed in Nm −1 .
dF Stiffness dL = (4) where dF and dL are changes in load and displacement, respectively. Figure 2. Normal deformation as a result of tensile and compressive forces applied to cylindrical specimens along with the shear deformation as a result of shear force applied to a cubic specimen.

Principles of elastic solid material
Elasticity is the ability of a material to resist force and return to its original shape when the force is removed. Elastic solid materials are divided into two main groups: Hookean and non-Hookean [13].
In the Hookean material where stress increases linearly with increase in strain, the slope of the stress-strain graph is defined as Young's modulus or modulus of elasticity (E) expressed in Pa.
where σ is stress and ε is the strain.
However, in the non-Hookean material stress is not linearly proportional to strain. The relationship between stress and strain changes during different stages of loading.

Principles of viscous fluid materials
A fluid material is defined as a material that bears shear deformation and consists of liquids and solids [14]. Liquids consist of atoms with interatomic connections and molecules with weak intermolecular connections. The shear force can be applied to break the weak intermolecular bonds to allow the material to flow [15].
In continuum mechanics, a Newtonian fluid is a fluid in which the arising viscous stresses are linearly proportional to the local strain rate [16]. In other words, the shear stress is proportional to the rate of change of the fluid's velocity vector.
. dv dy t m = (6) where τ represents the shear stress, μ represents viscosity and dv/dy represents the velocity gradient in the direction perpendicular to the velocity.
The simplest mathematical models that take viscosity into account can be applied for Newtonian fluids. Although there is no real fluid that fits this definition properly, many common liquids and gases, such as water and air, are assumed to be Newtonian under ordinary conditions.
On the other hand, a non-Newtonian fluid is the one with different properties from a Newtonian fluid. To be more specific, the viscosity, which is the quantity of a fluid's ability to resist gradual deformation by shear or tensile stresses, depends on shear rate or shear rate history [16].

Viscoelastic behaviour
A viscoelastic material combines properties of the elastic solid with viscous fluid. The elastic solid can be Hookean or non-Hookean and the viscous fluid can be Newtonian or non-Newtonian [17]. There is a variety of behaviour within different viscoelastic materials which ranges from completely elastic solid behaviour to completely viscous fluid behaviour [17]. In addition, the viscoelastic material has the distinctive characteristics which show both viscous and elastic behaviour when exposed to loading. The soft tissue exhibits a viscoelastic behaviour under compression which means that the force-deformation behaviour depends on the amount of deformation and deformation rate [18]. Viscoelastic materials behave in different ways under various types of loading exhibiting differences in deformation/force rate dependency, creep, stress relaxation and hysteresis [18]. For example, a viscoelastic material under cyclic loading behaves differently during loading and unloading. In a sense, the stiffness of the material decreases during unloading when compared to the stiffness that material shows during loading. The area between force-deformation graph in loading and unloading is called hysteresis that represents dissipated energy [18] (Figure 3). The load deformation graph of a heel pad tested using ultrasound indentation technique [19]. The different colours represent data gathered at different deformation rates as presented in the legend. The area surrounded between the loading and unloading curves represents hysteresis.
The force-deformation behaviour of the plantar soft tissue like other biological materials is influenced by the loading velocity [18] (Figure 3). Stress relaxation characteristics indicate that when a viscoelastic material is deformed suddenly and the deformation is kept constant for a specific time, the force decreases with time [18] (Figure 4). The creep characteristics indicate that when a specific amount of load is suddenly applied to the viscoelastic material and kept constant for a specific time, there would be an increase in deformation over time [18] (Figure 5).

Mechanics of the heel pad
Because of the liquid content of the heel pad tissue along with the arrangement of solid components which has an influence on regulating the fluid movement, the heel pad behaviour is mostly assumed as that of a nearly incompressible material [20].

The structure
The fat pad consists of two layers: microchambers and macrochambers. Microchambers shape the plantar layer of heel pad, which protects the fat pad from excessive bulging during loading [12]. The deeper layer is composed of sparser, fibro-adipose structure called a macrochamber [21]. The microchambers is the thinner layer of the two and contains mainly elastic fibres, but the thicker layer (macrochambers) contains roughly equal amount of elastic fibres and collagen [22]. Therefore, a different behaviour between two layers is expected [12], which is discussed under the mechanical behaviour section (Figure 6). The heel pad atrophy is a clinical condition, which is usually linked with diabetes, collagen disorders and peripheral neuropathy [22]. The results of histological studies have revealed that the adipocytes in the subcutaneous layers closer to the foot surface of the fat pad were 25% smaller in the mean cell area and 10% smaller in mean maximum diameter in an atrophic heel compared to a normal heel [22]. Additionally, the adipocytes in the deep subcutaneous layers are 45% smaller in the mean cell area and 25% smaller in mean maximum diameter in an atrophic heel compared to a normal heel [22]. Septa in both deep and superficial subcutaneous layers in the atrophic heel are 25% thicker than in a normal heel and include more percentage of elastic tissue which seems to be uneven in some cases [23].
Furthermore, in atrophic feet, collagen septa are found to be thicker, and also the adipose cells are smaller than healthy feet [22,24]. The amount of internal stress and strain depends on the structure of the heel (geometry) as well as on the material properties of the tissue and also the value and direction of force. Hence, the above-mentioned changes in the structure of the heel pad can increase the internal stress and strain that are claimed to be the main causes of tissue injuries [25][26][27]. In addition, the interface between the soft tissue structure and underlying bony prominence can be the area of high stress concentration and it is claimed that tissue damage starts in deep tissue close to the bony prominences and then develops up to the skin surface [25,26,28,29].

The function
The foot supports stabilising the whole body during standing and is the interface between the body and the ground during walking. The heel pad as the first part of the foot which has contact with the ground during locomotion acts as a shock absorber and shock reducer to protect the foot from local stress [2]. The strain and pressure applied to the heel during gait can be withstood by the honeycomb structure of the heel pad. The heel pad as a structure shows viscoelastic behaviour and provides cushioning during heel strike and absorbs shocks by dissipating energy [2]. This mechanical energy dissipates as heat just after heel strike and decreases the possibility of mechanical trauma to the foot [2].
The heel pad can absorb the impact shock during heel strike by carrying out deformation under loading and by distributing the force over a wider area of the skin to prevent stress concentration [30]. The chambered structure helps to spread the compressive force over the whole area of plantar surface of the bone and prevents any injuries on the calcaneus bone during heel strike [2]. Under compression, the heel pad expands easily as a result of low stiffness, but later the tension on the collagen fibres of the fat pad and skin limits the movement of the tissue and increases the stiffness of the heel pad gradually in the loading direction [31]. This justifies the nonlinear force-deformation behaviour of the heel pad and causes the strain stiffening inherent in the heel pad's mechanical behaviour.
Additionally, it was reported that the mechanical behaviour of microchambers and macrochambers is different under compression [12]. The microchamber layer experiences less strain compared to the macrochamber layer [12]. Hsu et al. [12] indicated that the stiffness of the microchambers is 10 times greater than macrochambers, and concluded that the observed difference plays an important role in the heel pad mechanical behaviour. It appears that the macrochamber layer is responsible for large deformation and cushioning behaviour of the heel pad during gait; however, the microchamber layer is responsible in preventing the heel pad from excessive bulging [12].
The main role of the heel pad is to decrease the impact shock during heel strike and to distribute pressure during the foot-ground contact through undergoing deformation. The deformability of the heel pad may reduce through tearing the fibrous septa or atrophy of heel pad due to a trauma or ageing [32,33]. After severe injuries such as tearing or breaking of the honeycomb structure, the heel pad does not have the ability to remodel itself [32,33]. As the fat pad is a semi-liquid structure with hydrostatic properties of fluids [33], a decrease in water content of the heel as well as a decrease in the elastic fibrous tissue and loss of collagen are the main reasons for the gradual weakening of the tissue due to ageing [32]. Overall, it can be concluded that the loss of soft tissue substances due to ageing, atrophy or any previous injury prevents the tissue from responding to load in an optimal way [33].

Mechanical behaviour of heel pad
Plantar soft tissue like other biological tissues has a nonlinear elastic behaviour. Initially by applying small amount of load, the tissue deforms easily (low stiffness) and by increasing the loads the stiffness increases gradually [34]. The heel pad expands easily as a result of low stiffness but afterwards the tension on the collagen fibres of the fat pad and skin limits the movement of the tissue and increases the stiffness of the heel pad in the loading direction [31].
Pathological changes in foot may not be detectable from the structure of the foot but normally correlate with the alteration in the mechanical behaviour of the tissue [12]. Therefore, quanti-fying the mechanical properties of the plantar soft tissue is important to assess the risk of mechanical trauma to the foot.

Method of assessing mechanical behaviour of heel pad
Several methods have been used to extract the mechanical behaviour of the heel pad that can be divided into three main groups: in vitro, in situ and in vivo.

In vitro tests
In some studies, the heel pad behaviour was characterised using an in vitro or in situ method to quantify the material properties of the plantar soft tissue [32,[35][36][37][38][39][40]. Miller-Young et al. [35] performed a series of tests on the cylindrical samples of plantar soft tissue extracted from cadaveric feet. Three series of tests that were performed on the samples included the quasistatic test to obtain the hyperelastic mechanical properties of the soft tissue; the stressrelaxation tests to calculate the viscoelastic time constants; and the dynamic compression test to extract the viscoelastic relaxation material coefficients. The reported results clearly showed the time dependency and viscoelastic behaviour of the heel pad [35].
Some other studies [39,41] used 2 × 2 cm samples from different sites of plantar soft tissue of cadaveric feet to calculate the elastic and viscoelastic coefficients of a mathematical model for plantar fat pad to compare the properties of the plantar soft tissue between a healthy and a diabetic foot. The results showed that stiffness and energy dissipation increase with loading frequency [39,41]. This was completely in contrast with Bennet and Ker's results, which indicated no changes in energy dissipation and stiffness with different testing frequencies [36]. The results also showed that frequency dependency is higher in younger subjects and this was attributed to the difference in the heel pad's water content in younger versus older tissue specimen [32,37]. The water content of the soft tissue can be considered as the main reason for viscoelastic behaviour of the soft tissue. Therefore, a decrease in water content of the soft tissue in older subjects can lead to a decrease in the viscoelastic characteristic of the soft tissue such as loading frequency dependency.
While there was no significant difference in the values of viscoelastic coefficients between diabetic and non-diabetic feet, it was claimed that changes in plantar soft tissue were mainly recognised at the structural level, and not reflected effectively at the material level [38,39].

In situ tests
Aerts et al. [40] compared the results of in vitro tests with the in situ tests in which the heel bone was fixed to the wall and a pendulum was used to impact the heel region.
In this study, differences between the test results in two conditions were observed and they were attributed to the differences in soft tissue behaviour at structural and material level [40]. In another study, Bennet and Ker [36] compared the results of in vitro versus in situ tests. The researchers performed two series of tests in which one group of specimens was tested while attached to the calcaneus and the surrounding tissues, while the other group of specimens were tested after complete removal from the calcaneus [36]. In this study, the specimens were tested using a dynamic loading machine and load-displacement data were recorded during the test. The energy dissipation ratios and stiffness were higher in isolated heel pads compared to the heel pads attached to the calcaneus [36]. Therefore, it was concluded that the results of the in vitro test show the properties of the isolated heel pad material, while the in situ test results represent the behaviour of the heel pad structure [40]. The fact that the mechanical properties of heel pad extracted from the in vitro and in situ results are different can be explained by the indications that structural factors such as heel pad skin and geometry of the calcaneus have an influence on the heel pad behaviour.
Although the in vitro test can provide more repeatable data due to eliminating geometrical complexities of the plantar soft tissue, it cannot provide a realistic assessment of the mechanical behaviour of the plantar soft tissue during weight-bearing activities [38,39]. Furthermore, the in situ tests cannot be used to assess the mechanical behaviour of the heel pad in different individuals.

In vivo assessments
In a number of studies, the human heel pad was characterised using the force-deformation data extracted from in vivo experiments.
Investigating the mechanical behaviour of the heel pad during walking can enhance our knowledge about the heel pad in realistic loading conditions.
In one study, radiographic fluoroscopy was utilised to measure the thickness of the heel pad during walking [42] while the plantar pressure was also measured using optical display method. Their results showed that maximum strain of the heel pad during walking was 40% and the absorbed energy ratio was estimated as 17.8% (SD 0.8) for different velocities (0.5-0.9 m/s) [43], which is considerably lower than 35% ratio that was reported for in vitro studies of the heel pad [40,41].
While the radiographic method allows measuring the deformation of the heel pad during actual gait, the control on the direction of loading affects the results and this may lead to variation in the observed results.
The device commonly consists of linear array ultrasound probe, load cell, motor and mechanical body, which is composed of a foot place perpendicular to the axial of loading (same axis to the probe head). The ultrasound probe is necessary to measure the deformation of the soft tissue. The force that is applied to the foot to compress/indent the foot can be measured by a load cell. It can be mounted at the back of the probe to measure the axial force. The mechanical part should be mounted in such a way that the probe can be compressed/indented to different feet sizes. The motor can generate uniform movement of the probe. A number of studies which conducted experiments with ultrasound indentation device [52,53] used a custom loading device that consists of a linear array ultrasound probe which was in series with a dynamometer (load cell) and mounted on a rigid frame (Figure 7). Figure 7. The ultrasound indentation device and a schematic representation of the procedure followed to create the tissue's force/deformation curve [52].
In addition to the use of indentation device in experimental analyses of the plantar soft tissue, this device has been commonly used to quantify the mathematical and finite element (FE) models that govern the behaviour of soft tissue during loading. This is discussed in the next section under mechanical behaviour model. Ultrasound strain elastography has been used for assessment of plantar soft tissue stiffness in patients with diabetic neuropathy and was recognised to have potentials for diagnosing tissue mechanical malfunction in clinical setting [54].

Plantar soft tissue stiffness and measurement methods
The mechanical properties of the plantar soft tissue show a high level of dependency on the measurement method [40,55]. For example, it was reported that the stiffness of the heel pad of humans using in vitro method is almost six times higher compared to the stiffness measured during in vivo tests, while the absorbed energy ratio is about three times lower using the in vitro method [40]. This can be explained by the indications that structural factors such as heel pad thickness and geometry of the calcaneus have a significant influence on the heel pad mechanical behaviour [3,23,56]. Aerts et al. [40] compared the energy dissipation and stiffness of the soft tissue in an amputated leg and in an isolated heel pad. They showed that the whole lower leg which is involved during in vivo test affects the test results in different ways and influences stiffness and energy dissipation in in situ test. The presence of lower leg makes a difference in terms of limiting expansion in some directions and dissipating the energy [40].
Furthermore, while indentation seems to provide more realistic and reliable assessment of the mechanical behaviour of the plantar soft tissue during in vivo conditions, the effects of indenter's shape (i.e. probe's head geometry) together with the effects of calcaneus bone geometry need to be taken into account in analysing the results from these type of studies.

Changes to mechanical behaviour
The mechanical behaviour of the plantar soft tissue can be changed by ageing, heel pain and other pathologies; some of which are described below.

Effect of ageing on the mechanical behaviour
Hsu et al. [57] compared the average unloaded heel pad thickness in a young group and an elderly group of participants. The average unloaded heel pad thickness was 20.1 (±2.4) mm in the elderly group and 17.6 (±2.0) mm in the young group. These results are in line with the findings by Kwan et al. [58] who showed that the thickness of the soft tissue was higher in the elderly group and that the stiffness increases with age [58]. The stiffening of the soft tissue by ageing may reduce the adaptability of the tissue in responding to the stress, which may lead to foot diseases in elderly people [59].
The shock absorption of the heel pad was determined in two different age groups for two different impact velocities by Kinoshita et al. [60]. It was reported that the absorbed energy in younger adults is significantly higher than in elderly [60]. The energy absorbed density depends on the viscoelastic properties of the plantar soft tissue and can be calculated by subtracting the energy return density from the energy input density [53]. The energy dissipation of soft tissue is directly linked to its viscoelastic characteristic. Therefore, it can be concluded that the viscoelastic properties of the soft tissue alter with age.

Effect of heel pain on mechanical behaviour
Physical activity and repetitive high-impact loading during sports activity can cause some micro-damage in the heel pad and consequently lead to heel pain. The sports that involve running or jumping apply repetitive impact loading on the heel pad and may cause collapse of the heel pad and finally cause a greater force impact on the calcaneus.
The inflammatory oedema may be a sign of changes in the structure of the heel pad that can lead to decrease in shock absorption capability during heel strike [61]. To compare the mechanical behaviour of the heel pad, a parameter known as heel pad compressibility index (HPCI) was defined as the ratio of the loaded tissue thickness to unloaded tissue thickness in percentage [46]. This parameter was utilised by Prichasuk et al. [62] to compare between normal and painful heel and revealed that the compressibility index increases in the patients with heel pain.
Tong et al. [50] compared the plantar tissue thickness and HPCI in normal, plantar heel pain and participants with diabetes. The results showed that compressibility in patients with diabetes and heel pain was less than the corresponding value in healthy volunteers. While the findings of Tong et al. [50] on the effect of ageing on HPCI contradict Prichasuk et al.'s findings [62], in another study Ozdemir et al. [63] reported the increase in heel pad thickness and HPCI with ageing and increase in body weight. This was attributed to the gradual loss of collagen, water content and elastic fibrous tissue of the heel pad as result of ageing [60], all of which can lead to a change in the viscoelastic behaviour of the heel pad.

Effect of other pathologies on mechanical behaviour
The structural changes in the soft tissue of diabetic patients will cause some changes in the macroscopic and microscopic behaviour of plantar soft tissue and make it more vulnerable to mechanical trauma which can lead to ulceration [54]. Less elastic tissue and impaired ability in distributing pressure are the other changes in diabetic tissues which can lead to a weakened cushioning effect [54]. Increase in energy dissipation during weight bearing is the other factor which can increase the risk of ulceration [54].
A tissue with increased viscosity and decreased elasticity may provide the similar amount of stiffness during loading, but during the swing phase of gait and as the tissue unloads the increased viscosity may not allow the tissue to completely return to its original thickness. This may cause the next loading cycle to start with a tissue that is partially deformed, which may increase the likelihood of tissue bottoming out.
In addition, generally the stiffening and thinning of the fat pad would make the tissue more fragile that makes it more likely to damage compared to a healthy tissue [64]. This also applies to the behaviour of skin in a diabetic foot, which can be less flexible and more brittle when it becomes drier [50]. It was also found that the deformability of the heel pad is less in participants with diabetes compared to their healthy counterparts [50]. However Hsu et al. [46] measured the elastic modulus, compressibility index and unloaded tissue thickness in people with diabetes and found no difference with the respective tissue characteristics in healthy participants. Hsu et al. [48] compared the strain and elastic modulus in macrochambers and microchambers of the heel pad in people with diabetes and healthy subjects. Strain in microchambers in people with diabetes was significantly greater than that of the healthy subjects [48]. In healthy subjects, macrochambers' strain was significantly greater than that of the people with diabetes, which can be a result of uneven distribution of the collagen fibrils in a diabetic heel [65].
Furthermore, Hsu and co-workers concluded that the heel pad tissue properties are altered heterogeneously in people with diabetes, indicated by an increased stiffness in macrochambers but a decreased stiffness in microchambers, which were attributed to a diminished cushioning capacity in diabetic heels [48]. The energy dissipation or hysteresis (area between loading and unloading in force-deformation graph) was also shown to be significantly higher in people with diabetes [38,65]. Furthermore, in people with diabetes the ability of recovering the shape of the heel pad after unloading reduces which can be linked to the increase in the amount of energy dissipation [62,65].
Chatzistergos et al. [66] investigated the correlation between the mechanical behaviour of the heel pad in type-2 diabetes and blood biochemical parameters such as triglycerides and fasting blood sugar (FBS). A medium strength positive correlation was found between the stiffness of the heel pad and FBS and a strong correlation was found between the triglycerides and stiffness. In addition, a strong negative correlation was found between the triglycerides and energy absorption during loading [66]. Chatzistergos and co-workers [66] concluded that people with type-2 diabetes and high levels of triglycerides and FBS are more likely to have stiffer heel pads that may hinder the tissue's ability to evenly distribute loads that makes the tissue more vulnerable to trauma and ulceration.

Quantifying mechanical behaviour
As mentioned before, heel pad provides a cushioning interface between calcaneus bone and the ground and has a natural function of shock absorption. The mechanical properties of the heel pad govern the force-deformation behaviour of the heel pad during heel strike and therefore these properties affect the loading on musculoskeletal system [67]. In order to investigate the force-deformation behaviour of the heel pad in in vitro, in situ and in vivo conditions, several mathematical models have been developed and utilised [36,40,42,53,65,[68][69][70][71][72][73][74]. Additionally, a number of FE analyses were used to quantify the mechanical behaviour of the heel pad [45,52,[75][76][77][78][79][80].

Mathematical models
A number of studies developed mathematical models to quantify the mechanical behaviour of the heel pad in vitro [35,41], while others utilised mathematical models to describe the in vivo mechanical behaviour of the heel pad [42,44,49,51,53,65,70,81,82]. Most of the studies measured the elastic behaviour of the heel pad and just a few of them represented the viscoelastic behaviour of the heel during dynamic loading.
Ledoux and Belvis [41] performed in vitro test on a freshly frozen cadaver's plantar soft tissue and used the following equation to characterise the mechanical behaviour of the plantar soft tissue: A e e s = - (7) The elastic modulus and absorbed energy of different plantar sites were calculated in which the subcalcaneal heel pad showed higher elastic modulus and absorbed energy compared to other sites of the plantar surface. This can be partly related to the fact that the subcalcaneal fat pad is designed to absorb energy during heel strike, while the fat pad underneath other plantar sites is mainly adapted to provide functions such as even pressure distribution.
In another study by Pai and Ledoux [39], quasilinear viscoelastic (QLV) model was used in two approaches: traditional frequency-sensitive and indirect frequency-sensitive: where G is the time-dependent function, τ is the relaxation time, σ e is the elastic stress, ε is the strain, ∂ σ e (ε) ∂ ε is the derivative of elastic strain over strain and ∂ ε ∂ τ is the derivative of strain over time. QLV theory normally assumes that the elastic and time-dependent properties can be separated and a linear combination of the elastic behaviour and viscous behaviour can describe the mechanical behaviour of the system.
The coefficients of stress-relaxation response for diabetic and non-diabetic subjects were compared and significant differences were found in the value of B between diabetes versus healthy subjects. However, no significant differences were found in viscous coefficients between diabetic and non-diabetic specimens. The lack of differences between diabetic and non-diabetic tissue was attributed to the changes at the structural level that have not been reflected effectively at the material level [38,39].
A number of studies characterise the mechanical behaviour of the heel pad during in situ tests [40,83,84]. Ker [83] employed a nonlinear equation to characterise the force-deformation behaviour of the heel pad. However, the stiffness values determined from this equation were dependent on the stage of loading cycle.
Although a number of studies utilised the in vitro and in situ data for mathematical modelling of the plantar soft tissue, in vivo assessment of the biomechanical behaviour of the plantar soft tissue has been a more common method for obtaining data for mathematical model.
In order to quantify the in vivo mechanical properties of the heel pad, the following function was utilised to represent the force-deformation data from in vivo test on the heel pad [49]: where a and b are the constants that are calculated from fitting the function to the in vivo data. The parameters were extracted from two groups of subjects with and without heel pain. Although the value of b was significantly lower in the group with heel pain compared to the group without heel pain, there was no significant difference in the value of a between the two groups [49]. While the a value that represents the slope of the force-deformation graph is more related to the stiffness, the value of b relates to the rate of the changes in the stiffness by loading and represents the curvature of the force-deformation graph.
Challis et al. [70] used the same formula as proposed by Ledoux and Belvis [41] for modelling the in vivo force deformation relationship of the heel pad during indentation. They compared the thickness, strain, energy loss and stiffness of the heel pad in cyclists versus runners. Although the heel pad stiffness was found to be significantly less in runners compared to cyclists, the percentage of absorbed energy was not found to be significantly different.
A general formula with separate terms for geometry and material parameters was utilised in a number of other studies in order to introduce coefficients which reflect the material characteristic of the heel pad.
where E is the elastic modulus of the heel pad, a is the radius of the indenter, ν is the Poisson ratio, h is the thickness of the soft tissue and K is the function of ratio of the radius of the indenter to the thickness of the tissue. Zheng et al. [51] calculated Young's modulus for different regions of the plantar soft tissue which was different from 40 to 50 kPa in healthy subjects, while the values were 160% more in average for the same sites of the plantar soft tissue in diabetic feet [51]. While there was no indication about how the K value is different between the two groups, limiting the maximum deformation to 10% of the initial thickness of soft tissue may have an influence on the calculated coefficients. As the heel pad shows nonlinear elastic behaviour, the coefficients may be different with higher deformation.
Chao et al. [44] used the same formulation (Eq. 10) to compare the elastic modulus between two different age groups. The air-jet indentation system was used along with non-contact optical coherence tomography in four loading-unloading cycles with 0.4 mm/s deformation rate. It was found that the modulus of elasticity under the second metatarsal head is significantly higher in older group compared to their younger counterparts.
Most studies concentrated on representing the force-deformation behaviour during loading, while one of the characteristics of viscoelastic behaviour is having a different force-deformation behaviour during loading and unloading. Hsu (11) where σ is stress, ε represents the strain, σ max is the maximum stress, ε max represents the maximum strain and α is the curvature constant that is different for loading and unloading graphs. Hsu et al. [65] compared in vivo data from diabetic and healthy subjects in which the curvature constant was significantly higher in diabetics compared to healthy participants during unloading; however, there was no significant difference in the α value during loading.
In addition to the importance of elasticity in the mechanical role of the heel pad, viscosity also plays an important role in dissipating energy and hysteresis. In order to identify the dissipating energy ratio, a number of studies added viscous term to the force-deformation formulas.
One of the in vivo models representing both elasticity and viscosity of the heel pad was proposed by Gefen et al. [42]. They proposed Kelvin-Voigt model in which elasticity behaviour was characterised by a linear spring and viscous behaviour was represented by a nonlinear damper.
. . E s e h e e = + (12) where σ represents stress, ε and ε • represents the strain and strain rate respectively, E is the elastic modulus and η is the viscosity parameter of the soft tissue. Digital fluoroscopy along with pressure plate was utilised to measure the pressure and deformation during walking.
The proposed model improved the predicted force-deformation behaviour of the soft tissue significantly by adding viscosity; however, the assumption of linear spring is an oversimplification of the behaviour of the soft tissue due to the fact that the quasi-static behaviour of the soft tissue still shows nonlinear force-deformation behaviour [20,34].
In a more comprehensive approach, Natali et al. [20] developed a constitutive mathematical model for the mechanical behaviour of the heel pad based on in vitro and in vivo tests. They considered hyperelastic and viscoelastic behaviour of the plantar soft tissue and developed stress-strain curve based on the second Piola-Kirchhoff stress tensor and used Miller-Young et al. [35] and Ledoux and Belvins's [41] in vitro data and also Zheng et al. [85] and Erdemir et al.'s [45] in vivo data in order to adapt the formula for the plantar soft tissue: where S was the second Piola-Kirchhoff stress tensor, S ∞ was the elastic stress when viscous condition is completely relaxed, q i was the viscous stress, C was right Cauchy-Green strain tensor, k v and C 1 relate to initial bulk modulus and shear stiffness, respectively, whereas r and α are hyperelastic coefficients of the soft tissue. γ ∞ and γ i relate to stress-strain history and τ i represents the relaxation time.
On the other hand, Sciume et al. [86] used a mathematical model, which was based on thermodynamically constrained averaging theory (TCAT) [87]. The soft tissue was modelled as a porous medium filled by an interstitial fluid.
Ultrasound indentation has also been frequently used to inverse engineer the material coefficients of heel pad and also to measure heel pad stiffness and energy dissipation [20,45,88]. However, there have been only a few studies that have used mathematical modelling to characterise the elastic and viscous behaviour of the plantar soft tissue. Naemi et al. [53] developed a mathematical model, which considered both elastic and viscous behaviour of the heel pad. They also modelled the nonlinear behaviour of the heel pad using nonlinear spring and nonlinear damper. They claimed that during quasi-static tests in which only the elastic component of the heel pad plays role in resistance to compression, strain-stiffening can be observed. Therefore, they used power function for depicting the elasticity component proposed by Scott and Winter [89]. This study proposed a method to quantify the force-deformation behaviour of the heel pad under compression in cyclic loading and took into account the nonlinear viscoelastic behaviour of the heel pad.
Energy input and energy return were derived as follows: By fitting the elastic and viscous energy densities to the data, Naemi et al. [53] showed that elastic energy density was much higher than absorbed energy density and elastic stress was significantly more than viscous stress. Elastic scaling factor and exponent were also 1.9 times and 14 times higher than viscous scaling factor, respectively. Despite the findings, they also reported that the deformation rates at which the tests were performed were much lower than the realistic deformations of the heel pad achieved during walking and they recommended testing the heel pad at more realistic strain rates to achieve realistic coefficients.
Although a significant number of studies utilised mathematical modelling in order to model elasticity and viscosity behaviour of the soft tissue at the heel, there has been a paucity of studies in which mathematical models were utilised for quantifying the viscoelastic characteristics of the plantar soft tissue at the forefoot. Although Hsu et al. [12,48,65] found the differences between macrochamber and microchamber behaviour of the soft tissue at the heel between healthy and diabetic subjects, no study has investigated the model parameters of different layers of this soft tissue.

FE models
As shown earlier, in pathological conditions such as diabetes there is an increased interest in investigating the mechanical characteristics of the human plantar soft tissue. In vivo observations indicating that plantar soft tissue properties change as a result of tissue damage [49] or diabetes [12,46,65,90,91] have highlighted the clinical relevance of plantar soft tissue biomechanics.
In this context, inverse FE analysis can be employed for the deterministic assessment of plantar soft tissue mechanical properties. FE analysis is a powerful numerical method that enables solving problems with complicated geometry, material properties or loading that cannot be approached using analytical solutions. In its direct application, FE analysis enables the calculation of the mechanical response to loading (e.g. internal stresses/strains) of deformable bodies that have known geometry and material properties. However, in inverse FE analysis the values of the mechanical properties that minimise the difference between the in vivo (experimentally measured) response to loading and the simulated (FE) one are calculated.
Two main in vivo material testing techniques have been used to inform inverse FE analyses: indentation [45,52,[75][76][77][78] and compression [79,80]. In both cases, the plantar soft tissue is compressed between a rigid loading surface and a bony prominence but in the case of indentation, the loading surface is significantly smaller than the plantar surface of the foot. In both cases, the applied force is measured using a load sensor and tissue deformation is either inferred from the displacement of the loading surface [75,92] or directly measured using medical imaging techniques such as ultrasound [45,52] or MRI [80]. These measurements enable the calculation of an indentation/compression force-deformation graph that describes the macroscopic mechanical response to loading of the tissue.
As a next step, a FE model simulating the same loading scenario is designed and used to numerically calculate the same force/deformation graph. At first, the material coefficients of the tissue are assigned random values (within predefined range of values) and the difference between the experimental in vivo graph and the numerical one is calculated. An optimisation algorithm is used to update the material coefficients of the plantar soft tissue in search of those values that minimise the difference between experiment and FE simulation.
In this context, Erdemir et al. [45] combined ultrasound indentation with FE modelling to calculate the material properties of plantar soft tissue on a subject-specific basis using an axisymmetric model of the heel. This technique was used to find the hyperelastic (Ogden 1st order) coefficients of the heel pads of two age-matched groups of diabetic and non-diabetic volunteers without however revealing any statistically significant difference between them. In order to improve the level of subject specificity of the model, Chatzistergos et al. [52] developed a 2D plane stress model of the heel from frontal ultrasound images that took into account the subject-specific geometry of heel pad. This modelling technique was later enhanced with a method for the automated generation of 3D models of the heel pad from ultrasound images [19].
The aforementioned studies [19,52] assumed a bulk plantar soft tissue; however, in a more comprehensive approach Petre et al. [80] differentiated between three layers of soft tissue: skin, fat pad and muscle. For this purpose, a 3D model of the forefoot was developed based on loadbearing MRI and an optimisation-based method was used to inverse engineer the material coefficients for all three different layers. The case of a single layer of soft tissue (i.e. bulk plantar soft tissue) was also considered indicating that simulating different layers affects the value of the calculated peak pressures; however, it was noted that the location where they appear does not get affected.
At this point, it has to be emphasised that there is a limit to the number of coefficients that can be calculated from inverse FE based on indentation or compression alone. To solve this limitation, Fontanella et al. [79] combined in vitro data with the information gathered from in vivo tests. They utilised the in vitro tests to estimate the values of all twelve coefficients for their visco-hyperelastic material model and then used in vivo compression tests to modify six of the coefficients on subject-specific basis.
Even though FE analyses have shed new light on plantar soft tissue biomechanics [35,45,79,[93][94][95][96][97], their actual contribution for the improvement of the diagnosis and treatment of the diabetic foot or other foot-related pathologies is limited [98]. This is mainly attributed to the difficulty of using FE analysis outside the research domain and particularly within the context of clinical practice [98]. Developing reliable subject-specific FE modelling techniques that are easy to use and not computationally demanding remains the key barrier for clinically applicable FE modelling.

Concluding remarks
The mechanical behaviour of plantar soft tissue shows viscoelasticity characterised by the reaction force being affected by the amount of deformation and deformation rate. The behaviour of the plantar soft tissue is highly nonlinear and shows the strain stiffening effect, with the stress-strain graph showing difference between loading and unloading. In a sense, viscosity causes the reaction force at the same deformation to be less during loading compared to the force during unloading and the difference that is caused by hysteresis increases with an increase in deformation rate. As a complex structure, the plantar soft tissue's mechanical behaviour shows a high degree of dependency on the method of testing evidenced by the fact that the mechanical properties of heel pad extracted from the in vivo and in situ results were observed to be different. This, for example, can be explained by the indication that structural factors such as heel pad thickness and geometry of the calcaneus have an influence on the heel pad behaviour. The models that represent the viscoelastic behaviour of heel pad are scarce and there are few that can fully justify the features of the mechanical behaviour of the plantar soft tissue in different testing scenarios. Despite these limitations, the model parameters that show the viscoelastic behaviour of tissue under load have the potential to be used to assess the mechanical behaviour of soft tissue under load with implication in identifying its malfunction as a result of disease or injury.