Left Atrial trajectory impairment in Hypertrophic Cardiomyopathy disclosed by Geometric Morphometrics and Parallel Transport

The analysis of full Left Atrium (LA) deformation and whole LA deformational trajectory in time has been poorly investigated and, to the best of our knowledge, seldom discussed in patients with Hypertrophic Cardiomyopathy. Therefore, we considered 22 patients with Hypertrophic Cardiomyopathy (HCM) and 46 healthy subjects, investigated them by three–dimensional Speckle Tracking Echocardiography, and studied the derived landmark clouds via Geometric Morphometrics with Parallel Transport. Trajectory shape and trajectory size were different in Controls versus HCM and their classification powers had high AUC (Area Under the Receiving Operator Characteristic Curve) and accuracy. The two trajectories were much different at the transition between LA conduit and booster pump functions. Full shape and deformation analyses with trajectory analysis enabled a straightforward perception of pathophysiological consequences of HCM condition on LA functioning. It might be worthwhile to apply these techniques to look for novel pathophysiological approaches that may better define atrio–ventricular interaction.

Scientific RepoRts | 6:34906 | DOI: 10.1038/srep34906 of homology: as the anatomical homology is required to compare shapes, physiologically homologous times are needed in order to compare motion trajectories 11,12 . Moreover, these times have to be separately recognized within each individual and variables of interest (univariate such as 3DSTE parameters or multivariate such as shape) are interpolated at these times 11,12 . Following the above path, in this contribution we extend the acquired knowledge of LA functional impairment in HCM by performing a deformation analysis along with a trajectory analysis, i.e. an analysis of the shape of the motion path of LA, that we hypothesize different in HCM patients versus Controls. Figure 1 shows the results of volumetric LA performance. Our results agree with those of the Literature 4-10 . Control and HCM volume curves are dramatically shifted being HCM volumes more than twofold dilated. Control mean volumetric trajectory shows a typical flattening before P-peak while it is absent in HCM. Global, reservoir and conduit ejection fractions are lowered in HCM, as illustrated by box plots. While booster pump ejection fraction seems not impaired, the absolute booster pump is increased in HCM. Controls and HCM differ for several 3DSTE parameters as shown in Supplementary Table 1 presented in Supplementary Online Appendix as already noted 4 , especially in some of the sub-regions automatically computed by Artida device. Our study was done, however, at homologous times and does not compare absolute minima or maxima as it might occur, in different individuals, at different times of their cycles. Figure 2 shows the necessity to adopt a PT procedure when performing shape analysis via classic GPA + PCA. It is evident in fact that, using classic GPA + PCA, PC scores explain concomitantly both inter-and intra-individual variation due to LA contraction. For example, at the positive extreme of PC1 we find a pathological individual whose cycle is limited to PC1 positive values, while other individual trajectories are all limited at the other PC1 extreme. When performing GPA in SSS + LS in SSS + PCA, as shown in Fig. 3, we clearly see what we expect from a cyclic deformation as heart revolution: a series of approximately elliptical trajectories, centered, in this case on the GM of the entire sample. Table 1 reports the variance explained by the first 10 PC scores. Supplementary Table 2 shows that, among the first 10 PC scores, individual significances slightly change during LA revolution. PC1 and PC2 (the dominant dimensions) are significant, again, in correspondence of 3rd and 4th homologous times.

Results
The meaning of first 3 PC scores in terms of traditional 3DSTE parameters can be inferred from Supplementary Table 3 where univariate correlations between the first 3 PC scores and all 3DSTE parameters are reported. PC1 is, as expected, the most correlated PC with the majority of 3DSTE parameters, mainly with longitudinal and circumferential strains; however, PC2 and PC3 present correlations with torsional variables such as twist and rotation that, in at least some local atrial sub-regions, are significant under ANOVA (between   Table 2 shows that in correspondence of ED times (3rd and 4th) Controls and HCM are significantly different under MANOVA. Figure 4 shows the deformed states of Control and HCM LA in correspondence of 5th homologous time, thus at ED. The deformational difference is remarkable: Controls deform much more than HCM patients; this is evident not only on the purely homothetic component (i.e. the mere change in volume/size) but also in the rotational one (as also testified by 3DSTE ANOVAs in Supplementary Table 1) which is evident in the animation presented in Supplementary Video 1.
Classification performance of SVM applied to LS PC scores is shown in Fig. 5. Total accuracy reaches, at homologous time 3th about 97% with sensitivity always larger than specificity. LOOCV is almost constantly (except for homologous time 8th) at about 18%. AUC reaches 0.99 at homologous time 3rd. Figure 6 shows accuracy shows that PCA performed after GPA mixes intraand inter-individual variation, thus impeding to appreciate pure deformations occurring within individuals. Deformations associated to PC axes extremes are shown. These are relative to the GM that is shown in light grey together with coloured shapes corresponding to axes extremes. The colormap ranges from blue (min) to red (max) and refers to |x M − x|, with x M the position of a point in the GM and x its position at the specified axis extreme (positive or negative).  Fig. 2) the appreciation of pure deformation patterns without the noise due to inter-individual variation. Importantly, centering deformations on GM allows also to consider both systole and diastole as deformed states. Deformations relative to the GM (shown in light grey) are depicted. The colormap ranges from blue (min) to red (max) and refers to |x M − x|, with x M the position of a point in the GM and x its position at the specified axis extreme (positive or negative).
performance depicted in Fig. 5 when compared with the same analysis performed on classical 3DSTE global parameters. LS PC scores accuracy performs much better than global 3DSTE variables and than volume itself. The mean trajectories of Controls and HCM patients identified by the first 3 PC scores coming from GPA in SSS + LS in SSS + PCA (whose PC1/PC2 space was shown in Fig. 2) and treated as landmarks are shown in Fig. 7. The difference in both size and shape of the 2 trajectories can be easily seen. In fact, trajectory analysis shows that this difference is statistically significant. PC1/PC2 space of trajectory shape analysis is shown in Fig. 8. The sole PC1 is significant under individual ANOVA (R 2 : 0.21; p-value: 0.001); however, MANOVA using the first 10 PC scores (that explain, collectively about 90% of total variance) remains significant (Multivariate R 2 : 0.07; p-value: 0.002).     shows that PC1 explains a specific trajectory morphology: the bending of the trajectory in correspondence of landmarks 6-9; these correspond to deformations occurring around the P peak that represents a crucial step in LA function (see discussion). PC1/PC2 and PC1/PC3 angles of trajectories orientation were not significant under ANOVA. We correlated trajectory's shape PC1 and trajectory's size with volumetric indicators. Results (Spearman) are shown in Table 3. Classification via SVM was applied to PC1 of trajectory shape analysis (the sole PC significant under ANOVA) and trajectory size (significant under ANOVA). SVM results are shown in Fig. 10. Accuracies of trajectory's shape PC1 and size are 83% and 79% while AUCs are 0.86 and 0.84, respectively. These parameters have, on average, a better classification performance than the majority of volumetric indicators.

Discussion
Shape and deformation analyses here presented were not only useful to characterise and discriminate HCM patients from Controls, but they also had a much higher accuracy and/or diagnostic performance than standard volumetric evaluation. While that LA volume and function are impaired in HCM patients is not a new   Scientific RepoRts | 6:34906 | DOI: 10.1038/srep34906 notion 4-10,13 , shaping the motion via trajectory analysis provided, for the first time, parameters with a high classification power. Apart from classification, the greatest advantage is anyhow represented by the pathophysiological interpretation of LA functional role that our approach might carry. This could also have important translational consequences. The most intriguing evidence of the trajectory shape analysis is the position of the bending occurring before the P peak that is present in Controls and absent in HCM (Figs 7 and 9). This electrical event is of particular importance in healthy LA functioning: at this time we see the shift from the conduit LA function to the booster pump function 14 . This meso-diastolic transition is crucial during the active contractile chamber that augments LV ventricular filling in late diastole and acts a suction source that refills LA in early systole. The difference between Controls and HCM may indicate that the conduit function is impaired not only in terms of extent of LA volumetric performance; indeed, the way the contraction is achieved is expressed by the trajectory shape that, as stated above, is particularly different between the 6th an 9th homologous times measured starting from R-peak. Looking at Figs 7 and 9 enables to see that Control and HCM trajectories shapes begin to diverge at the interval between the 6th and 7th landmarks. Controls bend there downward so that between the 7th and 9th landmarks the trajectory rises. In fact, as shown in Fig. 2, the Conduit function is significantly impaired in HCM and it is the volumetric indicator most correlated with trajectory's shape PC1. On the opposite, trajectory's size, that is smaller in HCM, is correlated mostly with reservoir function, the second volumetric indicator that separates Controls from HCM. It is important to note here that the shape of the trajectory as depicted in our Figs 7 and 9 was reconstructed in the present context by using the first 3 PCs after LS procedure. These explain collectively about 71% of total variance. However, the procedure we describe here can be effectively applied to more than 3 PCs. GPA on trajectory shapes can be performed on hyper-shapes identified by n dimensions, thus by using n PCs. A difficulty rises in this case as the hyper-shape cannot be illustrated and the interpretation of results requires a higher level of mathematical abstraction. This could limit the accessibility of this approach to a wide scientific community. It is for further investigations to look for classification performance of this more advanced approach.
These signals may be detected using LA volume, albeit LS deformation analysis better discriminates the pathology as shown in Fig. 6 and outperforms classic 3DSTE parameters. Indeed, these latter deliver information about LA deformation; however, 3DSTE deformation is evaluated at each time frame (in general, not a homologous time) by comparison with respect to the LA systolic state, a sort of blank that is 0 for all, by definition. As

Volumetric indicators
Trajectory size p-val PC1 trajectory shape p-val  we centered shapes on GM, diastolic and systolic states are both recognized as deformed states and both are put in contrast with a real shape. Our entire procedure makes sense solely upon the identification and interpolation of homologous times 11,12 . Indeed, only temporal homology can guarantee the rationality of comparing deformed states. Interestingly, a cross correlation with LV deformational change might be done in order to explore the morphological and deformational covariance between LA and LV. It was recently shown 15 that in HCM the transmitral flow velocity is much reduced in correspondence of P-peak, just the point where we observed the absence of bending in LA trajectory shape. We hypothesize, therefore, that the functional/morphological covariation between LA and LV express there its maximum potential to be further investigated in a conjoint analysis to derive potentially important translational information. In this way the hypothesis could be tested that a break in the natural deformational covariance between the 2 chambers is present in diseased conditions such as HCM and other diseases [14][15][16] .

Methods
Subjects and Ethic statements. The study was conducted after the approval of the "Dipartimento di Scienze Cardiovascolari, Respiratorie, Nefrologiche, Anestesiologiche e Geriatriche, Sapienza-Università di Roma" after submitting the protocol, as an observational spontaneous investigation, to the pertinent local Ethical Committee. All methods were carried out in accordance with the ethical guidelines of the Declaration of Helsinki. Written informed consent was obtained from each subject and 68 subjects were studied. For 46 healthy subjects we assessed, based on an accurate cardiological visit, the absence of any type of known heart disease. There were 22 non-obstructive HCM patients. Table 4 shows descriptive statistics of the sample under study.
3D data acquisition. The acquisition of LA geometry is achieved by means of 3DSTE, a technique that is today a consolidated standard in cardiological diagnosis, allowing the real time evaluation of LA and LV motion and the measurements of several diagnostic parameters as twist, torsion and longitudinal strain among others. LA geometrical data are collected during the cardiac revolution by means of PST-25SX Artida device, Toshiba Medical Systems Corp., Tokyo, Japan. LA geometries are reconstructed starting from a set of 6 homologous landmarks indentified on the endocardial layer, manually detected by the operator for all subjects under study. The same operator (AE) is involved in LA geometry reconstruction. Additional details and discussion about the comparison between fully automated or manual reconstruction can be found elsewhere 11,12 . The final dataset of any subject is a time-sequence of shapes, each composed by 1297 landmarks (assumed to be homologous) for endocardial surfaces, positioned along 36 horizontal circles, each composed of 36 landmarks, plus the apex. We get the landmark cloud (upon which the standard rotational, torsional and strain parameters are normally computed and outputted by any Artida device) due to an unlocked version of the software equipping our PST-25SX Artida, thanks to a special opportunity provided in the context of an official research and development agreement between the Dipartimento di Scienze Cardiovascolari, Respiratorie, Nefrologiche Anestesiologiche e Geriatriche, Sapienza Università di Roma and Toshiba Medical System Europe, Zoetermeer, The Netherland.
Temporal homology and Interpolation. The first step of our method aims to introduce the concept of homologous times 11,12 . As we used Ecocardiographic data, we had the possibility to use the electrocardiogram associated to it in order to identify homologous times. The same cannot be applied when using MRI and peculiar spatiotemporal registration techniques are needed in order to standardize different individuals at comparable temporal intervals 17,18 . Among others, a spatio-temporal registration method was proposed, based on some specific physiological events 19 . Identifying and defining homologous times is a crucial step for correctly comparing any type of variable among different individuals acquired at different temporal resolutions and thus composed by different number of shapes in any cycle. Temporal homology fully extends in 4D the concept of anatomical homology and allows comparing not only shapes varying in time but also the shape of the trajectories themselves. Figure 11 shows the 12 homologous times chosen to interpolate LA univariate (3DSTE parameters) and multivariate (shapes) descriptors. Specifically, we have 3 strictly homologous times representing specific electromechanical events, i.e. R peak, LA end-diastolic volume, and P-peak; their exact time (in ms) were visually identified Mutation not known -6 Not investigated -2 Table 4. Descriptive statistics ± standard deviation for the populations under study.
on each individual 3DSTE video clip. These 3 times are critical in LA functioning as they separate the 3 main LA functions occurring in the cardiac revolution, as Fig. 11 shows: reservoir (grey), conduit (pink) and active booster pump (green) 14 . The other times are obtained by sampling 3 equally spaced times between successive and strictly homologous times. As shown in Fig. 11, this approach guarantees to have a complete representation of LA trajectory during cardiac revolution. In order to reconstruct shapes, we performed separate Generalized Procrustes Analysis (GPA) in the Size and Shape Space (SSS), i.e. without scaling at unit size, followed by a Principal Component Analysis (PCA); then, we interpolated the full set of PC scores for each individual. GPA centers and optimally rotates shapes, optionally scaling to unit size, in order to remove non-shape attributes 20 . Typically, GPA is followed by PCA performed on aligned coordinates; PCA can be linear or non-linear and allows to get a ranking of the main shape-change modes and to visualize main shape-change directions. We used GPA in SSS rather than in the Shape Space (SS), thus without scaling shapes at unit size, as LA size is an important attribute of LA function. Individuals PC scores and PC loadings got from PCA are used to reconstruct LA shapes and sizes at any homologous time, thus obtaining a set of 12 shapes (each composed of 1297 landmarks) for each of the 68 subjects under study. The collection of these shapes is the object of the successive deformation analysis. The SSS is less used than the SS but its geometrical structure is simpler. Full details can be found in elsewhere refs 20-22 apex. Here we can only summarize its geometrical genesis and some of its important features. We recall that the SS is built starting from the Euclidean space of k × m (k = number of landmarks, m = number of dimensions) configuration matrices via three operations that we can consider (without entering in more subtle mathematical distinctions), as projections: centering via translation, alignment via rotations, scaling to unit centroid size. While the first projection does not change the Euclidean structure of the space, the second and third make the SS a Non Euclidean Space. In particular, the third operation projects each configuration on an hyper-sphere of unit radius. The major consequence of this construction is that, after aligning two configurations, the geodesic line joining it, is an arch of circumference. In order to approximate the curved space around a given point (shape) one can project back each shape to the tangent space to the hyper-sphere on that point. The SSS can be built in the same way, but performing only the first two geometrical operations (translation and alignment) without projecting on the hyper-sphere. It is easy to understand that this procedure, leads, roughly speaking, to a less Non Euclidean space. In particular, in the SSS the geodesics joining two aligned configurations are straight lines: i.e. the same of the Euclidean configuration space. The consequence is that no projection is necessary toward the tangent space because a segment joining two aligned configuration coincides directly with a vector of the tangent space Figure 11. LA volume curve along 2 cardiac cycles. Small dots identify the times at which 3DSTE acquisition is done (observed times). Among these, the 3 strictly homologous times are identified by a triangle, whereas the others homologous times by a red dot. All in all, a fine reconstruction of the actual LA trajectory is obtained. to the SSS on one of the two configurations. Practically this means that, once GPA is performed, the resulting configurations lie on the tangent space at the Grand Mean. This allows us to perform directly a PCA.
Volumetric and deformation analysis. Albeit we do not focus on LA volume impairments, as done by previous investigators 4-10 , it is essential to extract LA volume indicators in order to correlate them with morphometric descriptors that we introduce and describe below. We thus computed the main LA functions using volumes corresponding to the various LA phases as described in Table 4, and showed them in Fig. 1; therein, the difference between absolute booster pump function and the corresponding ejection fraction can be appreciated. A crucial step in a deformation analysis is the management of inter-individual differences, which have to be eliminated; indeed, when different LA geometries moving in time are compared, classic GPA and PCA mix intra-and inter-individual variation, thus inhibiting the appreciation of pure deformation patterns. To this end, we follow a Linear Shift (LS) strategy 11,12,23 . LS consists in the application of the Parallel Transport (PT) tool in the Euclidean tangent plane to the Riemannian space of the shapes, instead of applying it in the Riemannian space itself 22,24 ; it introduces an approximation which was however justified with specific reference to the echocardiographic data 11,12 .
LS strategy allows to estimate, hierarchically, the deformations occurring within each individual and to apply them to a mannequin. Specifically, we identify as our mannequin the Grand Mean (GM) of the entire sample, and use LS to apply individual deformation occurring within each cycle to the GM. It is worth noting that, as 3DSTE strain data are measured for each individual with respect to the end-systolic (ES) frame, typically 3DSTE mannequin is the ES shape and ES state is not a deformed state. Using GM as mannequin, instead of ES, allows to appreciate both ES and end-diastole (ED) as deformed states. Thus the sequence of deformation analysis consists of (GPA in SSS) + (LS in SSS) + PCA. Additional details about PT and LS can be found elsewhere 11,12,22,23 . The outcome of this procedure is a series of deformations of GM (as standard 3DSTE procedure delivers a series of deformations with respect to the ES state). Close points in the resulting PCA plot do not correspond to individuals that have the same shape but indicate individuals that are experiencing the same deformation.
Trajectory analysis. As 3DSTE landmarks identify the shape of LA through a cloud of points, likewise PC scores values of LS analysis represent new landmarks which identify the shape (approximately elliptical) of each individual trajectory 11,12,25,26 . In this sense, we say that trajectory analysis aims at shaping LA motion. Thanks to the notion of homologous times and to the elimination of inter-individual differences via LS, we have a set of deformation trajectories, each composed of an equal number of individual shapes in correspondence of homologous electrophysiological times. We use the first 3 PC scores of LS analysis as landmarks which identify a 3D configuration for each trajectory, and study these latters by exploring their shapes, orientations and sizes. It is worth noting that these new landmarks are still homologous, due to the homologous time interpolation procedure described above.
We study trajectory shapes by performing a second-order GPA in SS + PCA on these shapes. During the GPA only landmarks 1, 5, 9, corresponding to the strictly homologous times identified in Fig. 11, were used to estimate optimal alignment, while the other homologous times are excluded during Procrustes distance minimization process, and are passively appended to GPA transformations (translation, scaling and rotation). This strategy allows to compare relatively complete shapes of trajectories without adding noise due to non perfectly homologous physiologically-based event estimation. PC scores deriving from such analysis are not referred to LA shape but to LA trajectory shape. Trajectory size can be evaluated by performing an ANOVA (between Controls and HCM) on the Centroid Size (CS) of trajectories' configurations. Trajectory orientation is estimated by measuring the direction on PC1/PC2 and PC1/PC3 planes of the vectors identified by the first landmark, i.e. the PC scores values corresponding to R peak and the fifth landmark, i.e. PC scores values corresponding to LA ED. This strategy permits a full evaluation of the LA motion attributes.
Linear models and classification. We used a series of ANOVAs on descriptive 3DSTE parameters in order to assess differences between Controls and HCM at any homologous time. The same was done for the first 10 PC scores of deformation analysis; on these we also performed the multivariate test via MANOVAs.
The first 3 PC scores of deformation analysis were also correlated with per-individual mean centered 3DSTE variables in order to assess their significance in terms of traditional 3DSTE parameters. ANOVAs were also performed on PC1/PC2 and PC1/PC3 angles and on trajectory size, while MANOVA was performed on the first 5 PC scores of trajectory shape (albeit reporting them also for univariate significances). Deformation PC scores (at any homologous time) and trajectories' attributes were also subjected to classification analysis in order to assess the potential of these indicators in recognizing HCM condition. As shown previously 27 , Support Vector Machine (SVM) outperforms classic classification procedures (logistic regression, discriminant analysis and random forest, among others) in discriminating healthy LV from infarcted patients using data very similar to those used here. For this reason, we used SVM with radial basis Gaussian kernel function, hyperparameter C = 1, and with reference to the first 50 LS PC scores evaluated at any homologous time, in order to classify our sample. At any time a Univariate Association Filtering (UAF) with p-value cut-off at 0.05 (thus pretty conservative in a classification problem) was performed as feature selection. The same was done for trajectory attributes. In the trajectory attributes' analysis there are no homologous times as they are not related to any single homologous time, being representative of the entire LA revolution. Any SVM run was evaluated in terms of total accuracy, sensitivity, specificity, Area Under the Curve (AUC) and Leave One Out Cross Validation error (LOOCV).

Conclusions
This study sought to test the hypotheses that LA deformation path is impaired in HCM, and that recognizable abnormalities are present not only in deformation but also in motion trajectory. Indeed, whereas recent studies of LA dysfunction in HCM suggest that LA dilation is augmented while LA strains are reduced [4][5][6][7][8][9][10]13 , no investigations were performed on full LA deformation and whole LA deformational trajectory in time.
We enrolled 22 HCM patients and 46 Healthy subjects. We used 3D Speckle Tracking Echocardiography (3DSTE) derived landmark clouds in order to perform full shape and deformation analyses. This was achieved using Geometric Morphometrics and Parallel Transport that allow to appreciate solely the deformation without the noise due to inter-individual differences 11,12 . We also used the notion of temporal homology in order to interpolate shapes at physiologically comparable times for each individual. Thus not only deformation per se but also the shape of the deformation trajectories during LA revolution were assessed. Both deformation analysis scores and trajectory shape analysis scores underwent classification via SVM in order to test the discriminating power of these parameters in characterizing HCM.
We showed that deformation analysis was highly effective in recognizing HCM: AUC during meso-diastole reached 0.99 and total accuracy 0.81. Trajectory shape and trajectory size were much different in Controls versus HCM and their classification power was high in terms of AUC (0.87 and 0.81 respectively). The two trajectories were particularly different just before the transition between LA conduit and booster pump functions. So, full shape and deformation analyses accompanied by a trajectory analysis allow a straightforward perception of pathophysiological consequences of HCM condition.
It is provocative that by shape analysis traslational differences might be appreciated of a peculiar pathophysiological time frame in LA function that we showed abnormal in HCM, confirming results obtained by different methods 15 , but whose approach might be extended to different conditions such as congestive heart failure whereby LA function might so be disclosed as a potentially important target for individualised therapeutical interventions 16 . It might be worthwhile to apply the complex techniques that we described to look for novel pathophysiological approaches in order to better define atrio-ventricular interaction.