Stiffness of the human foot and evolution of the transverse arch

TL;DR

This study reveals that the transverse arch primarily determines foot stiffness, traces its evolution in hominins, and suggests it predates other bipedal adaptations, highlighting its importance in human locomotion evolution.

Contribution

It introduces a curvature-based measure of transverse arch stiffness and demonstrates its key role in the evolution of human-like foot stiffness before Homo.

Findings

Transverse arch accounts for most stiffness difference between humans and primates.

Human-like transverse arch appeared around 3.4 million years ago.

Foot stiffness predates other bipedal adaptations in human evolution.

Abstract

Foot stiffness underlies its mechanical function, and is central to the evolution of human bipedal locomotion. The stiff and propulsive human foot has two distinct arches, the longitudinal and transverse. By contrast, the feet of non-human primates are flat and softer. Current understanding of foot stiffness is based on studies that focus solely on the longitudinal arch, and little is known about the mechanical function of the transverse arch. However, common experience suggests that transverse curvature dominates the stiffness; a drooping dollar bill stiffens significantly upon curling it along the transverse direction, not the longitudinal. We derive a normalized curvature parameter that encapsulates the geometric principle underlying the transverse curvature-induced stiffness. We show that the transverse arch accounts for almost all the difference in stiffness between human and monkey feet (vervet monkeys and pig-tailed macaques) by comparing transverse curvature-based predictions against published data on foot stiffness. Using this functional interpretation of the transverse arch, we trace the evolution of hominin feet and show that a human-like stiff foot likely predates Homo by $\sim 1.5$ million years, and appears in the $\sim 3.4$ million year old fossil from Burtele. A distinctly human-like transverse arch is also present in early members of Homo, including Homo naledi, Homo habilis, and Homo erectus. However, the $\sim 3.2$ million year old Australopithecus afarensis is estimated to have possessed a transitional foot, softer than humans and stiffer than other extant primates. A foot with human-like stiffness probably evolved around the same time as other lower limb adaptations for regular bipedality, and well before the emergence of Homo, the longitudinal arch, and other adaptations for endurance running.

Figures

Tables

Table 1. Table 1: Table of estimated normalized curvature, and the estimated and predicted normalized stiffnesses. Statistics for all estimates were obtained using Monte Carlo simulations (see methods). Normalized curvature for human feet ( c ^ h subscript ^ 𝑐 ℎ \hat{c}_{h} ) are derived from published morphometric data (Fig.  4 and table  3 ). Based on these, predictions are generated for the normalized stiffness of the human foot ( K ^ h , pred subscript ^ 𝐾 ℎ pred \hat{K}_{h,{\rm pred}} ) using equation ( M.7 )). Independently, published data 6 , 7 are used to estimate the ratio of the normalized human foot stiffness ( K ^ h subscript ^ 𝐾 ℎ \hat{K}_{h} ) to that of C. aethiops ( K ^ c subscript ^ 𝐾 𝑐 \hat{K}_{c} ) or M. nemestrima ( K ^ m subscript ^ 𝐾 𝑚 \hat{K}_{m} ). The stiffness of the human foot with the plantar fascia transected ( K ^ h , pf − subscript ^ 𝐾 ℎ limit-from pf \hat{K}_{h,{\rm pf}-} ) is compared against that of M. nemestrima with the plantar fascia transected ( K ^ m , pf − subscript ^ 𝐾 𝑚 limit-from pf \hat{K}_{m,{\rm pf}-} ). Dimensional foot stiffnesses from published data 6 , 7 are listed in the table  2 .

Variable	mean	median	SD	1st quartile	3rd quartile
${\hat{K}}_h/{\hat{K}}_c$	2.61	2.38	1.16	1.80	3.16
${\hat{K}}_h/{\hat{K}}_m$	2.59	2.34	1.22	1.74	3.16
${\hat{K}}_{h,\mathrm{pf}-}/{\hat{K}}_{m,\mathrm{pf}-}$	2.06	1.87	0.98	1.39	2.51
${\hat{c}}_h$	15.58	15.11	5.63	11.65	19.00
${\hat{K}}_{h,\mathrm{pred}}$	2.04	1.86	1.09	1.26	2.62

Table 2. Table 2: Table of estimated stiffnesses from published load versus displacement data for humans 6 , and for C. aethiops and M. nemestrima . 7 In addition to stiffness of the intact foot (shown here as K □ subscript 𝐾 □ K_{\square} ) for the three species, we use stiffness of feet with the plantar fascia transected (shown here as K □ , pf − subscript 𝐾 □ limit-from pf K_{\square,{\rm pf}-} ). These estimates were obtained by digitizing the published plots 6 , 7 of load versus displacement.

Species	Foot condition	Variable	Estimated value
Homo sapiens	intact	$K_h$ (N/mm)	481
Homo sapiens	plantar fascia removed	$K_{h,\mathrm{pf}-}$ (N/mm)	369
Chlorocebus aethiops	intact	$K_c$ (N/mm)	132
Macaca nemestrina	intact	$K_m$ (N/mm)	165
Macaca nemestrina	plantar fascia removed	$K_{m,\mathrm{pf}-}$ (N/mm)	159

Table 3. Table 3: Foot morphometrics for extant species are modeled as Gaussian random variables. Mean values and standard deviations were obtained from reported values in the literature (methods for details). The important ratios ( w / L 𝑤 𝐿 w/L and t / L 𝑡 𝐿 t/L ) were verified to be accurate using measurements of representative samples. These Gaussian random variables are used in the Monte Carlo simulations. The morphometric variables used are the lever-length of the foot L 𝐿 L , width of the tarso-metatarsal articular region w 𝑤 w , dorso-plantar thickness of the third metatarsal t 𝑡 t , and torsion of the fourth metatarsal θ MT4 subscript 𝜃 MT4 \theta_{\rm MT4} ).

Species	$L$ (mm)		$w$ (mm)		$t$ (mm)		${\theta }_{\mathrm{MT4}}$ (deg)
	mean	SD	mean	SD	mean	SD	mean	SD
Homo sapiens	200	14.0	60.0	5.4	18.0	1.6	23.6	7.1
Chlorocebus aethiops	85.0	4.3	24.0	1.2	9.0	0.45	0	2.5
Macaca nemestrina	100	6.0	35.0	2.1	10.0	0.6	0	2.5
Pan troglodytes	130	13.0	52.0	5.2	13.0	1.3	0	2.5
Gorilla gorilla	176	17.6	72.5	7.3	16.0	1.6	2.2	1.5

Table 4. Table 4: Values of L 𝐿 L , w 𝑤 w , t 𝑡 t and θ MT4 subscript 𝜃 MT4 \theta_{\rm MT4} used in estimating the normalized curvature c ^ ^ 𝑐 \hat{c} of fossil samples. Variable names with a subscript h ℎ h refer to human values (e.g.  t h subscript 𝑡 ℎ t_{h} ), subscript p 𝑝 p to the chimpanzee (e.g.  w p subscript 𝑤 𝑝 w_{p} ), and subscript g 𝑔 g to the gorilla (e.g.  L g subscript 𝐿 𝑔 L_{g} ). These values are represented by Normal distributions as shown in table  3 . Variables in angled brackets, such as ⟨ t h ⟩ delimited-⟨⟩ subscript 𝑡 ℎ \langle t_{h}\rangle , refer to the mean value shown in table  3 . See methods for details of source materials.

Species	Specimen	${\theta }_{MT4}$ (deg)	$L$ (mm)	$w$ (mm)	$t$ (mm)	$w/L$	$t/L$
H. naledi	UW 101-1456	38.0	137.0	38.0	16.0	0.277	0.117
H. erectus	D2669, D4165	28.0, 29.0	–	–	17.0	$\frac{w_h}{L_h}$	$\frac{t}{⟨t_h⟩}\frac{t_h}{L_h}$
H. habilis	OH 8	25.0	112	44.0	–	0.393	$\frac{t_h}{L_h}$
Burtele	BRT VP2/73	26.5	–	–	13.3	$\frac{w_g}{L_g}$	$\frac{t}{⟨t_g⟩}\frac{t_g}{L_g}$
A. afarensis	AL 333-160	17.0	–	–	17.1	$\frac{w_g}{L_g}$	$\frac{t}{⟨t_g⟩}\frac{t_g}{L_g}$

Equations

$$\hat{K}=\frac{K}{K_{\mathrm{plate}}}.$$

$$\frac{{\hat{K}}_1}{{\hat{K}}_2}=\frac{K_1}{K_2}{\left(\frac{L_1}{L_2}\right)}^3{\left(\frac{t_2}{t_1}\right)}^3\frac{w_2}{w_1}.$$

$$\hat{c}=\frac{L^2}{{\ell }^2}=\frac{c{L}^2}{t}.$$

$$\hat{c}=\frac{{\theta }_{MT4}}{\left(w/L\right)\left(t/L\right)}.$$

$${\mathbf{S}}_{\mathrm{T}}(x,y)=\left(x,R_{\mathrm{T}}\sin {\theta }_y,R_{\mathrm{T}}\cos {\theta }_y\right)$$

$${\mathbf{S}}_{\mathrm{L}}(x,y)=\left(R_{\mathrm{L}}\sin {\theta }_x,y,R_{\mathrm{L}}\cos {\theta }_x\right)$$

$$\log {\hat{K}}_{\mathrm{pred}}=\{\begin{array}{ll}0,& \hat{c}<10,\\ \frac{3}{2}\log \left(\frac{\hat{c}}{10}\right),& \hat{c}>10.\end{array}$$

$$\frac{t}{L}=\frac{t}{⟨t_g⟩}\frac{t_g}{L_g},$$

$${\bar{a}}_{\alpha \beta }$$

$$a_{\alpha \beta }$$

$${\bar{b}}_{\alpha \beta }$$

$$b_{\alpha \beta }$$

$${\epsilon }_{\alpha \beta }$$

$${\kappa }_{\alpha \beta }$$

$${\epsilon }_{\alpha \beta }$$

$${\kappa }_{\alpha \beta }$$

$$\mathcal{E}=\frac{Et}{2\left(1-{\nu }^2\right)}\int \mathrm{d}A\left[\left(1-\nu \right)\mathrm{tr}\left({\mathbf{\epsilon }}^2\right)+\nu {\left(\mathrm{tr}\mathbf{\epsilon }\right)}^2+\frac{t^2}{12}\left(\left(1-\nu \right)\mathrm{tr}\left({\mathbf{\kappa }}^2\right)+\nu {\left(\mathrm{tr}\mathbf{\kappa }\right)}^2\right)\right],$$

$${\sigma }^{\alpha \beta }$$

$${\mu }^{\alpha \beta }$$

$$\delta \mathcal{E}$$

$$\mathbf{f}=f^{\alpha }{\partial }_{\alpha }\mathbf{S}+p\bar{\mathbf{N}}.$$

$$\mathbf{T}=T^{\alpha }{\partial }_{\alpha }\mathbf{S}+\mathcal{T}\bar{\mathbf{N}}.$$

$$\mathbf{m}=m^{\alpha }{\partial }_{\alpha }\mathbf{S}=\left(m_{\perp }{n}^{\alpha }-m_{\parallel }{\mathrm{l}}^{\alpha }\right){\partial }_{\alpha }\mathbf{S},$$

$${\bar{\nabla }}_{\alpha }{\sigma }^{\alpha \beta }+2{\bar{b}}^{\beta }{}_{\gamma }{\bar{\nabla }}_{\alpha }{\mu }^{\gamma \alpha }+{\mu }^{\gamma \alpha }{\bar{\nabla }}_{\alpha }{\bar{b}}^{\beta }{}_{\gamma }+f^{\beta }$$

$${\bar{\nabla }}_{\alpha }{\bar{\nabla }}_{\beta }{\mu }^{\alpha \beta }-{\bar{b}}_{\alpha \gamma }{\bar{b}}^{\gamma }{}_{\beta }{\mu }^{\alpha \beta }-{\bar{b}}_{\alpha \beta }{\sigma }^{\alpha \beta }-p$$

$${{\mu }^{\alpha \beta }{n}_{\alpha }{n}_{\beta }|}_{\partial \bar{\mathcal{S}}}$$

$${-\left({\bar{\nabla }}_{\alpha }{\mu }^{\alpha \beta }\right)n_{\beta }-{\partial }_s\left({\mu }^{\alpha \beta }{n}_{\alpha }{\mathrm{l}}_{\beta }\right)|}_{\partial \mathcal{S}}$$

$${\left[{\sigma }^{\alpha \beta }+\left(2{\bar{b}}^{\alpha }{}_{\gamma }-{\bar{b}}^{\alpha }{}_{\rho }{n}^{\rho }{n}_{\gamma }\right){\mu }^{\gamma \beta }\right]n_{\beta }|}_{\partial \bar{\mathcal{S}}}$$

$$K_{\text{mbox}plate}=\frac{F}{\delta }=\frac{Ew{t}^3}{4L^3(1-{\nu }^2)}.$$

$${\bar{\nabla }}_{\alpha }{\sigma }^{\alpha \beta }\sim \mathcal{O}\left[\frac{Et\delta }{R\ell }\right],2{\bar{b}}^{\beta }{}_{\gamma }{\bar{\nabla }}_{\alpha }{\mu }^{\gamma \alpha }\sim {\mu }^{\gamma \alpha }{\bar{\nabla }}_{\alpha }{\bar{b}}^{\beta }{}_{\gamma }\sim \mathcal{O}\left[\frac{E{t}^3\delta }{R{\ell }^3}\right],$$

$$r\sim \mathcal{O}\left[\frac{t^2}{{\ell }^2}\right].$$

$${\bar{b}}_{\alpha \beta }{\sigma }^{\alpha \beta }\sim \mathcal{O}\left[\frac{Et\delta }{R^2}\right]{\bar{b}}_{\alpha \gamma }{\bar{b}}^{\gamma }{}_{\beta }{\mu }^{\alpha \beta }\sim \mathcal{O}\left[\frac{E{t}^3\delta }{R^2{\ell }^2}\right],$$

$${\bar{\nabla }}_{\alpha }{\sigma }^{\alpha \beta }=0,$$

$${\bar{\nabla }}_{\alpha }{\bar{\nabla }}_{\beta }{\mu }^{\alpha \beta }-{\bar{b}}_{\alpha \beta }{\sigma }^{\alpha \beta }=0.$$

$$\ell =\sqrt{Rt}.$$

Full text

Stiffness of the human foot and evolution of the transverse arch

Madhusudhan Venkadesan m.venkadesan@yale.edu Department of Mechanical Engineering & Materials Science, Yale University, New Haven, CT 06520, USA

Marcelo A. Dias Present address: Department of Engineering, Aarhus University, Inge Lehmanns Gade 10, 8000 Aarhus C, Denmark School of Science, Aalto University, Espoo, FI-02150, Finland

Nordic Institute for Theoretical Physics (NORDITA), Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden

Department of Physics and Astronomy, James Madison University, Harrisonburg, VA 22807, USA

Dhiraj K. Singh Present address: Engineering Mechanics Unit, Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore, 560064, India. Collective Interactions Unit, OIST Graduate University, 1919-1 Tancha, Onna-son, Okinawa 904-0495, Japan

Mahesh M. Bandi bandi@oist.jp Collective Interactions Unit, OIST Graduate University, 1919-1 Tancha, Onna-son, Okinawa 904-0495, Japan

Shreyas Mandre shreyas_mandre@brown.edu School of Engineering, Brown University, Providence, RI 02912, USA

Abstract

Foot stiffness underlies its mechanical function, and is central to the evolution of human bipedal locomotion.1, 2, 3, 4, 5 The stiff and propulsive human foot has two distinct arches, the longitudinal and transverse.3, 4, 5 By contrast, the feet of non-human primates are flat and softer.6, 7, 8 Current understanding of foot stiffness is based on studies that focus solely on the longitudinal arch,9, 10, 11, 12, 13, 14 and little is known about the mechanical function of the transverse arch. However, common experience suggests that transverse curvature dominates the stiffness; a drooping dollar bill stiffens significantly upon curling it along the transverse direction, not the longitudinal. We derive a normalized curvature parameter that encapsulates the geometric principle15 underlying the transverse curvature-induced stiffness. We show that the transverse arch accounts for almost all the difference in stiffness between human and monkey feet (vervet monkeys and pig-tailed macaques) by comparing transverse curvature-based predictions against published data on foot stiffness.6, 7 Using this functional interpretation of the transverse arch, we trace the evolution of hominin feet16, 17, 18, 19, 20 and show that a human-like stiff foot likely predates Homo by $\sim 1.5$ million years, and appears in the $\sim 3.4$ million year old fossil from Burtele.19 A distinctly human-like transverse arch is also present in early members of Homo, including Homo naledi,20 Homo habilis,16 and Homo erectus.17 However, the $\sim 3.2$ million year old Australopithecus afarensis18 is estimated to have possessed a transitional foot, softer than humans and stiffer than other extant primates. A foot with human-like stiffness probably evolved around the same time as other lower limb adaptations for regular bipedality,3, 21, 18, 22 and well before the emergence of Homo, the longitudinal arch, and other adaptations for endurance running.2

The arched human midfoot acts as a stiff propulsive lever,10, 12, 23 and effectively transmits forces from the ankle to the forefoot (Fig. 1 a, b). The flat feet of non-human primates are significantly softer,8 and severely bend at the midfoot when the forefoot is loaded (midfoot break).10, 1, 4, 8 Similarly, studies using 3D X-ray imaging of healthy and flatfooted humans24 suggest that healthy feet are over 100% stiffer than flatfooted ones; when subjected to body weight loading, the midfoot deformation was two-fold greater for flatfooted subjects, and registered as an increased dorsiflexion of the hallux relative to the cuneiforms.24 The plantar fascia (Fig. 2) and longitudinal arch are hypothesized to underlie the higher stiffness of the arched human foot.9, 10, 11, 12, 13, 14 However, cadaver studies that transect the plantar fascia show a decrease of less than 25% in foot stiffness (table 2)6. The transverse arch emerges as an alternative source of stiffness. Even slightly curling a thin elastic structure in the transverse direction causes it to substantially stiffen (Fig. 1c); a principle that is employed from the engineering design of measuring tape to the handling of pizza. A similar effect is evident in discrete mechanical mimics of the foot (Fig. 1 d),25 and rayed fish fins.26

Direct stiffness measurement using cadaveric feet (summarized in table 2) show that human feet are nearly three-fold stiffer than that of the semi-terrestrial vervet monkey (Chlorocebus aethiops) or the pig-tailed macacque (Macaca nemestrina). Such a comparison between human feet and other primates is however confounded by overall size differences (Fig. 6). For example, a curved sheet of paper may be stiffer than a flat sheet simply by virtue of difference in thicknesses, and not necessarily due to its curvature. First consider how size affects the bending stiffness of a flat elastic plate made of a homogeneous material. The stiffness of a cantilevered flat plate15 of thickness $t$, width $w$, length $L$, Young’s modulus $E$ and Poisson’s ratio $\nu$ is $K_{\rm plate}=Ewt^{3}/(4L^{3}(1-\nu^{2}))$ (appendix A). Normalizing the stiffness of a cantilevered curved shell (Fig. 1 b) by the stiffness of a flat plate with identical material, length, width and thickness isolates the contribution of curvature. This normalized stiffness ($\hat{K}$) is given by,

[TABLE]

The ratio of the normalized stiffnesses of two curved shells of identical materials, but different lengths ($L_{1}$, $L_{2}$), widths ($w_{1}$, $w_{2}$) and thicknesses ($t_{1}$, $t_{2}$), quantifies the individual contribution of curvature differences, and is given by,

[TABLE]

By isolating the contribution of curvature, this ratio helps compare feet of different sizes by factoring the influence of overall size, under the hypothesis that the material properties are conserved. The human foot still remains over $2.5\times$ stiffer than those of C. aethiops and of M. nemestrina ($\hat{K}_{h}/\hat{K}_{c}$ and $\hat{K}_{h}/\hat{K}_{m}$ in table 1). Transecting the plantar fascia softens the human foot by $\sim 23$%, but only $\sim 4$% for the macacque, and therefore the human foot remains over $2\times$ stiffer than the macaque ($\hat{K}_{h,{\rm pf}-}/\hat{K}_{m,{\rm pf}-}$ in table 1). The contribution of the arched structure of the human foot to its stiffness is therefore substantially greater than can be accounted for by the longitudinal arch and the plantar fascia. How much of this difference is attributable to the transverse curvature versus differences in their material properties arising from soft tissues and muscle activity?

We first analyze the mechanics of shells to delineate the contribution of the transverse arch. Because transverse curvature $c$ couples longitudinal bending and in-plane stretching (supplement A), the mechanical response is governed by a balance in the elastic energy stored in bending versus stretching. Bending and stretching energies scale differently with thickness, as $t^{3}$ and $t$, respectively. Therefore, a new length-scale $\ell=\sqrt{t/c}$ emerges from this bending-stretching trade-off induced by the curvature $c$ (supplement A.2). When the shell experiences appreciable stretching, its stiffness scales as $K\sim Ewt^{3}/\ell^{3}$, unlike a flat plate that scales as $K_{\rm plate}\sim Ewt^{3}/L^{3}$ (appendix A.2). Their ratio, the normalized stiffness $\hat{K}$, depends only on the normalized curvature $\hat{c}$,

[TABLE]

In both experiments and computer simulations (Fig. 3a and appendices A–B), the normalized stiffness $\hat{K}$ is mostly governed by $\hat{c}$ and no other parameter, as evidenced by the nearly perfect collapse onto a single master curve (Fig. 3b).

There are two distinct regimes to the behavior of a shell; “low curvature” where $\hat{K}\approx 1$ is nearly independent of $\hat{c}$, and “high curvature” where $\hat{K}\sim\hat{c}^{3/2}$ (Fig. 3b). The transition between these regimes occurs around $\hat{c}=10$. The longitudinal arch however has no such effect on stiffness (Fig. 3b), highlighting the dominant role of the transverse arch. We note that continuum shell models are not an exact model of the foot (see appendix C.1 for validity of this model). Rather, the normalized curvature $\hat{c}$ parameterizes the strength of the coupling between out-of-plane bending and in-plane stretching, a principle that is purely geometric and independent of material differences.

How well does the normalized transverse curvature predict differences in the normalized stiffness between humans and monkeys? Estimating $\hat{c}$ requires morphometric measurements of all three dimensions of the foot (Fig. 4a, b). However, direct measurement in fossil specimens is often impossible because of their condition (e.g. Fig. 4d). We therefore derive a curvature estimate that relies on torsion of the fourth metatarsal $\theta_{\rm MT4}$, a widely reported measurement for both extant and fossil hominin feet,17, 18, 19 and a known quantitative measure of the transverse arch.3, 17, 18, 19 When the proximal metatarsal heads form a transverse arch and the distal metatarsal heads rest on the ground, the lateral metatarsals increasingly acquire torsion about their long axis (Fig. 4b, c). The torsion $\theta_{\rm MT4}$ arises from the curvature $c$ over the width $w$ of the tarso-metatarsal articulation, and therefore the curvature is given by $c=\theta_{\rm MT4}/w$. Using equation (3) we find,

[TABLE]

Normalized transverse curvature of human feet $\hat{c}_{h}$ is $15.6\pm 5.6$ (mean$\pm$standard deviation, table 1). The mean predicted normalized stiffness (equation M.7) for this normalized transverse curvature of the human foot is $K_{h,{\rm pred}}=2.04$ (table 1). Because monkey feet are flat, we set their normalized stiffness to 1, i.e. the hypothesis that the feet of monkeys would be no stiffer than a perfectly flat structure that is identical in all other regards. Under this hypothesis, experimentally estimated $\hat{K}_{h}=\hat{K}_{h}/\hat{K}_{c}=\hat{K}_{h}/\hat{K}_{m}\approx 2.6$. However, this estimate of $\hat{K}$ with intact feet includes the contribution of both the arches. To assess the individual contribution of the transverse arch, we compare stiffnesses after transecting the plantar fascia (table 2), and find that the experimental mean of $\hat{K}_{h,{\rm pf}-}=2.06$ agrees with the predicted mean $\hat{K}_{h,{\rm pred}}=2.04$ (table 1 and the middle 50th percentile in Fig. 3 b, pink band). We conclude that the transverse arch, characterized by its normalized curvature $\hat{c}$, dominates foot stiffness.

The estimated $\hat{c}$ for extinct and extant hominin feet suggest a gradual evolution of a human-like foot stiffness (Fig. 4d, e). As expected,9, 7, 5 the vervet monkey, macacaque, chimpanzee and gorilla feet are substantially flatter than humans, with no appreciable transverse arch, and well below the stiffness transition at $\hat{c}\approx 10$. Also as expected,2, 2, 3 the feet of fossils identified as genus Homo have a transverse arch that falls within normal human range, including Homo naledi,20 Homo habilis (OH-8),16 and Homo erectus (Dmanisi).17 Surprisingly however, the transverse arch of the $\sim 3.4$ million year old Burtele foot (unidentified species)19 also falls within normal human variation. It was previously hypothesized to possess a gorilla-like foot and gait, primarily based on the proximal articulation of the hallux, and the apparent lack of a longitudinal arch.19 Also surprisingly, $\hat{c}$ of the $\sim 3.2$ million year old Australopithecus afarensis (AL-333) falls outside human range, despite the fourth metatarsal torsion being well within human range.18 This is despite using a gorilla template14 to estimate overall foot length and width; a chimpanzee template would further separate it from humans.

Although the transverse arch underlies the ability of the foot to act as a rigid propulsive lever in walking, it does not imply all human-like capabilities such as endurance running. By virtue of its transverse arch, the Burtele hominin was likely to have possessed a foot capable of a human-like walking gait. The absence of a longitudinal arch (lesser elastic energy storage), long toes and an abducted hallux (stress injury risk) however indicate poor endurance running capabilities.2 The A. afarensis foot is likely much softer than Homo and Burtele, yet slightly stiffer than the gorilla by virtue of being close to, but below the transition region of $\hat{c}\approx 10$. It probably possessed a transitional gait; partly human-like and striding, but with greater midfoot flexibility. This is consistent with the 3.66 million year old Laetoli G footprints thought to have been made by A. afarensis, which indicate partially human-like foot function.21, 22 Our findings contest the current interpretation of the H. naledi foot.20 A weak or absent longitudinal arch has raised debate on the gait of Homo naledi although it has been previously noted that its metatarsal torsion is within human variation (table 4), as are several other pelvic, leg and foot adaptations.20 Based on its transverse arch, but a weak or absent longitudinal arch, we propose a human-like walking gait for H. naledi, but poor endurance running abilities.

The talotarsal joints form an elastic foundation in series with the stiff transverse arch in humans. Therefore, foot muscles that affect structural stiffness of the inter-tarsal joints may significantly modulate foot stiffness; an essential aspect of healthy foot function in humans.27 Cadaveric and in vivo evidence suggest such a role for the tibialis posterior muscle,28 which is phasically activated in mid-stance during walking.29 The tibialis posterior suddenly relaxes29 in late stance, thereby abruptly increasing foot flexibility,23 and recovering the elastic energy stored in the plantar fascia during touch-down.6 The transverse arch may therefore act as a mechanical clutch (activated by tibialis posterior) that engages or disengages the elastic recoil from the longitudinal arch. Such a role in elastic energy storage and recovery is particularly important for reducing the metabolic cost of running, and supports the hypothesis that the longitudinal arch is a later adaptation for endurance running.2 Finally, our study calls for reorienting clinical practice on foot function in walking towards the transverse arch, instead of the current focus on the longitudinal arch.

Acknowledgments

Funding support from the Human Frontier Science Program. Access to skeletal specimens were provided by Gary Aronsen, Kristof Zyskowski, Eric Sargis, Yale Biological Anthropology Labs and the Yale Peabody Museum. K.J. Meacham III provided experimental support.

Author contributions

M.V., M.M.B., and S.M. conceived of the study; M.A.D. and S.M. performed the mathematical modeling in consultation with M.V.; D.K.S. and M.M.B. performed the shell experiments; M.V. performed the biological data collection and analyses; M.V. wrote the main paper, S.M., M.A.D., M.M.B., and M.V. wrote the supplement, and all authors jointly edited the entire paper.

Methods

Numerical simulations

We simulated the elastic response of arched shells using the Shell interface in the 3D Structural Mechanics module of COMSOL Multiphysics v5.1 (COMSOL AB, Stockholm, Sweden). The transverse arch is represented by the map for the central plane of the shell given by

[TABLE]

where $\theta_{y}=y/{R_{\mathrm{T}}}$, and the longitudinal by

[TABLE]

where $\theta_{x}={x}/{R_{\mathrm{L}}}$, and in both cases $x\in\left[-L/2,L/2\right]$, $y\in\left[-w/2,w/2\right]$. For all the simulations, we set $L=0.1$ m and $w=0.05$ m. The material was modeled as linearly elastic with Young’s modulus $E=3.5$ MPa, Poisson’s ratio $\nu=0.49$, and mass density $\rho=965$ kg/m3 to match the material used in the experiments.

The boundary at $x=-L/2$ is clamped, i.e. zero displacements and rotations. The conditions at the other boundary $x=L/2$, are a uniform shear load $\mathcal{T}$, zero bending moment along $z$, and zero in-plane traction so that the displacements are free (see Figure 5 for axes orientations).

We solve this model for a range of thicknesses $t$, from 3 to 9 mm in steps of 1 mm, and transverse curvature radii $R_{\mathrm{T}}=0.03,~{}0.05,~{}0.07,~{}0.1,~{}0.3,~{}0.5,~{}0.7,~{}1,~{}3$m. For each combination of $t$ and $R_{\mathrm{T}}$, shear $\mathcal{T}$ ranging from 0 to 1 N/m is applied in increments of 0.1 N/m. The resulting out-of-plane displacement $\delta z$ is measured (Fig. 5(b)), and plotted against $\mathcal{T}$. The slope of these curves extrapolated to $\mathcal{T}=0$ yield the stiffness defined as $k\equiv w\,\mathcal{T}/\delta z$. The above process is repeated for geometry in Eq. (M.6).

Mechanical experiments

We fabricated and measured the stiffness of shells with a transverse arch or a longitudinal arch, and of a flat plate. These were all fabricated using polymer moulding techniques with PDMS (Poly dimethyl siloxane). The mould was fabricated using additive manufacturing (3D printed using ProJet 460Plus, 3D Systems). The printed mould was few millimeters in thickness, with one side left open. PDMS silicone elastomer (Sylgard 184, Dow Corning) was employed to cast the arch in the mould. Because the volume ratio of the base polymer to the curing agent controls the material bulk modulus for PDMS, the same ratio of 5 parts base polymer to 1 part of curing agent by weight was consistently maintained across all fabricated arches (appendix B). During an experiment, the fabricated arch was mounted on the experimental rig with help of clamps that were custom fabricated to exactly match the arch curvature. The clamps were additively manufactured (Stratasys Dimension 1200es) with ABSPlus (Acrylonitrile butadiene styrene) thermoplastic material (glass transition temperature 108*∘*C). One end of the clamped arch was affixed to a rigid plate attached to a free-moving horizontal translation stage to maintain a nearly zero horizontal load during the test. The other side of the clamped arch was connected to a force sensor, which was affixed to a vertical translation stage. The forces were measured using a data acquisition system (LabView, National Instruments) at 2 KHz for 1 second duration. The load test was performed under quasi-static loading of the arch sample by providing small displacements (quasi-static steps) of $5\times 10^{-5}$ m ($50~{}\mu$m) per step for a total of 10 quasi-static steps ($5\times 10^{-4}$ m or $500~{}\mu$m). Forces were measured after each quasi-static displacement. The slope of the force-displacement curve is the bending stiffness $K$ for the arch sample. Three experimental runs were conducted for each arch and their force-displacement curves were reproducible to within measurement error.

Monte Carlo simulations

Anatomical variability in the size of feet (table 3) is incorporated using Monte Carlo simulations to generate statistics for normalized stiffnesses and curvatures (table 1). The histograms generated from the Monte Carlo simulations are non-Gaussian. Therefore, the median and quartiles are reported in addition to the mean and standard deviation (SD). We used 1 million random combinations of the anatomical dimensions, where each dimension was drawn from an independent Gaussian distribution with means and standard deviations according to tables 3 and 4. Increasing the size of the Monte Carlo beyond a million samples had no effect on the statistics of the estimated quantities, for the number of significant digits reported. Our simulations likely overestimate the variance of relevant ratios such as $w/L$ and $t/L$ in comparison to biological feet, because we do allow for independent variation of all dimensions and do not incorporate covariation that may exist.

Foot stiffness

Published force versus displacement measurements of cadaveric human feet,6 and from monkeys (Chlorocebus aethiops and Macaca nemestrina)7 were used to estimate their stiffness. These published measurements6, 7 mimic the numerical and experimental design of stiffness measurements; the heel is clamped, the ball of the foot is vertically displaced with no horizontal forces, and the vertical reaction force is measured. The ratio of total measured force to the total applied displacement is the foot’s stiffness $K$. The published plots, figure 3 in Ker et al.6 and figures 1–2 in Bennett et al.7, were digitized to extract values (using B. Tummers, DataThief III. 2006 $<$http://datathief.org/$>$).

Normalization of the dimensional stiffness estimates of the human foot are carried out by considering the ratio of human to monkey foot stiffness. We report the normalized stiffness relative to both C. aethiops and M. nemestrina in table 2 according to equation (2). In addition to the ratio of stiffnesses of intact feet, we also estimate the ratio for feet where the plantar fascia was transected. This allows the quantification of the transverse arch’s individual contribution to foot stiffness.

Curvature-based stiffness prediction

The normalized stiffness $\hat{K}$ scales as the $3/2$ power of $\hat{c}$ for continuum shells for high curvature, but is independent of $\hat{c}$ and equal to 1 for low curvature. The equation to predict $\hat{K}$ based on these asymptotic behaviors is given by,

[TABLE]

Feet of extant species

Morphometrics for all species were obtained from the published literature. Additionally, for the feet of humans, C. aethiops, and M. nemestrina, we carried out primary measurements using specimens that were most similar in their lever length $L$ (Fig. 4) to the mean value reported in the literature. We carried out these measurements using software-based photogrammetry1 of high resolution images, and cross-verified with measurements using a digital caliper (0.01 mm resolution). The human specimen is from a 72 year old Caucasian male (Skulls Unlimited, #11177), the C. aethiops foot from the Yale Biological Anthropology Laboratory (YBL.3032a) and the M. nemestrina specimen from the Yale Peabody Museum (YPM MAM 9621). These samples were used to verify that the ratios $w/L$ and $t/L$ obtained from the literature are indeed representative of a sample of similar size.

The lever length $L$ follows standard definition to be the distance from the posterior end of the calcaneus to the anterior end of the third metatarsal’s distal head. The width $w$ is measured at the tarsometatarsal joint, from the most medial aspect of the distal articular surface of the navicular to the most lateral aspect of the distal articular surface of the cuboid. The thickness $t$ is defined as the dorso-plantar thickness of the proximal head of the third metatarsal, or the average of the second and fourth, when the third metatarsal data are unavailable.

Mean and standard deviation (SD) of the lever length $L$ were estimated from published data for humans,2, 3, 4, 5, 6 chimpanzees,3, 5, 6 gorillas,3, 5, C. aethiops,7, 5 and M. nemestrina.8, 9, 5 Mean $w$ is based on 3D anatomical data from BodyParts3D for humans,10 estimated from reported $w/L$ or dorsal skeletal views for chimpanzees and gorillas,5, 6, and primary measurements for C. aethiops and M. nemestrina (Fig. 6). To estimate the SD of $w$, we used reported variability in the medio-lateral width of the proximal metatarsal heads for all species17, 18, 19 to estimate the coefficient of variation (SD/mean), and applied that to $w$. The mean and SD of $t$ were all obtained from published values,18, 19 and confirmed with primary measurements for available specimens. Torsion of the fourth metatarsal $\theta_{\rm MT4}$ is used to estimate the transverse curvature, and published values were used for all species included in this study 17, 18, 19, 20. For species where the feet are regarded as flat, we set $\theta_{\rm MT4}=0$ deg.

Fossil feet

We use photogrammetry1 on published images of fossil feet (Fig. 4 d), and also data tables that accompanied the publication of these fossil data to estimate necessary dimensions and ratios.16, 17, 20, 18, 19

Among the fossil feet, all but the foot of Homo naledi20 were incomplete in some regard. For those incomplete feet, an extant species was selected as a template by taking into consideration published analyses of other postcranial and cranial elements. Based on this, Homo sapiens was chosen as the template for Homo erectus (Dmanisi)17 and Homo erectus (Olduvai hominin),16 and Gorilla gorilla was chosen as the template for Australopithecus afarensis (AL 333)18 and the unknown hominin foot found in Burtele.19 For example, the sole fourth metatarsal of the A. afarensis does not permit the direct estimation of $w$. However, only the ratio $w/L$ is necessary for the analyses, and the ratio for the gorilla is used for the Monte Carlo analysis of the fossil. The metatarsal however provides a direct measurement of $t$, but not of $L$. Therefore, to estimate the ratio $t/L$, we incorporate the measured thickness $t$ and the gorilla’s ratio $t_{g}/L_{g}$ by using the formula,

[TABLE]

where $\langle t_{g}\rangle$ is the mean $t$ for the gorilla. This template-based estimation therefore incorporates direct measurements where available without assuming that the fossil exactly resembles the extant template.

Curvature of biological feet

Transverse curvature $c$ is estimated using the torsion of the fourth metatarsal17 $\theta_{\rm MT4}$ as given by, $c=\theta_{\rm MT4}/w$, where $w$ is the width of the foot as measured at the tarso-metatarsal joints. Equation (4) for estimating the normalized curvature $\hat{c}$ follows from equation (3).

References for methods

• 1

Schindelin, J. et al.

Fiji: an open-source platform for biological-image analysis.

Nat Methods 9, 676–82 (2012).

• 2

Hawes, M. R. & Sovak, D.

Quantitative morphology of the human foot in a north american population.

Ergonomics 37, 1213–1226 (1994).

• 3

Wang, W. J. & Crompton, R. H.

Analysis of the human and ape foot during bipedal standing with implications for the evolution of the foot.

Journal of Biomechanics 37, 1831–1836 (2004).

• 4

Witana, C. P., Xiong, S., Zhao, J. & Goonetilleke, R. S.

Foot measurements from three-dimensional scans: a comparison and evaluation of different methods.

International Journal of Industrial Ergonomics 36, 789–807 (2006).

• 5

Schultz, A. H.

Relations between the lengths of the main parts of the foot skeleton in primates.

Folia Primatologica 1, 150–171 (1963).

• 6

Gomberg, D. N.

Form and function of the hominoid foot (1981).

• 7

Anapol, F., Turner, T. R., Mott, C. S. & Jolly, C. J.

Comparative postcranial body shape and locomotion in chlorocebus aethiops and cercopithecus mitis.

Am J Phys Anthropol 127, 231–9 (2005).

• 8

Hamada, Y.

Standard growth patterns and variations in growth patterns of the japanese monkeys (macaca fuscata) based on an analysis by the spline function method.

Anthropological Science 102, 57–76 (1994).

• 9

Sirianni, J. E., Swindler, D. R. & Tarrant, L. H.

Somatometry of newborn macaca nemestrina.

Folia Primatol (Basel) 24, 16–23 (1975).

• 10

Mitsuhashi, N. et al.

Bodyparts3d: 3d structure database for anatomical concepts.

Nucleic Acids Res 37, D782–5 (2009).

Appendix

Appendix A Mathematical modeling

While detailed treatments on thin shell theory can be found in standard texts 2, 3, 5, 6, we present an abridged version here to define precisely the variables that follow. We use the resulting partial differential equations to derive the scaling relation that describe the essential features of the elastic response.

A.1 Continuum analysis of thin shells

Let $\bar{\mathcal{S}}$ represent a collection of material points in the reference (undeformed) state of the mid-surface of the shell. These material points in $\bar{\mathcal{S}}$ are mapped in 3D Euclidian space by a three component vector function $\mathbf{S}:\bar{\mathcal{S}}\longrightarrow\mathbb{R}^{3}$, where $\mathbf{S}(x^{1},x^{2})=\left(S_{1}(x^{1},x^{2}),S_{2}(x^{1},x^{2}),S_{3}(x^{1},x^{2})\right)$ and $\left\{x^{1},x^{2}\right\}$ are the local Lagrangian coordinates on $\bar{\mathcal{S}}$. Therefore, an induced prescribed metric for the surface $\bar{\mathcal{S}}$ is calculated to be111Throughout this appendix we adopt the Einstein convention, where summation is implicit for repetitive indices when multiplying components of vectors and tensors. $\bar{\mathbf{a}}=\left(\partial_{\alpha}\mathbf{S}\cdot\partial_{\beta}\mathbf{S}\right)\mathbf{\mathbf{d}}x^{\alpha}\otimes\mathbf{\mathbf{d}}x^{\beta}$. The deformation from $\bar{\mathcal{S}}$ to a target (deformed) mid-surface $\mathcal{S}$ is given by a deformation map $\mathbf{F}:\bar{\mathcal{S}}\longrightarrow\mathcal{S}$, thus inducing target metric given by $\mathbf{a}=\left(\partial_{\alpha}\mathbf{F}\cdot\partial_{\beta}\mathbf{F}\right)\mathbf{\mathbf{d}}x^{\alpha}\otimes\mathbf{\mathbf{d}}x^{\beta}$. Therefore, the reference and target components of the metric tensors, also known as first fundamental forms, are respectively given by

[TABLE]

The set of cartesian vectors $\partial_{1}\mathbf{F}$ and $\partial_{2}\mathbf{F}$ ($\partial_{1}\mathbf{S}$ and $\partial_{2}\mathbf{S}$) span the tangent space to the mid-surface of the shell $\mathcal{S}$ ($\bar{\mathcal{S}}$). Thus, we may define the unit normals to the target $\mathcal{S}$ and reference $\bar{\mathcal{S}}$ surfaces, respectively, by $\mathbf{N}=\left(\det\mathbf{a}\right)^{-1/2}\partial_{1}\mathbf{F}\times\partial_{2}\mathbf{F}$ and $\bar{\mathbf{N}}=\left(\det\bar{\mathbf{a}}\right)^{-1/2}\partial_{1}\mathbf{S}\times\partial_{2}\mathbf{S}$. The second derivative of the mapping of the material points of both reference and target surfaces, projected onto the direction of their respective normal fields, yields the curvature tensor, or second fundamental form. Their components, respectively for reference and target shells, are given by

[TABLE]

We define three components of the displacement field, pure in-plane $u^{\alpha}(x^{1},x^{2})$ and out-of-plane $v(x^{1},x^{2})$ displacements, and they give us complete information of the deformation of the mid-surface of the shell through the function $\mathbf{F}=\mathbf{S}+\mathbf{U}$, where $\mathbf{U}(x^{1},x^{2})=u^{\alpha}(x^{1},x^{2})\partial_{\alpha}\mathbf{S}+v(x^{1},x^{2})\bar{\mathbf{N}}$ is the displacement vector in normal coordinates. In order to relate stresses and moments on the shell with the displacement of material points, we define the two local measures of deformation, in-plane strain, $\varepsilon_{\alpha\beta}$, as a measure stretching and curvature strain, $\kappa_{\alpha\beta}$, as a measure of bending. Their general definitions are

[TABLE]

For small displacements, Eqs. (A.11) become

[TABLE]

where $\bar{\nabla}$ is the covariant derivative with respect the reference shell 7.

Once the kinematic variables have been defined in Eq. (A.11), we may write the total elastic energy as a surface integral of an energy density which is entirely expressed as a function of material constants and a linear combination of the first two invariants of the variables in Eq. (A.11)222Here, we shall only consider linear constitutive laws, which justifies dropping terms that go with the third invariants $\mathrm{I\!I\!I}_{\boldsymbol{\varepsilon}}=\mathrm{det}\left(\boldsymbol{\varepsilon}\right)$ and $\mathrm{I\!I\!I}_{\boldsymbol{\kappa}}=\mathrm{det}\left(\boldsymbol{\kappa}\right)$., namely $\mathrm{I}_{\boldsymbol{\varepsilon}}=\mathrm{tr}\left(\boldsymbol{\varepsilon}\right)$, $\mathrm{I\!I}_{\boldsymbol{\varepsilon}}=\left(\left(\mathrm{tr}\boldsymbol{\varepsilon}\right)^{2}-\mathrm{tr}\left(\boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon}\right)\right)/2$, $\mathrm{I}_{\boldsymbol{\kappa}}=\mathrm{tr}\left(\boldsymbol{\kappa}\right)$, and $\mathrm{I\!I}_{\boldsymbol{\kappa}}=\left(\left(\mathrm{tr}\boldsymbol{\kappa}\right)^{2}-\mathrm{tr}\left(\boldsymbol{\kappa}\cdot\boldsymbol{\kappa}\right)\right)/2$. It is crucial to note that the trace operation $\mathrm{tr}$ on any matrix, e.g. $\mathbf{M}$, is here taken in a specific way that depends on the metric of the reference state in the following way $\mathrm{tr}\left(\mathbf{M}\right)\equiv\bar{a}^{\alpha\beta}M_{\alpha\beta}$. For homogenous and isotropic elastic shells, the total elastic energy takes the following form

[TABLE]

where $E$ is the Young’s modulus, $\nu$ the Poisson ratio, and $t$ the thickness of the shell. The terms in Eq. (A.13) proportional to $t$ compose the stretching energy, whereas the ones proportional to $t^{3}$ the bending energy.

Symmetric stress and moment tensors are derived from Eq. (A.13), as conjugates of the in-plane and curvature strains, respectively. For simplicity, we choose the following linear constitutive laws for the components of stress, $\sigma^{\alpha\beta}$, and moment, $\mu^{\alpha\beta}$,

[TABLE]

where $\nu$ is the Poisson’s ratio and $\delta\mathcal{E}/\delta f(x)$ stands for the functional derivative of Eq. (A.13) with respect to $f(x)$.

In the classical theory of shells 2, 3, 5, 6, the principle of virtual work allows us to derive the equations of equilibrium and boundary conditions. This principle is stated through the following relation

[TABLE]

where $\delta\Omega^{\alpha}\equiv\mathrm{e}^{\alpha\beta}\left(b_{\beta\gamma}\delta u^{\gamma}+\partial_{\beta}\delta v\right)$ is the measure of rotation of a shell element and $\mathrm{e}^{\alpha\beta}$ is the permutation symbol ($\mathrm{e}^{12}=+1$, $\mathrm{e}^{21}=-1$, and $\mathrm{e}^{\alpha\beta}=0$ if $\alpha=\beta$). The second line in Eq. (A.15) contains the following external forces and torques acting on the shell that balance with the stresses and moments:

• •

Load $\mathbf{f}\mathrm{d}A$ that acts on an area element of the reference mid-surface:

[TABLE]

• •

Load $\mathbf{T}\mathrm{d}s$ that acts on a length element along the boundary of the reference mid-surface:

[TABLE]

• •

Moments acting on the boundary of the mid-surface are of two kinds, namely twisting moment perpendicular to the boundary, $m_{\perp}$, and bending moment parallel to the boundary, $m_{\parallel}$:

[TABLE]

where the directions $\hat{\mathbf{n}}$ and $\hat{\mathbf{l}}$ are defined as the perpendicular outward and tangent to the boundary edge, respectively.

From Eq. (A.15), after carefully expressing the variations $\delta\varepsilon_{\alpha\beta}$ and $\delta\kappa_{\alpha\beta}$ through a fundamental variation of the embedding $\delta\mathbf{F}=\delta u^{\alpha}\partial_{\alpha}\mathbf{S}+\delta v\,\bar{\mathbf{N}}$ (details of this derivation is found in 6), we find the balance equations

[TABLE]

and boundary conditions

[TABLE]

where $Q$ is an effective shear force and $M^{\alpha}$ is an effective membrane force component. From this point on we assume that the body surface forces are zero, $f^{\alpha}=p=0$.

A.2 Scaling analysis

We first discuss scaling laws for a flat plate of length $L$, width $w$, and thickness $t$, Young’s modulus $E$, and Poisson ratio $\nu$. The in-plane stretching may be neglected for a flat plate under purely bending loads 4. This neglect can be justified as the simplification of (A.12a) to $\varepsilon_{\alpha\beta}=(1/2)\left(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha}\right)+\partial_{\alpha}v\partial_{\beta}v$. In this case, the deformations are forced by out-of-plane bending, and therefore, if $v$ scales as the applied displacement $\delta$, then $u_{\beta}$ scales quadratically with $\delta$ as $\delta^{2}/L$ and $\varepsilon_{\alpha\beta}$ as $\delta^{2}/L^{2}$. In addition, the curvature in this case scales linearly with $\delta$ as $|b_{\alpha\beta}|=|\partial_{\alpha}\partial_{\beta}v|\sim\delta/L^{2}$.

The total elastic energy $\mathcal{E}$ for deforming a plate is given by the sum of two terms, stretching and bending energies. The bending energy scales as $Et^{3}Lw/12(1-\nu^{2})\times\delta^{2}/L^{4}$, while the stretching energy can be neglected because it scales a factor $(\delta/L)^{2}$ smaller as $EtLw/(1-\nu^{2})\times\delta^{4}/L^{4}$. Under this assumption, which decouples the stretching from the bending, we may estimate the stiffness of a plate when a shear vertical force, $F$, is applied at its edge. Adding the work done by the shear vertical force, $F\delta$, in order to deflect the plate by an amount $\delta$, energy minimization yields the following scaling for the deflection: $v\sim 4FL^{3}(1-\nu^{2})/\left(t^{3}wE\right)$. The stiffness of the plate is

[TABLE]

We next consider a curved shell, where the stretching energy can no longer be neglected. Let $\ell$ denote the characteristic length scale of the deformation resulting purely from an external force to bend the shell out of plane. Then the scales of in-plane and curvature strains that establish in response to this forcing are $\delta/R$ and $\delta/\ell^{2}$ respectively. This yields, according to the linear constitutive laws in Eqs. (A.14), the following scaling for the stress and moment components: $\sigma^{\alpha\beta}\sim\mathcal{O}\left[Et\delta/R\right]$ and $\mu^{\alpha\beta}\sim\mathcal{O}\left[Et^{3}\delta/\ell^{2}\right]$, where we let $R$ be the smallest principal radius of curvature. If we compare the orders of the terms in Eq. (A.19a-A.19b), we have

[TABLE]

and the ratio of the latter terms to the former is

[TABLE]

Similarly, for Eq. (A.19b) we have that

[TABLE]

and the ratio of the latter to the former is also $r$.

The length scale over which a shell deforms is commonly much larger than the thickness of the shell, implying $r\ll 1$, and therefore we can approximate Eq. (A.19a-A.19b) by

[TABLE]

Balancing the first term in (A.25b), which scales as $\mathcal{O}\left[Et^{3}\delta/\ell^{4}\right]$, with the second term yields

[TABLE]

The applied shear stress $Q$ in Eq. (A.20b) scales as $\mathcal{O}\left[Et^{3}\delta/\ell^{3}\right]$. The stiffness of the shell, defined as the total force applied on that edge divided by the displacement, scales as

[TABLE]

Eqs. (A.26-A.27) are the main results of the scaling analysis verified using computations and experiments on a thin shell.

Appendix B Experimental measurement of stiffness of shells

B.1 Fabrication of Arches

Fabrication of arches started with preparation of a Stereolithography (STL) file of the arch mould on a computer aided design (CAD) software tool (Solidworks, Dassault Systemes). The arch length ($L$), its width ($w$), thickness ($t$), and the longitudinal ($R_{L}$) and transversal ($R_{T}$) radii of curvature formed the relevant input parameters to the STL file. This STL file, in turn formed the input to the 3D printer (ProJet 460Plus, 3D Systems) to print the relevant mould for the arch. A component infiltrate (Colorbond) was applied to impart additional strength to the part printed with a powder composite (Visijet PXL Core) combined with a binding agent (Visijet PXL). The powder composite based printing method was adopted for reasons of convenience. It offered fast printing with reasonable layer resolution of 100 micrometers, while maintaining minimal corrugated step features that are inevitable when printing curved objects. Furthermore, with a melting point of 1450 *∘*C, powder composite based moulds proved very stable when baking the elastomer based arches as discussed in the following.

The printed mould served as a shell, few millimeters in thickness, with one side left open. Poly (dimethyl siloxane) (PDMS) silicone elastomer (Sylgard 184, Dow Corning) was employed to cast the arch in the mould. The PDMS base polymer and the curing agent were thoroughly mixed in a centrifuge (Thinky ARE-310). Since the ratio of base polymer to curing agent controls the elastic modulus for PDMS, the same ratio of 5 parts base polymer to 1 part of curing agent by weight was consistently maintained across all fabricated arches employed in experiments discussed here. This mixture was poured into the printed mould and degassed for a 30 - 40 minute duration, to remove all trapped air bubbles. The degassed PDMS mixture in the arch mould was then transferred to an oven and baked at 75 *∘*C for 12 hours to obtain the experimental arches as shown in fig. 7a.

During an experiment, the fabricated arch was mounted on the experimental rig with clamps that were custom fabricated to exactly match the arch curvature. The clamps were 3D printed (Stratasys Dimension 1200es) with ABSPlus (Acrylonitrile butadiene styrene) thermoplastic material (glass transition temperature 108*∘*C). The arch clamps were designed such that the clamp spacing exactly matched the arch thickness to ensure the arch shape (cross-sectional area and radius of curvature) remained unchanged when tightening the clamps. Glue was applied between the arch and clamp to maintain a strict no-slip condition during loading.

B.2 Experimental setup and protocol

The arch with clamps was mounted on to the experimental rig (fig. 7a). One end of the clamped arch was affixed to a rigid plate attached to a free-moving horizontal translation stage to maintain constant horizontal distance between the clamps during the load test. This prevented the arch from changing its length under loading. The other side of the clamped arch was connected to a force sensor (SMT-1 S-Type Load Cell, Interface Advanced Force Measurement) affixed to a translation stage that moved along the vertical direction. The force sensor was interfaced with a data acquisition system (LabView, National Instruments) to obtain the force measurements from the load test.

Unlike computer simulations, experimental arches must contend with gravitational force owing to mass of the arch sample. Circumvention of gravitational effects therefore requires that arch samples cannot be too long or too thick, lest they slacken under their own weight in the central region between the two end clamps. Arch slackening causes a systematic error in the experimental determination of bending stiffness $K$, for then the initial displacement must first straighten the arch, and only then proceed to load it. The arches were therefore designed to span the range of values for $L^{2}/Rt$ through a combination of the fabrication parameters, while ensuring the length $L$ and thickness $t$ stayed within reasonable bounds to avoid slackening.

Prior to each experiment, the vertical translation stage was positioned to ensure the arch sample was perfectly horizontal. This position was set as the zero displacement point for the vertical translation stage. An initial force reading was then taken to measure the force offset value arising from weight of the clamps and the rigid metallic connector between the clamp and force sensor. This offset was then subtracted from the measured force.

The load test was performed under quasi-static conditions of the arch sample by providing small displacements (quasi-static steps) of $5\times 10^{-5}$ m ($50~{}\mu$m) per step for a total of 10 quasi-static steps ($5\times 10^{-4}$ m or $500~{}\mu$m). A force measurement was conducted after each quasi-static displacement with the LabView interface at 2 KHz sampling frequency for 1 second duration. The arithmetic mean over the 2000 force data points was recorded as the measured force for a given displacement. The slope of the force-displacement curve obtained from a linear fit of experimental data then provided the bending stiffness $K$ for the arch sample. Three experimental runs were conducted for each arch and their force-displacement curves were measured to be reproducible to within measurement error.

B.3 Calibration of PDMS material constants

Unlike crystalline solids 8, there exists no microscopic theory that permits calculation of material constants for amorphous media which include polymers, alloys, granular solids etc. Whereas part of the difficulty arises from the structural disorder which imparts an effective modulus to amorphous solids, material constants also sensitively depend upon the preparation protocol. In particular, the elastic constants and poisson ratio of PDMS based soft solids depend upon the base polymer to cross-linker ratio, curing temperature and baking time. This calls for calibration of material constants of PDMS solids used for curved shells in this study.

Flat PDMS test specimens were prepared using same preparation protocols as the PDMS curved shells used in this study, viz. 5:1 base to curing agent ratio was thoroughly mixed in a centrifuge (Thinky ARE-310), followed by removal of trapped air bubbles via degassing for 30 - 40 minutes, and then baked at 75 *∘*C for 12 hours. The dimensions of the PDMS calibration specimens were $7\times 10^{-2}$ m (length), $5\times 10^{-2}$ m (width), and $5\times 10^{-3}$, $7\times 10^{-3}$, and $9\times 10^{-3}$ m (thickness) in order to match the dimensions of curved PDMS shells used in experiments.

The fabricated calibration specimens were mounted onto the experimental rig using a clamping mechanism (3D printed from ABS plastic) as detailed in section B.2. Clamps on either side of the calibration specimen were mounted onto the experimental rig with one minor change in the setup. Instead of vertical loading perpendicular to horizontal plane of the flat sheet, a horizontal loading in plane along the direction of longitudinal axis of the calibration sheet was performed. The tensile loading yielded a simple stress-strain curve (see Fig. 7c,d). The Young’s modulus obtained from a linear fit (Fig. 7c) through the experimental data was $1.07\times 10^{6}$ N/m2 for calibration sheets of thickness $5\times 10^{-3}$ m.

Whereas the same calibration specimens were used, the experimental setup was further modified for measurement of Poisson ratio ($\nu$). A camera (Pixelink PL-B778, 5 Megapixels) with a macro lens (single focal length 106 mm) was mounted atop the experimental setup to record the lateral strain $\epsilon_{\perp}$ acting normal to the direction of the applied tensile loading. A longitudinal strain $\epsilon_{||}$ (absolute strain of $1\times 10^{-5}$ m) was administered and the lateral strain ($\epsilon_{\perp}$) obtained from processing the images. The slope obtained from a linear fit (Fig. 7d) of the lateral strain ($\epsilon_{\perp}$) versus longitudinal strain ($\epsilon_{||}$) curve yielded a poisson ratio of $\nu=$0.48, 0.47 and 0.49 for calibration specimens with thickness $5\times 10^{-3}$, $7\times 10^{-3}$, and $9\times 10^{-3}$ m respectively.

B.4 Data Analysis: Determination of Bending Stiffness $K$

A typical force-displacement curve from an experimental load test of the arch is shown in Fig. 7b. The data collected in discrete steps of $5\times 10^{-4}$ m ($500~{}\mu$m) were fit to a continuous curve in order to ascertain the arch bending stiffness $K$ from the slope of the force-displacement curve. In the shell theory formulation used in our predictions, linearization of the kinematic relations in equations (A.11), as in equations (A.12), and constitutive laws in equations (A.14) are assumed. Therefore, the theoretically predicted scalings arise within the linear response regime for force-displacement curves. The loading experiments therefore had to be conducted within the linear response regime. A verification of this requirement was performed using polynomial fits of the experimental force-displacement curves. To this end, polynomial fits to linear (solid line in Fig. 7b) as well as quadratic (dashed curve in Fig. 7b) order were performed. The linear fit yielded $F=-4.4\times 10^{-4}+1.4\times 10^{3}x$ and the quadratic fit yielded $F=-1.2\times 10^{-3}+1.4\times 10^{3}x-1.4\times 10^{4}x^{2}$. The zeroth order coefficients for both linear and quadratic fits being below the force sensor noise floor ($1\times 10^{-2}$ N), they are treated equivalent to zero. The two fits are in agreement at linear order, with very slow variation in quadratic coefficient with displacement $x$. Hence, the experiments are considered to be operating within the linear response regime. The bending stiffness $K$ is computed from the slope of force-displacement curves using linear fits to experimental data.

Appendix C Biological feet

C.1 Validity of the shell approximation

The shell may be decomposed in two parts: a soft part of length scale $\ell$ near the loaded edge as described by the scaling analysis in section A.2, and the rest of the shell of length $L^{\prime}$ which acts as an Euler beam. The approximation of the shell as an Euler beam is supported by the results of the computation (Fig. 5d, e) of the distribution of in-plane longitudinal stress $\sigma_{xx}$ throughout the shell. The neutral plane for the cross-section of a shallow shell whose central plane is approximated by $z=y^{2}/2R_{T}$ may be determined as

[TABLE]

which is plotted in figure 5(c). This expression for the neutral plane agrees with the zero-contour of the longitudinal stress, as shown in Figure 5(e). The second moment of the cross section of such a shell is

[TABLE]

The bending rigidity of the corresponding Euler beam $B_{n}=\dfrac{Etw^{5}}{720R_{T}^{2}}$, and the contribution to the stiffness of the shell is $\dfrac{3B_{n}}{L^{\prime 3}}=\dfrac{Etw^{5}}{240R_{T}^{2}L^{\prime 3}}$.

This Euler beam model loses validity near the loaded edge. In this region, which we call the edge-boundary layer, the deformation departs from the simple pattern required for the Euler beam approximation. The scaling estimates we derived in the previous subsection hold in this region. Based on those scaling estimates, the stiffness contributed by this region to the shell scales as $K$ calculated in equation (A.27) and occurs in series with the stiffness of the Euler beam. Since the softest spring determines the stiffness of springs in series, the edge-boundary layer determines the elastic response of the shell if

[TABLE]

To estimate $L^{\prime}$, we subtract $2\ell$ each to account for the edge-boundary layers at the loaded and the clamped ends, implying $L^{\prime}=L-4\ell$. If the criterion in (C.30) is satisfied, the bending deformation is localized to the edge-boundary layer.

Using Monte Carlo simulations, as described in the methods, we determine the distribution of $\eta$ in human feet. The first, second and third quartiles for the distribution of $\eta$ are 0.0033, 0.030, and 0.14 respectively. The 95${}^{\text{th}}$ percentile for the distribution is 0.65, based on which we conclude that for the vast majority of feet, $\eta\ll 1$.

References for the appendix

• 1
• Koiter 1966

Warner Tjardus Koiter.

On the nonlinear theory of thin elastic shells. i, ii, & iii.

Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen, Series B, 69(1):1–54, 1966.

• Koiter 1970

Warner Tjardus Koiter.

On the foundations of the linear theory of thin elastic shells.

Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen, Series B, 73(3):169, 1970.

• Landau and Lifshitz 1959

L. D. Landau and E. M. Lifshitz.

Course of Theoretical Physics Vol 7: Theory and Elasticity.

Pergamon Press, 1959.

• Fluegge 1973

Siegfried Fluegge.

Stresses in shells.

Berlin: Springer, 2nd ed., 1973.

• Niordson 1985

Frithiof I Niordson.

Shell theory.

Elsevier, 1985.

• do Carmo 1976

Manfredo Perdigao do Carmo.

Differential geometry of curves and surfaces.

Prentice-Hall Englewood Cliffs, 1976.

• Kittel 2004

C. Kittel.

Introduction to Solid State Physics, Chapter 3, 8th Edition.

Wiley, 2004.