A publicly-available musculoskeletal model was expanded upon for this study (DeMers et al., 2017) (Figures 1A). Briefly, the whole-body model (49 muscles, 21 degree-of-freedom) contained elastic foundation foot-ground contact (Seth et al., 2018) and lumped passive ankle stiffness represented by an uncoupled, rotation-only bushing (Chen et al., 1988). Joint motion was driven by custom stretch-reflex feedback muscle controllers. For baseline validation, the generic model anthropometry was scaled to a single human subject (female, 68 kg, 180 cm), and muscle control parameters optimized to track whole-body joint kinematics measured via 3D motion capture during a 40 cm drop onto a level surface. The model also had a landing platform that could be inclined to induce virtual inversion sprains. The unmodified, subject-specific model was used as the mean baseline model for probabilistic simulations.

We also extended the baseline model with a representation of external, ankle-foot structural support (e.g., a brace, Figure 2). The rotation-only bushing representing passive anatomy was first duplicated, then scaled to represent adjustable magnitudes of bracing acting in parallel with passive anatomy. Specifically, angular displacement data was multiplied with a fixed factor (i.e., scaling flexibility (Wright et al., 2000a)), such that a brace with a 150% multiplier would be 50% less flexible than the model’s passive anatomy. This approach was chosen as ankle flexibility versus compliance (scaled torque at fixed displacements) –has been suggested to be more influential in sprain prevention (Wright et al., 2000a).

After verifying reproducibility of the baseline deterministic simulation (DeMers et al., 2017), we assigned probabilistic distributions to 10 input parameters expected to influence sprain occurrence (Table 1). Mean values were kept equal to the baseline model, while assumed variances were relatively minor perturbations aimed to preserve validity of the baseline model (Table 1). Monte Carlo (MC) simulations with 1000 trials were instantiated with random sampling (NESSUS, v9.8.0, SwRI, San Antonio, TX) and executed via the MATLAB-OpenSim API (version 3.3). Six separate probabilistic simulations were generated with a level or 30-degree inclined landing surface; each using a common set of 7 parameter distributions and a combination of muscular co-activation, reflexive gain, and external bushing stiffness (Table 1).

Rotational kinematics of the ankle were kept as defined in the baseline, generic OpenSim model. Briefly the model contained two body-fixed, anatomic axes (Inman, 1976; Delp et al., 1990). Talocrural motion (dorsiflexion (+)/plantarflexion) was defined as rotation of the talus around an oblique, talocrural axis fixed in the tibia (Figure 3). Subtalar motion (supination (+)/pronation) was defined as rotation of the calcaneus around an oblique, subtalar axis fixed in the talus (Inman, 1976). Kinematics were additionally resolved post-hoc to a joint coordinate system (JCS) for the combined ankle complex recommend by the International Society of Biomechanics (ISB) (Wu et al., 2002), which measures non-orthogonal rotations of the body-fixed calcaneus frame relative to the body-fixed tibia frame. The JCS sequence was per convention: dorsiflexion (+)/plantarflexion about a tibia-fixed Z-axis, which was identical to the talocrural anatomic axis; inversion (+)/eversion rotation around a floating axis; and internal (+)/external rotation around a body-fixed calcaneus y-axis (Grood and Suntay, 1983; Wu et al., 2002).

Peak subtalar supination angle was extracted as the primary outcome from probabilistic analysis, along with several complementary biomechanical metrics to aid interpretation and validation: peak talocrural angle, peak angular velocities, and moments borne by each passive bushing at the instant of peak supination. Bushing moments, expressed as vectors in the global reference frame by default, were orthogonally projected (i.e., dot products) onto each anatomic and JCS coordinate axis.

Cumulative distribution functions (CDF) of peak subtalar supination angle were computed from each MC simulation. The 5th, 50th, and 95th percentile values were extracted for comparison (Figures 1C). Sensitivity of the response to each independent parameter was characterized with standard Pearson correlation coefficients, and probabilistic sensitivity factors—which offer a unique assessment of how the mean and variance of the response are influenced by changes in either the mean or variance of each input (Easley et al., 2007).

We also tested the accuracy of the advanced mean value (AMV) probabilistic analysis method, a more computationally-efficient approach relative to traditional MC methods (Wu et al., 1990). Briefly, the AMV method approximates the numerical model (i.e., the sprain simulation) as a first-order Taylor series expansion to estimate the response gradient at specified probability levels. The AMV method requires only tens of trials to generate results, but can be sensitive to non-monotonic system behavior. Thus, we used our MC results as a gold standard to test if AMV predictions were similar at 5, 50, and 95% performance levels using the same inputs. Lastly, we evaluated AMV importance levels, which are the components of a unit vector pointing in the direction of the system response gradient. Thus, importance levels at each performance level have a root-mean-square sum of unity, and provide a more concise sensitivity analysis.

The original purpose of the validated deterministic model was to explore if muscular co-activation or reflexes could theoretically prevent a sprain (DeMers et al., 2017). Through uniform variation of preparatory invertor/evertor co-activation and reflex gain (20 total simulations), the authors of the prior study concluded that co-activation could prevent sprains, while reflexes could not; specifically, 60% co-activation reduced peak subtalar supination from 49 to 30°, while the strongest reflexes (gain of 10) produced an unsafe 41° supination. Our corresponding, probabilistic simulations - Study #3 with strong reflexes (gain 10 ± 1), Study #4 with strong co-activation (60 ± 5%) - aimed to verify if these conclusions held with added uncertainty across a breadth of sprain factors (Table 1). Median response of each MC simulation matched deterministic predictions, although 90th percentile ranges were relatively wide (16°, Study #4, Figure 4). With co-activation only, our sensitivity results showed response variance was most strongly driven by subtalar posture at contact (r = 0.74, Study #4, Table 2), with only weak associations across co-activation, drop height, and muscle strength (r = −0.37, 0.34, −0.21). With reflexes only, 95% of our trials produced 35° or greater peak supination (Study #3, Figure 4), with response variance showing strong associations to drop height and subtalar posture at contact (r = 0.61, 0.55, Study #3, Table 2), but not to reflex gain. These findings corroborate prior conclusions that—even with added intrinsic and extrinsic uncertainty—the strongest reflexes could not prevent sprains for this specific subject and movement (DeMers et al., 2017).

Our analysis of predicted ankle joint mechanics expressed in a JCS highlighted several points. Firstly, 35 ± 6° inversion has been suggested as a range within which lateral ligament damage beings to occur in cadaveric models (with foot 20° plantarflexed, 15° internally rotated), measured about an anterior-posterior axis parallel to the plane of the foot, through the talocrural joint (Aydogan et al., 2006). In this simulation, the generic subtalar joint axis is inclined 38° superior and 9° medial to anatomic directions of the foot (Figure 3). Thus, with 20° talocrural plantarflexion, 52° subtalar supination is required to yield 35° inversion in JCS coordinates, highlighting care needed when comparing across different conventions. As the ISB ankle complex JCS is more commonly used in human motion capture experiments, resolution of default OpenSim outputs to a JCS facilitates comparison. For example, one experiment had subjects drop land onto a 25° inclined platform, and captured 3D kinematics and kinetics of the JCS ankle complex during safe versus incidental sprain-causing trials in two female participants (Li et al., 2018). A range of 35–43° peak inversion in a JCS convention was measured during safe trials, and 55° during sprain trials. Our probabilistic simulations of a 30 cm drop onto a 30° incline predicted 11°—23° peak JCS inversion (Table 3); notably lower than observed injury ranges. Similarly, 600–900°/second inversion velocities were measured during sprain trials, while our simulations predicted a lower 500–700°/second (Li et al., 2018). More severe, virtual sprains could be achieved with a steeper platform incline, or higher drop height, in future applications requiring such.

Selected kinetics were also resolved to the non-orthogonal ankle complex JCS (Table 3). Simulated inversion torque resisted by the rotational bushing representing internal passive anatomy, peaked at 24.5 ± 4.6 Nm in Study #2 in the worst-case of no invertor/evertor muscular contributions or external bracing, and decreased to 4.8–5.3 Nm on average in MC Study’s #4-6 aimed to simulate near-injury scenarios. These latter predictions agree with measurements from the same cadaveric test data (Aydogan et al., 2006); which observed anterior talofibular ligament damage initiation at 2.7–4.3 Nm external inversion torque and 35 ± 6° inversion.

Our addition of a second passive, rotational bushing was aimed to provide a simple means of estimating load sharing between internal anatomy and the net structural effect of a worn brace. Such human-device interactions are difficult or impossible to measure experimentally, but could be useful to inform mechanical design requirements for bracing. Torque-angle data used in the generic bushing definitions can be safely measured on humans with/without worn braces, to provide subject-specific inputs of passive flexibility and external bracing (Siegler et al., 1997; Crabtree and Higginson, 2009). For example, MC Study #5 suggested that - in the absence of invertor/evertor muscular contributions - a parallel, external bushing with a 240% scale factor (Figure 2) would be needed to achieve a similar supination reduction as with 60 ± 5% muscular co-activation alone (Study #4 Table 3). However, the complete absence of muscular contributions is unlikely, and such a rigid brace would likely be uncomfortable during typical ankle motions. Thus, we tested a more realistic scenario including low co-activation, moderate reflex gain, and applied bracing in combination (Study #6, Figure 4 and Table 3). Predictions suggested a more flexible 150% scaled bushing could provide similar protective effect (Figure 4). Future applications may use this framework to estimate design requirements on magnitudes of structural ankle support needed for sprain prevention, to achieve an effective but minimally-restrictive brace design (Dierker et al., 2017). Conversely, structural flexibility of a prototype brace and of a user’s ankle complex can be measured (Siegler et al., 1997), and input to this framework to safely predict probability of a sprain.

An additional aim of this study was to characterize model sensitivities. Standard input-response Pearson correlations generally agreed with probabilistic sensitivity factors computed from MC results (Table 2); moderate and strong correlations corresponded with probabilistic sensitivities that suggested input parameter means drove peak subtalar supination. Increasing variance in correlated parameters also led to increased variance in the response. However, a unique feature of probabilistic sensitivity factors is identification of input-output relationships beyond linear correlation. For example, in MC studies #3 and #4, increased mean ground contact dissipation decreased variance in the response. Similarly, increased plantarflexion at contact increased variance in the response. These relationships are physically sensible and increase confidence in model validity. Specifically, in the absence of external bracing, increased foot-ground damping created a less variable subtalar response, while greater plantarflexion has potential to increase response variance through variance in location of foot-ground contact and moment arms of primary evertors (e.g., peroneus, extensor digitorum). Of note, invertor/evertor muscle strength showed no association with peak subtalar supination across all scenarios. This agrees with clinical studies that have found no predictive relationship between evertor weakness and lateral ankle sprain occurrence (de Noronha et al., 2006; Martin et al., 2021).

Lastly, while probabilistic analysis is a powerful supplement to forward dynamic musculoskeletal simulation, it is often not feasible due to time constraints of traditional techniques. For this reason, we tested the efficient AMV method on Study #6, which included all 10 probabilistic inputs, and found AMV accurately replicated MC predictions at the 5, 50, and 95% levels (Figures 5A). Only 14 simulations were required for the AMV analysis (4 min processing time, laptop Intel i7 CPU), in contrast to the high computational cost of MC (1000 trials, 13 h). Another added benefit of the AMV method relative to MC are probabilistic importance levels, which provide a more succinct, relative ranking of sensitivities at each probability level. For example, importance levels for Study #6 revealed that drop height was the most influential, with increasing effect between 5 and 95% response probability, followed by talocrural and subtalar joint postures, and muscle co-activation, the latter of which showed decreasing influence between 5 and 95% probability levels (Figures 5B).

The probabilistic framework and study findings must be considered with several limitations. The experimental data used to validate the baseline model was limited to one subject (female, 68 kg, 180 cm), and one motion (40 cm unilateral drop landing onto a level surface). Therefore, conclusions should not be generalized to substantially different contexts of use without further validation (Viceconti et al., 2020). Comparison to experimental data of persons with varied anthropometry, landing on inclined surfaces, with/without ankle instability could broaden future utility. An added challenge is the definition of injury thresholds to internal anatomic structures via externally-measured joint kinematics or kinetics. Therefore, we chose not to assume a fixed injury threshold, as the primary focus of this study was on incorporating appropriate input uncertainties and characterizing the most influential sensitivities.