Investigation of the dependence of joint contact forces on musculotendon parameters using a codiﬁed workﬂow for image-based modelling Journal of Biomechanics

The generation of subject-speciﬁc musculoskeletal models of the lower limb has become a feasible task thanks to improvements in medical imaging technology and musculoskeletal modelling software. Nevertheless, clinical use of these models in paediatric applications is still limited for what concerns the estimation of muscle and joint contact forces. Aiming to improve the current state of the art, a methodology to generate highly personalized subject-speciﬁc musculoskeletal models of the lower limb based on magnetic resonance imaging (MRI) scans was codiﬁed as a step-by-step procedure and applied to data from eight juvenile individuals. The generated musculoskeletal models were used to simulate 107 gait trials using stereophotogrammetric and force platform data as input. To ensure completeness of the modelling procedure, muscles’ architecture needs to be estimated. Four methods to estimate muscles’ maximum isometric force and two methods to estimate musculotendon parameters (optimal ﬁber length and tendon slack length) were assessed and compared, in order to quantify their inﬂuence on the models’ output. Reported results represent the ﬁrst comprehensive subject-speciﬁc model-based characterization of juvenile gait biomechanics, including proﬁles of joint kinematics and kinetics, muscle forces and joint contact forces. Our ﬁndings suggest that, when musculotendon parameters were linearly scaled from a reference model and the muscle force-length-velocity relationship was accounted for in the simulations, realistic knee contact forces could be estimated and these forces were not sensitive the method used to compute muscle maximum isometric force. (cid:1) 2018 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http:// creativecommons.org/licenses/by/4.0/).


Introduction
Musculoskeletal models are mathematical representations of the skeleton and musculature that through computer simulations allow the estimation of internal forces, such as muscle and joint articular forces, occurring within the human body during locomotion and not measurable without surgically-invasive procedures.
Various so-called ''generic" lower limb musculoskeletal models, developed mainly using cadaveric data, have been proposed (Carbone et al., 2015;Delp et al., 1990;Modenese et al., 2011). Easier access to medical images and the emergence of dedicated software has now made it possible to generate subject-specific models with a high level of anatomical personalization (Taddei et al., 2012). Valente et al. (2014), for instance, assessed the robustness of a methodology for creating image-based subject-specific models using a software called NMSBuilder, which was proposed as standard tool for implementing OpenSim model generation workflows (Valente et al., 2017). Scheys et al. (2009) validated software to map musculotendon geometries from atlases of children with cerebral palsy into new models. However, due to the complexity of the topic and the absence of a commonly accepted and established workflow, a comprehensive description of the methodologies used to generate image-based patient specific models is still lacking, potentially compromising reproducibility (Erdemir et al., 2016) and hindering adoption in clinical practice.
Since comprehensive datasets of patients with implanted instrumented prostheses were made available (Fregly et al., 2012), several image-based subject-specific models of adult patients have been created and used with established biomechanical computational tools (Damsgaard et al., 2006;Delp et al., 2007)  to estimate knee articular forces (Gerus et al., 2013;Marra et al., 2015). However, very limited studies have focused on the estimation of internal forces in children using a subject-specific approach. Previous pediatric subject-specific models have been mainly used for estimating musculotendon lengths and moment arms (Hainisch et al., 2012;Scheys et al., 2008), joint kinematics (Kainz et al., 2016;Scheys et al., 2011) and muscle function , whereas previous walking simulations employed models with limited personalization of musculotendon geometries  or estimated only the ankle joint contact forces (Hannah et al., 2017;Prinold et al., 2016). Nonetheless, image-based models seem highly desirable for investigating children gait biomechanics, due to the sensitivity of muscle force estimation to the muscle moment arms (Ait-Haddou et al., 2004;Herzog, 1992). The muscle moment arms significantly differ when estimated by scaling generic models using gait analysis data as opposed to subject-specific image-based models, both in healthy children and in children with cerebral palsy .
A critical aspect when creating image-based musculoskeletal models, is the identification of musculotendon parameters (namely, optimal fiber length l m o and tendon slack length l t s Þ. In fact, such parameters cannot be derived directly from magnetic resonance imaging (MRI), whereas in principle muscle strength can be estimated using segmented muscle volumes or regression equations (Handsfield et al., 2014). Since previous extensive sensitivity studies, performed using generic models (Carbone et al., 2016;Navacchia et al., 2015;Scovil and Ronsky, 2006) showed these parameters can strongly affect muscle force generation, their effects within a subject-specific modelling approach should be assessed, especially if a complete modelling workflow has to be codified. This paper aims to: (a) codify the generation of juvenile subjectspecific musculoskeletal dynamics models, starting from medical imaging and motion capture data; (b) assess different methods for the identification of musculotendon architectural parameters to identify the most suitable techniques to include in our workflow.

Participants
To investigate the effects of modelling the muscle dynamics without confounding factors coming from pathological conditions, we focused on healthy children. The complexity of the imaging and instrumental protocol, however, made it unethical to run it in a juvenile population without a clinical justification. Hence we used data from a group of children with a previous history of Juvenile Idiopathic Arthritis (JIA) (Ravelli and Martini, 2007), who were enrolled within a larger clinical study run as part of the FP7 research project MD-Paedigree in various European pediatric hospitals. The study design involved three visits six months apart (M00, M06, M12), when a child's health status was assessed using clinical regional MRI scans (Malattia et al., 2013), ultrasound joint examination (Lanni et al., 2016;Roth et al., 2016), clinical gait analysis (Davis et al., 1991), and disease-specific clinical activity scores, e.g. JADAS (Consolaro et al., 2009). A child was deemed healthy at the time of the assessment if they were pain-free, showed no clinical or ultrasound signs of JIA activity, and no radiological signs of bone erosion or cartilage structural damage. Following these criteria, data from eight children (11 visits, 8 MRI and 11 clinical gait analysis datasets; 4 males and 4 females, age and anthropometry at each visit reported in Table 1) collected within MD-Paedigree for other purposes, were re-used for this study, with the approval of hospitals' research ethics committees.

Imaging and gait analysis data
MRI was used to acquire images of the lower limbs of each subject at the M06 visit. A 3D T1-weighted fat-suppression sequence (e-THRIVE) scanning sequence with 1 mm in-plane resolution and a 1 mm slice thickness was used to scan the entire lower limb in 5 stacks. This sequence allowed a scanning duration of around 10 min and hence favored pediatric applicability. The imaging protocol included positioning some MRI-visible skin markers on welldefined bony landmarks to allow registration between MRI and stereophotogrammetric data ( Fig. 1 A and B). Gait analysis data were collected using a modified Vicon Plug-in-gait markers set ( Fig. 1 B), employing an 8-camera stereophotogrammetric system (Vicon Motion System Ltd, Oxford, UK; 200 Hz), and two force plates (AMTI OR6-6; 1000 Hz).

Musculoskeletal models
A statistical shape model approach developed after Steger et al. (2012) was used to segment the MRI images and produce geometrical models of the lower limb bones and external skin surface. Meshlab (Cignoni et al., 2008) (http://www.meshlab.net) and Netfabb Basic (https://www.autodesk.com) were then used to group the bones and soft tissues belonging to each body segment Table 1 Anthropometrics of participants as measured at the control visits where they were considered in healthy conditions. If the visit did not correspond with the full lower limb MRI data collection, variations of anthropometry, as recorded in the gait lab, are reported. For subjects whose mass varied by more than 3%, the inertial properties of the segments were updated. For subjects whose height varied by more than 3%, the model produced from MRI was uniformly scaled. Patient controls, i.e. the planned visits, were M00 (baseline), M06 (six month follow up) and M12 (12 months follow up).

Subj
Gender ( Fig. 1 B). Inertial properties were calculated for each body segment using NMSBuilder (http://www.nmsbuilder.org), assigning different densities to bone (1.42 g/cm 3 ) and soft tissue (1.02 g/ cm 3 for females, 1.03 g/cm 3 for males) (Dumas et al., 2005;White et al., 1987). The lower limb skeletal model was created combining the functionalities of NMSBuilder and OpenSim. Joint parameters were identified by fitting analytical shapes to articular surfaces identified on the bone geometries in Meshlab. A sphere was fitted to the femoral head to model the hip joint as a ball and socket joint (3 Degrees Of Freedom, DOF), while cylinders were fitted to the posterior articular surface of the femoral condyles (Yin et al., 2015) and to the talar trochlea (Siegler et al., 2014) to identify the axes of hinge joint (1 DOF) for the tibiofemoral and talocrural joints, respectively. The models included a patella rigidly attached to the tibia, rotating around the knee axis. The axis of the subtalar joint (1 DOF hinge) was estimated joining the centers of two spheres fitted to the talocalcaneal and talonavicular articular surfaces (Parr et al., 2012). All fitting operations were least-squares procedures implemented in MATLAB (v8.5, R2015a, Mathworks, USA); the cylinder fitting relied on the LSGE MATLAB library (http://www.eurometros.org). Including a free joint between pelvis and ground (6 DOF), the bilateral model presented 18 DOF in total.
A detailed step-by-step description of the procedure is made available (see Appendix A).

Musculotendon geometry
The paths of 42 musculotendon actuators were defined through a supervised atlas registration procedure. For each segment of the generic model gait2392 (Delp et al., 1990), related muscle attachments (including via points) were grouped into an atlas. A set of well-defined (van Sint Jan, 2007), repeatable (Ascani et al., 2015) bony landmarks, constituting a landmark cloud, were then identified on the bone surface. Finally, these landmark clouds were reg-istered onto the subject-specific bones using affine transformations (Horn, 1987) to match correspondent, virtually palpated  bony landmarks. The subjectspecific muscle attachments and via-points resulted from applying these transformations to the points in the muscle attachment atlases. Minor adjustments to the muscle paths, performed consulting sectional anatomy textbooks, e.g. Moller and Reif (2007), were usually necessary to finalize the muscle paths. The path points variability due to these adjustments was estimated in the order of 3.17 ± 2.16 mm in a pilot study where a single operator created, for three times, three MRI-based models from datasets not used in this study. The peroneus tertius was the only muscle from gait2392 not included in the model as not identifiable in the MRI. For six models built using the gait data collected when the MRI was not available, a uniform scaling procedure was applied using the anthropometric information available from the gait lab (Table 1).
All the above operations, performed in NMSBuilder (Valente et al., 2017), allowed the generation of an OpenSim model from each MRI dataset (Fig. 1 C).

Musculotendon architectural parameters
Musculotendon units were modelled using a dimensionless three element Hill-type muscle model (Thelen, 2003;Zajac, 1989), which requires the definition of five parameters: the opti-  (Delp et al., 1990), and M2, that aims to reproduce in the subject-specific model, for all combinations of joint angles in the sagittal plane, the same maximum isometric muscle contractile (C) OpenSim model including musculotendon actuators and markers used for marker placement and for calculating joint angles. A marker placer procedure is used to place HEE in the OpenSim model by registering the markers available in the MRI with the static trial. MFC and D5M are not used for calculating the joint angles in dynamic trials. Please note that a metatarsophalangeal joint was included in the model, defined as a hinge with axis aligned with the line connecting the 1st and 5th distal metatarsal heads, but was not employed in the simulations for this study.
conditions of the gait2392, i.e. the same normalized fiber and tendon lengths Winby et al., 2008 Weijs and Hillen (1985)). In Eq. (4), V SS is computed from the linear regressions of Handsfield et al. (2014) and used to estimate F iso from the muscle physiological crosssectional area (Lieber and Friden, 2000) neglecting the effect of muscle pennation angle, accounted for in the computational muscle model. The combinations of M1-M2 and F1-F4 led to eight different estimates of the muscle parameters. Finally, v max was set to 10 fibers per second (Zajac, 1989) and a o values were taken directly from the gait2392 model.

Biomechanical simulations
The biomechanical simulations were performed in OpenSim 3.3  leveraging the MATLAB API and using single limb subject-specific models (11 DOF). Joint angles and moments were computed with the inverse kinematics and inverse dynamics OpenSim tools, respectively. The body and joint reference systems, based on ISB conventions (Wu et al., 2002), were defined using the bone geometries, whereas joint kinematics was computed tracking skin markers. This introduced an angular offset that depended on the joint axis definition, the subject pose during the MRI scans, and the location of the virtual markers in the model and the experimental skin markers. In spite of the impossibility of isolating the above different sources of variation, to facilitate -results' interpretation and make data from different models directly comparable, we decided to consider as zero the static posture joint angles. Joint powers were computed in MATLAB using angular velocities and inverse dynamics output. Muscle forces were estimated using a static optimization approach that, considering the physiological muscle force-length-velocity (F-L-V) relationship that regulates muscle force generation (Zajac, 1989), minimizes the sum of muscle activations squared. Ideal moment generators (reserve actuators), providing joint torque when muscle forces could not balance the external moments, were included for each DOF. Previous literature showed that the F-L-V relationship does not affect the static optimization results for a generic model (Anderson and Pandy, 2001). To verify whether this also applies to a subjectspecific case, simulations were performed both with and without considering the F-L-V relationship. To prevent muscle recruitment alterations, reserve actuators were made unfavourable to recruit by assigning them a unitary maximum force, i.e. a provided joint torque of 1 Nm weighted on the static optimization objective function as much as a fully activated muscle. Finally, the ''JointReaction" analysis in OpenSim (Steele et al., 2012a) was used to calculate the joint contact forces.

Quality assurance of simulations
The quality of the registration between the stereophotogrammetric landmarks trajectories to the subject-specific anatomical model derived from medical imaging was assessed by the peak tracking error over the entire markers' set, which, in agreement with best practices (Hicks et al., 2014), was ensured to be <20 mm for all the models.
The static optimization approach was verified by ensuring the contribution of reserve actuators to the joint dynamic equilibrium never exceeded 10% of the total joint moment over the entire gait cycle. This threshold was consistent with a previous study on a pediatric population (Steele et al., 2012b). The percentage of simulations discarded because of exceeding this threshold was then used as a metric to compare different modelling methods.
Additional quality checks to be performed on the experimental data, surface fitting and generated models are reported in the supplementary materials (step-by-step guide).

Statistical analysis of the results
Groups of simulation results were compared using Statistical Parameter Mapping (SPM). SPM relies on the Random Vector Field theory to account for data variability, allowing the comparison of group means for time-varying quantities. All analyses were done with MATLAB, using the spm1d package (Pataky et al., 2013).
The effect of the different methods used to estimate the maximum isometric force on the variability of the joint contact forces was evaluated with the SPM equivalent of a one-way ANOVA; the significance of individual differences between the joint reaction estimates was tested using the SPM equivalent of an unpaired Students' t-test. In all analyses the significance level was set to a = 0.05.

Results
Two of the models had to be scaled to account for height changes and four to account for mass variations (details in Table 1) in the absence of MRI data at that particular control. In total, gait from 11 control visits were simulated (107 trials: speed: 1.08 ± 0.07 m/s, step length 1.18 ± 0.05 m, cadence: 55.1 ± 2.7 stride/ min). Sagittal plane joint kinematics, moments and powers are presented in Fig. 2, out of sagittal plane results are available as supplementary materials.
Musculotendon parameters estimated using M1 resulted in consistently longer l t s and shorter l m o than those obtained by M2 (Fig. 3). The first three estimates of F iso;SS (F1-F2-F3) did not depend on musculotendon parameters, with F1 estimating the smaller F iso;SS and F3 the larger. F4 yielded the most variable F iso;SS because, despite estimated muscle volumes depended on M and H, large differences in l m o between M1 and M2 propagated to muscle strengths (Fig. 4).
As expected, no method achieved a 100% success rate in the verification of the muscle static optimization. Percentage of successful simulations are reported in Table 2 (see supplementary material for reserve moments). Hip internal/external rotation moment required intervention of reserve actuators most frequently, causing the majority of failures, followed by the ankle moment. M1 led up to 63-88% of successful simulations, while M2 (5-35% of successful simulations) resulted in models generally unable to provide the required ankle joint moment and was not further analyzed. The percentage of successful simulations decreased with decreasing estimation of F iso;SS , while neglecting muscle F-L-V relationship entirely allowed more than 92% successful simulations. Results from M2 successful simulations are available in supplementary materials.
Joint reaction forces are presented for hip, knee and ankle joints in Fig. 5, combining all approaches to estimate musculotendon parameters and maximum isometric force, while peak values are summarized in Table 3. When the muscles' F-L-V relationship was not considered, the second peak of knee contact forces was larger than that computed accounting for F-L-V relationship (p < 0.003), reaching mean magnitudes of up to 3.2 times the subject's body weight (BW) for F1, F2 and F3 (see Table 3). Conversely, estimated joint contact forces at hip and ankle joints were found insensitive to the method used to estimate F iso;SS throughout the walking stance phase (no significant differences, p > 0.05, see web supplementary material).

Discussion
The first aim of this paper was to codify the generation of subject-specific musculoskeletal dynamics models, starting from medical imaging and motion capture data. We achieved this by proposing a methodology that combines various modelling approaches (Parr et al., 2012;Prinold et al., 2016;Valente, 2013;Valente et al., 2017;Yin et al., 2015) and enables a complete musculoskeletal analysis, making this the first study to present a comprehensive subject-specific model-based characterization of the gait biomechanics of a juvenile population, including profiles of internal forces computed from musculoskeletal models. To favor reproducibility and adoption of the proposed methodology, the atlas of muscle attachments and landmark clouds used in the registration, a detailed step-by-step description of the procedure, and the models used in this study are publically available (see Appendix A). It is worth noting that the described methodology relies entirely, with the exception of MATLAB, on freely available software. Whilst the examined cohort is too limited for the reported curves to be considered as normality curves, they can still be used by others as a reference to detect possible deviations from healthy patterns. In this respect, a shortcoming of the study is that participants were classified as healthy from a cohort of children previously affected by JIA. Although they were objectively assessed by clinical experts and none of them exhibited signs of persisting disease activity, minor motor compensation or recovering mechanisms could have still occurred.
In the proposed workflow, the knee and talocrural joints are modelled as hinges. More complex kinematic models (Parenti-Castelli and Sancisi, 2013) would have required identification of ligaments attachments, and therefore more detailed medical images to be personalized (Brito da Luz et al., 2017). The patellofemoral joint is also simplified compared to other approaches, e.g. Sancisi and Parenti-Castelli (2011). However, since no muscle path penetrated the bone geometries, the chosen simplification was considered satisfactory for the purposes of this study. Finally, the proposed multi-segmental ankle joint complex presents a stronger link to functional anatomy (Isman and Inman, 1968;Kapandji, 1987;Parr et al., 2012) than other subject-specific modelling approaches (Prinold et al., 2016;Saraswat et al., 2010) relying on models inspired to gait analysis (Stebbins et al., 2006), so enabling investigation of the subtalar joint kinematics (Montefiori et al., 2017).
The second aim of the paper was investigating the effect of different methods for the identification of the musculotendon architectural parameters on the models' predictions. Our findings confirmed that appropriate estimation of musculotendon parame-ters are essential for obtaining successful simulations. Method M2, which successfully predicted muscle parameters consistent with cadaveric measurements in a previous study , was too sensitive to the differences between the generic and subject-specific models in terms of geometry of musculotendon paths, joint models and joint angle offsets, which are all elements compromising the resulting muscle function in the subject-specific model. This led to decreased muscle moment gen-  (1)-(4)). The estimations for method F4 depend on optimal fiber length, so they have been calculated for both M1 and M2. eration capacity and a small percentage of successful simulations. M1, however, yielded up to 88% of successful simulations.
Interestingly, the method chosen for estimating muscle strength did not affect the estimated joint contact forces, presenting magnitudes aligned with previous subject-specific simulations (Prinold et al., 2016;Valente et al., 2014). Nonetheless, this choice can strongly influence the percentage of successful simulations.
Conversely, when neglecting the muscle F-L-V relationship almost all simulations were successful, but the computed second peak of knee contact forces was significantly larger than that estimated using M1. This peak reached mean values of 3.2 BW, i.e. around 0.5 BW larger than those measured in adults using instrumented prostheses (Fregly et al., 2012;Kutzner et al., 2010). This difference might be due to larger forces produced by the rectus femoris and Table 2 Percentage of simulations for which reserve actuators contribute to the joint balance with <10% of the peak moment, reported for each degree of freedom of the model. When reserve actuators were below the 10% threshold at all joints, the simulation was considered successful. Note that maximum isometric forces for method F4 can be estimated only if musculotendon parameters have been computed, so only results considering the force-length-velocity relationship are presented. gastrocnemii (Figs. 6 and 7) when the F-L-V relationship is neglected, i.e. neglecting the effect of the muscle contraction dynamics. The proposed modelling methodology has some limitations. First, it relies on segmented bone geometries, and MRI segmentation is a time-consuming processing bottleneck. Nonetheless, promising developments in statistical shape modelling techniques (Zhang et al., 2014) and atlas based segmentations (Kolk et al., 2015) have been published recently. Excluding segmentation, an expert operator required around 10 h to generate a complete bilateral lower limb model, against 8-10 h reported by Prinold et al. (2016) for a foot and ankle model using similar procedures. Second, the supervised atlas registration is intrinsically limited by the fact that the atlas bony landmark clouds are identified on the simplified bone geometries of the gait2392 model. Whereas Valente et al., (2014) demonstrated the robustness of model output to variation in muscle attachments, Hannah et al. (2017) reported up to 24% of inter-operator variability for a foot and ankle model. Future work is needed to further automate musculotendon path generation and assess the repeatability of the procedure and the sensitivity of the model output to these variations. Additionally, future studies might focus on full validation and quantification of the accuracy of the model, possibly relying on data from instrumented prosthesis. Third, a direct comparison of predicted muscle activation against electromyographic signals could not be performed with the available data. As a consequence, M1 was preferred to M2 based on percentage of successful simulations, although both methods produced realistic results, with M2 yielding hip joint reactions closer to values measured in vivo (Bergmann et al., 2001) (see web supplementary results). An enhanced formulation of M2, robust against musculotendon path variations, could be preferable to M1 in the future. Finally, in the presented simula- Table 3 Peaks of joint reaction forces at hip, knee and ankle joint (mean ± standard deviation) reported as function of the method used to estimate the musculotendon parameters and the maximum isometric force. Note that the number of successful simulations was not consistent between simulation techniques. Method M2  was not considered because of the small percentage of successful simulation. Hip joint Method M1 3.8 ± 0.7 3.8 ± 0.7 3.8 ± 0.7 3.7 ± 0.8 No F-L-V relationship 4.0 ± 0.9 4.0 ± 0.9 4.0 ± 0.9 NA Knee joint Method M1 2.7 ± 0.7 2.7 ± 0.5 2.8 ± 0.7 3.0 ± 0.6 No F-L-V relationship 3.2 ± 0.9 3.1 ± 0.8 3.2 ± 0.9 NA Ankle joint Method M1 6.1 ± 1.3 6.2 ± 1.4 6.4 ± 1.4 6.3 ± 1.4 No F-L-V relationship 6.3 ± 1.5 6.3 ± 1.5 6.4 ± 1.5 NA Fig. 6. Muscle forces of vastus medialis, vastus lateralis, vastus intermedius and rectus femoris (rows) computed for the different methods used to estimate maximum isometric forces (columns) considering the force-length-velocity relationship with musculotendon parameters from method M1 (in black) and ignoring the force-length-velocity relationship (in red). Note that maximum isometric forces for method F4 can be estimated only if musculotendon parameters have been computed, so only results from simulations considering the force-length-velocity relationship are presented. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) tions, the muscle tendon was considered infinitely stiff and the passive muscle force was neglected. Consequently, reported findings might not extend to models including a compliant model of the tendon and considering passive muscle force.
In conclusion, this paper presented a methodology to create subject-specific models from medical images and applied it to a group of typically developing children. A complete analysis of their gait biomechanics, including muscle and articular forces that could serve as comparison for future studies on juvenile populations has been provided. By assessing different methods available in the literature, it was shown that musculotendon parameters can be estimated maintaining the ratio of tendon and fiber lengths from a generic model, maximum isometric force can be computed in different ways with little influence on final muscle force predictions and muscle force-length-velocity relationship should be considered in static optimization simulations to predict realistic knee contact forces. In the future, we plan to adapt the methods presented here to clinical populations, including children and adults with neurodegenerative diseases and musculoskeletal conditions such as osteoarthritis. Fig. 7. Muscle forces of gastrocnemius medialis, gastrocnemius lateralis and soleus (rows) computed for the different methods used to estimate maximum isometric forces (columns) considering the force-length-velocity relationship with musculotendon parameters estimated from M1 (in black) and ignoring the force-length-velocity relationship (in red). Note that maximum isometric forces for method F4 can be estimated only if musculotendon parameters have been computed, so only results from simulations considering force-length-velocity relationship are presented. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)