Probabilistic biomechanical finite element simulations: whole-model classical hypothesis testing based on upcrossing geometry

Statistical analyses of biomechanical finite element (FE) simulations are frequently conducted on scalar metrics extracted from anatomically homologous regions, like maximum von Mises stresses from demarcated bone areas. Advantages of this approach are numerical tabulability and statistical simplicity, but disadvantages include region demarcation subjectivity, spatial resolution reduction, and results interpretation complexity when attempting to mentally map tabulated results to original anatomy. This study proposes a method which abandons the two aforementioned advantages to overcome these three limitations. The method is inspired by parametric random field theory (RFT), but instead uses a non-parametric analogue to RFT which permits flexible model-wide statistical analyses through non-parametrically constructed probability densities regarding volumetric upcrossing geometry. We illustrate method fundamentals using basic 1D and 2D models, then use a public model of hip cartilage compression to highlight how the concepts can extend to practical biomechanical modeling. The ultimate whole-volume results are easy to interpret, and for constant model geometry the method is simple to implement. Moreover, our analyses demonstrate that the method can yield biomechanical insights which are difficult to infer from single simulations or tabulated multi-simulation results. Generalizability to non-constant geometry including subject-specific anatomy is discussed. ABSTRACT 11 Statistical analyses of biomechanical ﬁnite element (FE) simulations are frequently conducted on scalar metrics extracted from anatomically homologous regions, like maximum von Mises stresses from demarcated bone areas. Advantages of this approach are numerical tabulability and statistical simplicity, but disadvantages include region demarcation subjectivity, spatial resolution reduction, and results interpretation complexity when attempting to mentally map tabulated results to original anatomy. This study proposes a method which abandons the two aforementioned advantages to overcome these three limitations. The method is inspired by parametric random ﬁeld theory (RFT), but instead uses a non-parametric analogue to RFT which permits ﬂexible model-wide statistical analyses through non-parametrically constructed probability densities regarding volumetric upcrossing geometry. We illustrate method fundamentals using basic 1D and 2D models, then use a public model of hip cartilage compression to highlight how the concepts can extend to practical biomechanical modeling. The ultimate whole-volume results are easy to interpret, and for constant model geometry the method is simple to implement. Moreover, our analyses demonstrate that the method can yield biomechanical insights which are difﬁcult to infer from single simulations or tabulated multi-simulation results. Generalizability to non-constant geometry including subject-speciﬁc anatomy is discussed.


48
Producing a probabilistic input-output mapping is conceptually simple: iteratively change input 49 parameters according to a particular distribution and assemble output parameters for each iteration to 50 yield an output distribution. The simplest method is Monte Carlo simulation which randomly generates 51 input parameters based on given mean and standard deviation values (Dar et al., 2002). More complex 52 methods like Markov Chain Monte Carlo can accelerate probabilistic output distribution convergence 53 (Boyaval, 2012). 54 Once probabilistic inputs / outputs are generated they may be probed using a variety of statistical 55 methods. A common technique is to extract scalars like maximum von Mises stress from anatomically 56 demarcated regions of interest (Radcliffe and Taylor, 2007). Other techniques include Taguchi global 57 model comparisons (Taguchi, 1987;Dar et al., 2002;Lin et al., 2007) to fuzzy set modeling (Babuska and 58 Silva, 2014) and probability density construction for specific model parameters (Easley et   The purpose of this paper is to propose an alternative method which conducts classical hypothesis 61 testing at the whole-model level using continuum upcrossing geometry. An 'upcrossing' is a portion of 62 the continuum that survives a threshold ( Fig.1) like an island above the water's surface or a mountain and 3D continua, respectively. Parametric solutions to upcrossing geometry probabilities exist for n- 66 dimensional Gaussian continua in the random field theory (RFT) literature (Adler and Taylor, 2007), and 67 non-parametric approximations have been shown to be equally effective (Nichols and Holmes, 2002). The 68 method we propose follows the latter, non-parametric permutation approach because it is ideally suited to 69 the iterative simulation which characterizes probabilistic FE analysis. 70 The method is inspired by hypothesis testing approaches in nonlinear modeling (Legay and Viswanatha, 71 2009) and in particular a label-based continuum permutation approach (Nichols and Holmes, 2002). It first 72 assembles a large number of element-or node-based test statistic volumes through iterative simulation, 73 then conducts inference using non-parametrically estimated upcrossing probabilities. These upcrossing Here a is the hyperelastic parameter, k is the elasticity volume modulus, I is the deformation tensor's 100 first deviatoric invariant, and J is the deformation Jacobian. The parameter a was set to 100 and eight k 101 values (800, 817, 834, 851, 869, 886, 903, 920) were compared to a datum case of k=820.    Random Field Theory (Adler and Taylor, 2007). The method is described below and is depicted in Fig.6.

122
All permutations described below were applied to pre-simulated FEA results. The datum Young's modulus (E=14 GPa) was subtracted from the eight 1D Young's modulus continua 125 (Fig.2b), and the resulting difference continua were sign-permuted ( Fig.6a) to generate a number of 126 artificial data samples. For each sample the t continuum was computed according to the typical one-127 sample t statistic definition: where y is the sample mean, µ is the datum, s is the sample standard deviation, N is sample size and q 129 is continuum position. Repeating for all permutation samples produced a distribution of 1D t continua 130 ( Fig.6b), whose maxima formed a 'primary' probability density function (PDF) (Fig.6c). This primary 131 PDF represents the expected maximum difference (from the datum case of E = 14 GPa) that smooth, 132 purely random continua would be expected to produce if there were truly no effect. 133 We conducted classical hypothesis testing at α=0.05 using the primary PDF's 95th percentile (t * )     distributions in each group (Table 1) was tested using a slight modification of the permutation approach 171 described above (Fig.6). The only differences were that (i) the two-sample t statistic was computed instead 172 of the one-sample t statistic, and (ii) group permutations were conducted instead of sign permutations.

173
Group permutations were performed by randomly assigning each of the 20 continuum observations to 174 one of the two groups, with ten observations in each group, then repeating for a total of 10,000 random and that different continuum regions can respond in opposite ways to probabilistic inputs.

192
Although stiffness increased non-uniformly as a Gaussian pulse (Fig.7a) the test statistic magnitude 193 was effectively uniform across that region (elements 60 -80; Fig.7d). This suggests that mechanical and 194 statistical magnitudes are not directly related, and thus that statistical conclusions mustn't be limited to 195 areas of large mechanical signal unless one's hypothesis pertains specifically to those areas. , with absolute t values increasing near element #70 but decreasing elsewhere (Fig.8c), re-emphasizing 201 the complex relation amongst different field variables' response to probabilistic model features.

202
The t distributions for stress and strain were not qualitatively affected by the number of additional 203 FE simulations; 16 simulations, or one extra simulation per observation (Fig.8a,c) yielded essentially the 204 same results as 400 simulations (Fig.8b,d). The reason is that permutation leverages variability in small

Model B 214
The mean stress distributions associated with the three indenter faces (Fig.10) closely followed indenter 215 face geometry (Fig.3). Variation in material parameters was associated with stress distribution variability 216 (Fig.11a). Nevertheless, t values were effectively constant across all elements and all three models 217 (Fig.11b). This suggests that test statistic continua are more robust to model geometry imperfections than 218 are stress/strain continua.

220
A two-sample t test regarding the material parameters (Table 1)

Manuscript to be reviewed
Computer Science mean stresses which were generally higher in Group B vs. Group A (Fig.12), where a stress distribution 223 difference plot clarified that inter-group differences were generally confined to areas of large stress 224 (Fig.13). The inter-group statistical differences were much broader, covering essentially the entire femoral 225 cartilage (Fig.14). Moreover, relatively broad regions of the cartilage exhibited significant stress decreases, 226 similar to the result observed in the simple bone model (Fig.7f).

227
These results reiterate many of the aforementioned methodological points. In particular, changes in    (Fig.7b,e and Fig.7c,f), which eliminates the need to separately tabulate statistical results. effects on stress/strain continua, but have comparably little-to-no effect on test statistic continua 257 (Fig.??), implying that continuum-level hypothesis testing may be more robust than commonly 258 employed procedures which analyze local maxima. This potential danger is highlighted in the 259 more realistic Model C, whose mean differences (Fig.13) exhibited high focal stresses whereas the 260 statistical continuum was much more constant across the contact surface (Fig.14). Manuscript to be reviewed Computer Science

278
In the literature, FE-based classical hypothesis testing is typically conducted via scalar analysis of local 279 extrema (Radcliffe and Taylor, 2007). Applying that approach to the local mechanical change extrema in 280 Model A (Fig.7a-c)  While the test statistic magnitudes are the same for both the proposed whole-model approach (Fig.7) 284 and these local extremum analyses, the critical threshold at α=0.05 is different because the spatial scope 285 is different. The broader the spatial scope of the hypothesis, the higher the threshold must be to avoid 286 false positives (Friston et al., 2007); in other words, random processes operating in a larger volume have a 287 greater chance of reaching an arbitrary threshold.

288
The proposed model-wide approach (Fig.7)  Historically in biomechanical FEA, low sample sizes (frequently n = 1 for each model) permitted  this approach to its perhaps logical extreme -the unknown parameter is allowed to vary randomly within 325 defined limits over a large number of iterations (usually on the order of 10,000). These iterations produce 326 a distribution of results that can be statistically compared with other such distributions.

327
A final advantage is that statistical continua may be less sensitive to geometric mesh peculiarities than 328 stress / strain continua. In Fig.10 and Fig.13, for example, it is clear from the oddly shaped regions of 329 stress difference that these effects were likely caused by mesh irregularities and that remeshing would 330 likely smooth out these areas of highly localized stress changes. The test statistic continuum, on the other 331 hand, appeared to be considerably less sensitive to localization effects (Fig.11) and (Fig.14). This may 332 imply that one needn't necessarily develop an ideal mesh, because statistical analysis may be able to 333 mitigate mesh peculiarity-induced stress distribution irregularities.

344
A second limitation is computational feasibility. Although our results suggest that incorporating a 345 single additional uncertain parameter into the model may not greatly increase computational demand, 346 this may not be true for higher dimensional parameter spaces. In particular, given N experimental 347 measurements, our results show that 2N simulations are sufficient to achieve probabilistic convergence 348 (Fig.9). However, this result may be limited to cases where the uncertainty is sufficiently small so 349 that it fails to produce large qualitative changes in the underlying stress/strain continua. Moreover, the 350 feasibility for higher-dimensional parameter spaces is unclear. In particular, a sample of N observations is 351 likely unsuitable for an N-dimensional parameter space, or even an N/2-dimensional parameter space.

352
The relation between uncertainty magnitude, number of uncertain parameters, the sample size and the 353 minimum number of FE simulations required to achieve probabilistic convergence is an important topic 354 that we leave for future work.

355
A third potential limitation is that both upcrossing features and the test statistic continuum can be 356 arbitrary. In this paper we restricted analyses to the upcrossing maximum and integral due to the robustness 357 of these metrics with respect to other geometric features (Zhang et al., 2009). Other upcrossing metrics 358 and even arbitrary test statistic continua could be submitted to a non-parametric permutation routine.

359
This is partly advantageous because arbitrary smoothing can be applied to the continuum data, and in

369
This paper has proposed a probabilistic finite element simulation method for conducting classical hypoth-