Towards automated ligamentous injury evaluation in 1 syndesmotic ankle lesions

Forced external rotation is hypothesized as the key mechanism of syndesmotic ankle injuries. 37 This complex trauma pattern ruptures the syndesmotic ligaments, inducing a three-dimensional 38 deviation from the normal distal tibiofibular joint configuration. However, current diagnostic 39 imaging modalities are impeded by a two-dimensionalassessment, without taking into account 40 ligamentous stabilizers. Therefore, our aim is two-fold: (1) to construct an articulated statistical 41 shape model of the normal ankle with inclusion of ligamentous morphometry and (2) to apply 42 this model in the assement of a clinical cohort of paient with syndesmotic ankle injuries. Three-dimensional models of the distal tibiofibular years),patients syndesmotic 35 +/- contralateral the statiscal shape model was generated after aligning all ankles based on the distal tibia. The position of the syndesmotic ligaments was predicted based on previously validated iterative shortest path calculation methodology. Evaluation of the model was described by means of accuracy, compactness and generalization. Canonical Correlation Analysis was performed to assess the influence of syndesmotic lesions on the distal tibiofibular joint congruency. Our presented model contained an accuracy of 0.23 +/- 0.028 mm. Mean prediction accuracy of ligament insertions was 0.53 +/- 12 mm. A statistically significant difference in anterior 60 syndesmotic distance was found between ankles with syndesmotic lesions and healthy controls 61 (95% CI [ 0.32 , 3.29], p = 0.017). There was a significant correlation between presence of 62 syndesmotic injury and the morphological distal tibiofibular configuration (r = 0.873, p 63 <0,001). In this study, we constructed a bony and ligamentous statistical model representing the distal tibiofibular joint Furthermore, the presented model was able to detect an elongation injury of the anterior inferior tibiofibular ligament after traumatic syndesmotic lesions in a clinical patient cohort.


78
Ankle sprains are among the most common sport injuries and compromise 5% of all 79 emergency department visits (1, 2). Up to 24% involve injury to the ligaments stabilizing the 80 ankle syndesmosis (3,4). This anatomical complex is comprised by the distal tibiofibular 81 joint (DTFJ) and the syndesmotic ligaments. Its intrinsic stability is provided by the inherent 82 osseous congruence between the convex distal fibula and concave incisura fibularis. 83 The extrinsic stability is is ensured by distinct ligamentous restraints. The anterior inferior 84 tibiofibular ligament is the primary restraint to fibular external rotation, whereas the posterior 85 inferior tibiofibular ligament restricts posterior translation. The interosseous ligament merges 86 with the interosseous membrane proximally and provides resistance to lateral translation of 87 the fibula (5-7). Injuries to these restraints are generally caused by forced dorsiflexion-88 external rotation during foot pronation. As the rotation moment of the talus within the mortise 89 progresses, its wide anterior dome directly drives the distal fibula into supraphysiological 90 external rotation and posterolateral translation, sequentially injuring the anterior inferior 91 tibiofibular ligament , interosseous ligament and finally the posterior inferior tibiofibular 92 ligament (8,9). These lesions, especially when subtle, present a deceitful diagnostic 93 challenge. Specific tests during clinical examination lack sufficient sensitivity, whereas plain 94 radiographs are highly affected by foot rotation and do not correlate with presence of 95 syndesmotic injury (10-12). Due to these flaws of traditional diagnostic modalities, 96 syndesmotic lesions have a high propensity for misdiagnosis and subsequent mistreatment, 97 potentially leading to long-term sequelae including chronic syndesmotic instability and post-98 traumatic osteoarthritis (8,9). Conventional Computed Tomography (CT) imaging provides 99 three-dimensional (3D) assessment of bony syndesmotic geometry, which allows for precise 100 evaluation of DTFJ configuration. Magnetic Resonance Imaging (MRI) boasts both 101 outstanding sensitivity and specificity in diagnosing specific ligamentous damage (13,14).
However, these traditional imaging modalities do not represent for weight-bearing conditions, 103 potentially underestimating lesion severity. The recent advent of weightbearing cone-beam 104 CT (WBCT) overcomes these drawbacks by imaging both ankles while in bipedal stance. 105 Using WBCT, novel noninvasive 3D imaging methods have already shown to identify 106 patients with a history of syndesmotic injuries (15,16). However, the manner of assessment 107 has been highly variable. Various methods have attempted to objectify DTFJ configuration, 108 mostly by use of manual measurements on 2D axial slices (17). These 2D metrics, however, 109 correlate poorly with the actual 3D deviation of the fibula. (18) . Indicated by the drawbacks 110 of previously described imaging modalities, an accurate diagnosis of syndesmotic lesions 111 should involve 3D weightbearing osseous imaging in combination withpreferably 112 automated-inclusion of patient-specific ligamentous information. This gap can be bridged by 113 creating idealized generic musculoskeletal models, which allows for scaling, morphing and 114 fitting into patient-specific ankle anatomy. 115 Acclaimed for their ability to reduce data dimensionality, statistical shape models (SSM) are 116 well-suited for 3D shape correlation to predict patient-specific missing anatomical 117 information (19). In this context, SSM has already been efficiently utilized for computer-118 assisted Anterior Cruciate Ligament (ACL) reconstruction (20), virtual reconstruction of 119 acetabular bone defects (21) and planning of corrective osteotomies for malunited forearm 120 bones (22). 121 In regard of ligament modelling, earliest techniques represent deformable soft-tissue 122 structures as straight segments connecting commensurate points (23). This has been shown to 123 yield inaccurate results, which can be explained by the fact that soft-tissues follow smooth 124 curved paths due to bony prominences. Recent techniques (24-26) overcame these limitations 125 by describing a novel methodology based on discrete element rigid body spring models to 126 describe muscle and ligament wrapping paths. 127 Combining geometric statistical models and soft-tissue wrapping methodology has not yet 129 been described in the ankle joint. Nevertheless, this could open the door for ligamentous 130 injury prediction in ankle sprains. Therefore, we aimto build an articulated, multi-object 131 statistical shape model representing ligamentous anatomy of the DTFJ. In addition, we 132 hypothesize that by applying this model to the clinical entity of traumatic syndesmotic 133 lesions, automatic ligament prediction can identify patterns in ligamentous injury. predicting ligament insertions on new, unseen cases is presented in Figure 1d. 159 160

Ligament Injury Evaluation 161
Mean length of the Anterior inferior tibiofibular ligament was 12.26 ± 1.89mm for the control 162 cases, 12.32 ± 1.58 mm for the contralateral cases and 14.13 ± 1.48 mm for the cases with 163 syndesmotic lesions. Statistically significant differences were found between the latter two 164 (95% CI [ 0.323.29], p = 0.017) (Fig. 2a). Mean length of the Posterior inferior tibiofibular 165 ligament was 13.44 ± 0.84 mm for the control cases, 16.03 ± 0.90 mm for the contralateral 166 cases and 15.99 ± 0.75 mm for the cases with syndesmotic lesions. No statistically significant 167 differences were found between the latter two (95% CI [-1.42, 1.35], p = 0.96) (Fig. 2b). 168 Mean length of the IOD was 3.78 ± 0.89 mm for the control cases, 4.97 ± 0.94 mm for the 169 contralateral cases and 5.18 ± 0.92 mm for the cases with syndesmotic lesions. No statistically 170 significant differences were found between the latter two (95% CI [-0.71 , 1.12], p = 0.65) 171 ( Fig. 2c). Mean lengths of each fiber as well as p-values are listed in Table 1. 172 173

DTFJ configuration in syndesmotic lesions 174
As defined by use of canonical correlation analysis (CCA), there was a significant correlation 176 between sustaining a high ankle sprain and the DTFJ configuration (r = 0.873, p <0,001). 177 Figure 3 provides a visual representation of the reduced data, mapped into 2D space by use t-178 SNE. It demonstrates that t-SNE was able to cluster WBCT syndesmotic ankles as a separate 179 entity, apart from their contralateral counterpart and non-weightbearing controls. 180 Additionally, t-SNE was able to cluster the contralateral (weightbearing) controls on the 181 periphery of the control group. This is presumably due to a difference in posterior inferior 182 tibiofibular ligament distance between the weightbearing and non-weightbearing controls, as 183 focused mainly on anatomical variance of the separate tibia and fibula rather than the 203 morphological variance of their articular congruence. Regarding compactness, these previous 204 models show a high amount of variance in the first components. This did not apply to our 205 study since the bones were scaled, which eliminated size differences. Also, only the distal part 206 of the tibia and fibula were used, further reducing the influence of size. Prominent shape 207 variance in these previous models were attributed to differences in shape of the proximal or 208 middle third of the tibia and fibula. Consequently, these shape variations could not contribute 209 to the variance of our model. Evidently, we suspect that both the scaling and cutting 210 procedure are the cause of the different compactness curve of our models compared to 211 previous shape models including the tibia/fibula. distance. This represents an elongation/rupture of the Anterior inferior tibiofibular ligament, 220 resulting in an 'anterior open-book injury', which could also be confirmed during arthroscopy 221 (Fig. 4). Additionally, these findings are supported by previous cadaveric studies. Xenos et al. 222 (33) found the fibula to rotate externally 2.7 degrees on average after sectioningthe Anterior 223 inferior tibiofibular ligament. The observation of anterior diastasis in patients with 224 syndesmotic instability are further reinforced by the investigations of Hagemeijer et al. (34) Our study additionally showed a tendency of increased anterior inferior tibiofibular 226 ligamentlength, posterior inferior tibiofibular ligament length and interosseous diastasis of 227 healthy ankles after weight-bearing. These findings are similar to these of Malhotra et al (35), 228 describing external rotation, lateral and posterior translation of the loaded fibula. 229 Further analyzing our data, t-SNE was performed to investigate whether this algorithm is able 230 to distinguish syndesmotic DTFJ anatomy from the control anatomy. As demonstrated, the 2D 231 scatter plot depicts a distinct clustering of the syndesmotic lesion group, compared to the 232 controls, indicating an underlying anatomical dissimilarity. 233 We acknowledge that our study has several important limitations. Firstly, control patients 234 were collected from a database of angio-CT images, ordered for clinical suspicion of vascular 235 disease and presumably do not represent the general population. Secondly, our clinical cohort 236 consisted of a small amount of patients with syndesmotic lesions included, which might 237 impair the generalizability of the findings. However, all patients were carefully selected after 238 a detailed arthroscopic examination, which is currently considered as the gold standard in 239 detecting syndesmotic ankle injuries (36, 37). With regard to confirm the findings of our 240 study, future research should therefore focus on the inclusion of more patients with high ankle 241 sprains, as well as inclusion of healthy asymptomatic control patients from similar age 242 groups. Thirdly, ligamentous insertion sites were based on description in anatomical atlases, 243 and our ligament modelling algorithm was not validated to the true subject-specific anatomy. 244 However, a previous anatomical study by Williams et al. (38) has shown limited variance in 245 ligament insertions between subjects. Finally, t-SNE clusters represent non-linear reductions 246 in dimensionality, which remain difficult to interpret (39). Nevertheless, rapid evolving 247 technology in this field has high potential for future algorithms interpreting these low-248 dimensional data clusters. (40, 41) 249 Thepresent tibiofibular shape model has a high propensity for application in various settings. 250 Our ligament modelling algorithm could be used for computer-aided ligament reconstruction 251 anatomy, e.g. anatomical lateral ankle ligament reconstruction or ACL reconstruction. Since 252 contralateral ankle anatomy is not always available to template asymptomatic anatomy, our homogenously arranged ankle meshes was performed to account for undesired scaling and 324 rotational and translation information. Since the specific hypothesis to analyze changes in 325 DTFJ configuration, all ankles were aligned based on the tibia. The Procrustes transformation 326 matrix of each tibia was subsequently used on the corresponding fibula, in order to retain the 327 subject-specific anatomical configuration of the DTFJ. In doing so, variance in fibular 328 positioning relative to the tibia could be visualized. Although the dataset consisted of mixed 329 right and a left ankles, Procrustes analysis resulted in a perfectly 'mirrored' configuration. 330 Model evaluation was then performed following Principal Component Analysis (PCA). 331 Theoretically, PCA decomposes a multivariate dataset into its mean and corresponding 332 covariance matrix. The eigenvectors of the covariance matrix are usually referred to as 333 principal components or eigenmodes, whereas the eigenvalues indicate there relative 334 importance. An overview of the SSM construction is depicted in Figure 5. and each arbitrarily divided into 11 equidistant nodes (Fig 6 A-B). These elastic lines are 346 progressively released to a position of lesser potential energy without penetration of the 347 surrounding bony contours resulting in a piecewise linear approximation to the true path. By 348 doing so, a computational efficient solution can be obtained. In case the algorithm results in node penetration of any local obstacle, the surface point on the obstacle closest to the local 350 minimum is withheld. Penetrating nodes are consequently penalized to return to the closest 351 point on the surface (Fig. 6 E-F). After calculation of each elastic line, nodes are 352 interconnected forming a flat mesh representing the ligament. For transition to a volumetric 353 model, each ligament was assigned a ligament-specific thickness, based on the literature (7). 354 interosseous ligament distance was modelled differently. Five insertion points were defined 355 on the anatomical location within the incisura. Each point was subsequently connected with 356 its geometrical closest neighbour on the fibula (Fig 6C). Ligaments could then automatically 357 be modeled on newly segmented cases by non-rigidly morphometric fitting its bony structures 358 onto their closest neighbor in the statistical shape model (Fig 7). Model accuracy is defined as an average error between the original bone samples and their 364 model reconstruction of the bone ′ using the first principal components that describe 95% of 365 variance within the dataset. This is calculated as the Root Mean Squared Error (RMSE) of the 366 distances between the original and the reconstructed shape, averaged over all samples. 367 To evaluate the accuracy of ligament modelling, RMSE was separately calculated for each of 368 the ligaments insertion sites. 369 370

Compactness 371
A shape model should explain as much variance in the dataset as possible while keeping the 372 model relatively compact. Evaluation of the model's compactness was done according to the 373 amount of components that account for a 95% of the accumulated variance. Considering the variable impact of size of bones on the variance in the ankle, as described by Audenaert et al. 375 (27), size dominance was reported as the percent variance described by the first principal 376 component. 377 378

Generalization 379
The model generalization is the capacity of the model to accurately describe samples outside 380 of the dataset. This can be calculated by the means of leave-one-out analysis. Shape models 381 were built using an increasing numbers of samples (K = 1, 5, 10, … , 95, 99). For each 382 increment, the corresponding shape model was used to reconstruct the excluded samples. 383 Subsequently, for each increment, the RMSE of the difference between the original shape 384 and the reconstructions ′ was calculated and averaged. To assess the algorithm's capability 385 to consistently model ligament insertions on new samples, Leave-one-out analysis was 386 performed by iteratively excluding one sample in the shape model. Of the excluded model, 387 the nearest neighbor in the shape model was non-rigidly registered, after which shape-specific 388 ligaments were modelled. Hereafter, RMSE was calculated between corresponding insertion 389 points of the actual ligament and its modeled equivalent. RMSE was subsequently averaged 390 over all iterations. 391 392

Ligament injury prediction 393
Ligaments were modelled for each of the 76 controls, 13 contralateral controls and 13 ankles 394 with syndesmotic injury. Subsequently, the length of each ligament fiber was calculated. 395 Difference in length of healthy contralateral versus syndesmotic injured ligaments were 396 statistically analyzed by use of two-tailed Two-Sample student's t-test. 397 398

DTFJ configuration in syndesmotic lesions 399
Since initial alignment was based solely on tibia anatomy and the corresponding fibula was 401 transformed along the same transformation matrix, positional changes in the fibula between 402 all meshes could be visualized. Aligned data were classified as control (+1) or syndesmotic 403 lesion (-1). Statistical analysis of the relationship between syndesmotic injury and DTFJ 404 anatomy was performed by means of canonical correlation analysis in the Euclidian subspace. 405 CCA serves as a technique to find commonalities between two sets of multivariate data (46).

Authors' contributions
All Authors contributed equally to this article. MP designed, coordinated and drafted the 421 manuscript. TH aided with the data collection and processing. AB co-designed the research 422 hypothesis, helped with data collection and co-drafted the manuscript. SDM supervised 423 drafting of the manuscript. MVW helped drafting the manusript KB collected patient data, co-424 designed the research hypothesis and aided drafting the manuscript. JV supervised the final 425 manuscript. EA supervised the research hypothesis, data processing, result analysis and 426 drafting of the manuscript. All authors read and approved the final manuscript. 427 428