Hakim Benkirane
email: hakim.benkirane@centralesupelec.fr

Yoann Pradat
email: yoann.pradat@centralesupelec.fr

Stefan Michiels
email: stefan.michiels@gustaveroussy.fr

Paul-Henry Cournède
email: paul-henry.cournede@centralesupelec.fr

CUSTOMICS: A VERSATILE DEEP-LEARNING BASED STRATEGY FOR MULTI-OMICS INTEGRATION A PREPRINT

Keywords: Multi-omics, Variational Autoencoders, Representation Learning

Recent advances in high-throughput sequencing technologies have enabled the extraction of multiple features that depict patient samples at diverse and complementary molecular levels. The generation of such data has led to new challenges in computational biology regarding the integration of highdimensional and heterogeneous datasets that capture the interrelationships between multiple genes and their functions. Thanks to their versatility and ability to learn synthetic latent representations of complex data, deep learning methods offer promising perspectives for integrating multi-omics data. These methods have led to the conception of many original architectures that are primarily based on autoencoder models. However, due to the difficulty of the task, the integration strategy is fundamental to take full advantage of the sources' particularities without losing the global trends. This paper presents a novel strategy to build a customizable autoencoder model that adapts to the dataset used in the case of high-dimensional multi-source integration. We will assess the impact of integration strategies on the latent representation and combine the best strategies to propose a new method, CustOmics (https://github.com/HakimBenkirane/CustOmics). We focus here on the integration of data from multiple omics sources and demonstrate the performance of the proposed method on test cases for several tasks such as classification and survival analysis.

Introduction

With the advent of high-throughput technologies, multiple omics data are more and more available to characterize the molecular portraits of patients, notably in oncology. Analyzing this mix of data sources and leveraging their information to improve our understanding of biological phenomena remain challenging.

Several studies and projects have made available cohorts of data that are characterized by several molecular sources. The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) has profiled thousands of tumor samples for multiple molecular assays and made available several types of data such as genome sequencing, RNA sequencing, DNA methylation, proteomics, etc. Integrating such diversified data is crucial to reduce the uncertainty over a biological process due to various experimental setups and to unravel interactions that a single source cannot extract. This high diversity of data comes with two significant challenges. The first one is linked to the high-dimensionality of the data. Due to the genetic complexity of the human molecular profile, omics data generally suffers from the 'curse of dimensionality' resulting from the high number of features compared to the smaller number of samples [START_REF] Misra | Integrated omics: Tools, advances, and future approaches[END_REF]. This high-dimensional space often contains correlated features that result in high redundancy, reducing the prediction performance of algorithms [START_REF] James | An Introduction to Statistical Learning: with Applications in R[END_REF]. The second issue lies in data heterogeneity: having its origins from different sources and expressing different phenomena in humans' biological system, omics data is very diverse [START_REF] Bersanelli | Methods for the integration of multi-omics data: Mathematical aspects[END_REF]. For example, transcriptomics and proteomics are normalized differently from other omics data and use different scaling before their analysis, leading to different ranges and data distributions. Finally, omics data like metabolomics can also generate sparsity, as some variables can be below the detection limit and hence be assigned null values.

To overcome those challenges and alleviate overfitting in prediction tasks, standard methods either manually select a small subset of molecular features based on domain knowledge [START_REF] William S Bush | Biofilter: a knowledge-integration system for the multi-locus analysis of genome-wide association studies[END_REF] or use dimensionality reduction techniques before downstream analysis. However, due to the high heterogeneity of the sources, those methods can overlook some genome-wide hidden patterns that justify the need for more flexible methods.

In the past few years, multi-omics integration has been a very active research subject in health science and precision medicine [START_REF] Zitnik | Machine learning for integrating data in biology and medicine: Principles, practice, and opportunities[END_REF]. Multiple statistical learning methods have been proposed from the desire to investigate various complex molecular systems behind cancer. One of the most well-known methods engraved in the statistical field is the Principal Component Analysis (PCA) (extensively presented, for example, in [START_REF] Jolliffe | most relevant genes responsible for the discrimination between each subtype against the others using CustOmics for both integration phases[END_REF]). From this standard model, several variations have been explored, which include Multiple Factor Analysis [Sánchez et al.], Consensus PCA, and multi-block PCA [START_REF] Hassani | Deflation strategies for multi-block principal component analysis revisited[END_REF]. All those methods try to extend the Principal Component approach to a multi-source framework by considering each source as a distinct block. Another variation of the PCA is the Non-negative Matrix Factorization (NMF) first introduced in [START_REF] Zhang | A novel computational framework for simultaneous integration of multiple types of genomic data to identify microrna-gene regulatory modules[END_REF] which uses the same philosophy as PCA but considers a non-negative constraint instead of an orthogonality one. Additionally, we can consider the joint Dimensionality Reduction methods (jDR) as an extension of factorial methods studied in a comprehensive benchmark in [START_REF] Cantini | Benchmarking joint multi-omics dimensionality reduction approaches for the study of cancer[END_REF]. Such methods can be divided into several categories including Bayesian methods [START_REF] Lock | Bayesian consensus clustering[END_REF]] that use assumptions on data distribution and dependencies to build a statistical model, as well as network-based methods that rely on a network representation and can identify modules of the disease-associated mechanisms. The genes are represented by nodes, and edges represent the association between those genes [START_REF] Vaske | Inference of patient-specific pathway activities from multi-dimensional cancer genomics data using paradigm[END_REF][START_REF] Bonnet | Integrative multi-omics module network inference with lemon-tree[END_REF]. Other more specific methods have been explored, mainly for clustering. For example [START_REF] Newman | Autosome: A clustering method for identifying gene expression modules without prior knowledge of cluster number[END_REF] that presented AutoSOME which is a framework using the Self-organizing Map approach and couples it with a density equalizer; this equalizer reduces dimension and preserves the local topology of the gene expressions. Another well-known multi-omics clustering method is iCluster [START_REF] Shen | Integrative subtype discovery in glioblastoma using icluster[END_REF], which uses the same idea as NMF but without the non-negative constraint and allows more diversity in data distributions.

With the rise of deep learning as the new state-of-the-art approach in many medical applications, several studies have explored the interest of these methods for multi-omics data integration. The primary use of these methods has been with autoencoders [START_REF] Rumelhart | Learning Internal Representations by Error Propagation[END_REF], notably by [START_REF] Chaudhary | Deep learning-based multi-omics integration robustly predicts survival in liver cancer[END_REF] in the context of survival prediction. The framework OmiVAE [START_REF] Zhang | Integrated multi-omics analysis using variational autoencoders: Application to pan-cancer classification[END_REF] introduced new improvements for standard representation learning and OmiEmbed [START_REF] Zhang | Omiembed: A unified multi-task deep learning framework for multi-omics data[END_REF] for multitask learning. [START_REF] Simidjievski | Variational autoencoders for cancer data integration: Design principles and computational practice[END_REF]] also proposed multiple deep learning architectures based on variational frameworks for multiple tasks. Huang et al. [START_REF] Huang | Salmon: Survival analysis learning with multi-omics neural networks on breast cancer[END_REF] developed a neural network architecture to integrate multi-omics data for survival by feeding the eigengene matrix instead of the raw data. Based on these methods, several application studies on specific cancer cases were presented, like [START_REF] Hira | Integrated multi-omics analysis of ovarian cancer using variational autoencoders[END_REF] that used the OmiVAE framework to conduct an in-depth study of the Ovarian Cancer's mechanisms. [START_REF] Tong | Deep learning based feature-level integration of multi-omics data for breast cancer patients survival analysis[END_REF] also used standard autoencoders to explore multi-view learning in survival analysis for Breast Cancer. To our knowledge, no benchmarking study has explored and compared the different deep learning approaches and strategies for multi-omics data integration in the context of multitask learning.

In this work, we discuss different strategies for the integration of high-dimensional multi-source data to learn lowdimensional latent representation from multi-omics datasets. Moreover, we introduce a new architecture called CustOmics that combines the advantages of the different strategies and alleviates some of their limitations. We evaluate the impact of this new method on different test cases: we first apply it to a pan-cancer dataset and then study how it handles smaller datasets by evaluating it on specific cohorts. We also provide a new package to help bioinormaticians and computational biologists build deep learning architectures that allow an easy switch between the multiple strategies to adapt to specific use cases.

Materials and methods

Representation Learning for Multi-Omics Integration

Representation Learning is a field of statistical learning that aims to automatically discover relevant representations of the input data [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF]. In the field of multiple omics data integration, it can help synthesize the heterogeneous distributions into a common space, revealing the interaction between the different sources. From a biological point of view, each omics source represents a specific view of the processes involved. Therefore, the main goal behind integration is to uncover a broader view of the phenomenon of interest and to better understand the underlying processes.

Each patient is characterized by K omics vectors (x k ) 1≤kK ∈ R M k with M k the number of features of the k-th source. The models presented in this paper aim at mapping the set of omics vectors to a vector z ∈ R m with m K k=1 M k , the latent representation. The latent features are the components of the vector z. We will use factorial methods as a basis of comparison as they have proved to be the reference method for mapping input data into a latent space of lower dimension, notably for multi-omics integration, as illustrated by the recent review of [START_REF] Cantini | Benchmarking joint multi-omics dimensionality reduction approaches for the study of cancer[END_REF].

For factorial analysis, we consider approaches based on the study of principal components using linear operations through Multiple Factor Analysis (MFA). Deep learning methods use autoencoder architectures to build the latent representation by jointly training encoder and decoder functions. The types of architecture can be divided into three main categories:

• Early Integration (EI): Refers to the set of methods that aim at merging the different sources before the dimensionality reduction. Instead of giving as model inputs a set of separate vectors for each source, we first concatenate the vectors of all the sources. It is characterized by its simplicity, but finds its limitations when some sources bear more significant signals than the others, which makes it difficult for the model to learn interactions between sources.

• Joint Integration (JI): In this setting, sub-representations are created inside the same model for each source before learning the output. This is the most promising strategy theoretically and is the most popular in the literature. However, it is challenged by the sources' heterogeneity, as they do not necessarily follow the same learning dynamics and may need different approaches with specific losses.

• Late Integration (LI): Consists in learning an output for each source separately, with its own representation model. It can adapt well to the specificities of each source, but does not retrieve any cross-modality interaction.

To alleviate the issues brought by those different integration methods, we propose a hybrid strategy, named mixedintegration, that comes halfway between the joint and late integration. It consists of two learning phases: the first phase independently trains the network of each source with an adapted loss to create sub-representations. Those specific models will then be jointly trained in a second phase to build a global output representation. That way, we can adapt our training to each source that may have different learning dynamics in the first phase, while still learing cross-modality interactions in the second phase.

Variational Autoencoders

A Variational Autoencoder (VAE) is a deep generative model which can learn meaningful data representations from high-dimensional input data. It is an extension of standard autoencoders in which the encoder encodes the input as a distribution over a latent space instead of a single point.

In this case, the encoding function q represents a variational distribution (known as encoding distribution) q φ (z|x) in which φ is the parameter to estimate and the decoding function p represents the posterior p θ (x|z).

The particularity of a VAE is its ability to encode a distribution. After the encoding phase, there is a sampling phase in which we sample points from the distribution q φ (z|x).

Traditionally, the distributions in the VAE architecture are supposed Gaussian: the encoder function will learn the two parameter vectors µ and σ (respectively the means and standard deviations) that are used to generate samples in q φ (z|x) using the reparametrization trick z = µ + σ where ∼ N (0, I) and denotes the element-wise product.

The loss function for this architecture can be written as the sum of two distinct losses. First, a reconstruction loss that focuses on the autoencoder's ability to reconstruct the data: Representation of the mixed integration approach. In phase 1, a specific model is trained for each source independently and embeds a sub-representation adapted to the source's specificities. In phase 2, the specific models are trained jointly to create the final output.

L recon = E q φ (z|x) [p θ (x|z)] (1)
It can be interpreted as the conditional entropy of x over z which quantifies the amount of uncertainty that one has over the joint distribution (x, z), knowing z. In a more practical way, it is related to the quantity of information of x that is retained by z (thus qualifying a reconstruction potential).

Second, a regularization loss aims at getting the encoding distribution as close as possible to the theoretical distribution of the latent vector. Traditionally, a Kullback-Leibler divergence is used:

L reg = D KL (q φ (z|x)||p θ (z)).
The total loss will benefit from the framework introduced in Higgins et al. [2017], L = L recon + βL reg , adding a weight on the regularization term to balance between the two parts of the loss function, with a hyperparameter β to optimize.

Early Integration

An early integration autoencoder is an autoencoder applied to the concatenated data from various sources. This type of architecture was reported in [START_REF] Viaud | Representation learning for the clustering of multi-omics data[END_REF] as the SDAE (Standard Deep Autoencoder). The drawback of such an integration is that the resulting concatenation of every dataset results in a noisier and more complex high-dimensional matrix that makes the training more difficult. Moreover, this type of integration does not consider the specific data distribution of each source separately and thus leads the model to learn irrelevant patterns that can, in the worst case, mask the essential signals. For those reasons, even though this strategy has been introduced in the benchmark, it is not recommended in a real setting for complex multi-source datasets such as omics data. For comparison purposes, we adapt the SDAE architecture to a Variational Autoencoder instead of a standard autoencoder as described in Figure 2 Joint Integration This is a standard way to approach multi-source integration using Variational Autoencoders. This model has been used in multiple designs in literature [START_REF] Zhang | Integrated multi-omics analysis using variational autoencoders: Application to pan-cancer classification[END_REF]; [START_REF] Zhang | Omiembed: A unified multi-task deep learning framework for multi-omics data[END_REF]; [START_REF] Simidjievski | Variational autoencoders for cancer data integration: Design principles and computational practice[END_REF]) and is probably the most popular as it is relatively straightforward. The idea is to reduce through a series of the hidden layer the dimension of each source separately, then concatenate them in a central layer before sampling. In this architecture, we would learn the mean and standard deviation of the concatenated data such that:

• z = N (µ, σ 2 I)

• µ = q µ (T (q (1) (x 1 ; θ e,1 ), ..., q (K) (x (K) ; θ e,K )))

• σ = q σ (T (q (1) (x (1) ; θ e,1 ), ..., q (K) (x (K) ; θ e,K )))

Late Integration

A late integration approach for variational autoencoders is a variation of the architecture proposed in Simidjievski et al.

[2019] with the Hierarchical VAE architecture. The idea is to use variational inference on each layer separately and input the resulting representations in an aggregation network learned independently of the other VAEs, which is the main difference from the joint integration approach. This architecture is described in Figure 2.

Mixed Integration / CustOmics

This integration strategy will be the foundation for building the CustOmics architecture. The proposed method is a hierarchical mixed-integration that consists of an autoencoder for each source that creates a sub-representation that will then be fed to a central variational autoencoder. This new integration strategy benefits from two training phases. The first phase will act as a normalization process: each source will train separately to learn a more compact representation that synthesizes its information with less noise. This will help the integration as we will lose all imbalance issues between the sources and avoid loss of focus when a source has an inferior dimensionality or weaker signal than the others. The second phase will constitute a simple joint integration between the learned sub-representations, while still training all the encoders to fine-tune those representations as some signals are enhanced in the presence of other sources.

Regarding the regularization loss for the central layer, the KL divergence can be an obstacle to generalization. As stated in [START_REF] Zhao | Infovae: Information maximizing variational autoencoders[END_REF], the KL divergence suffers from various problems. The first is that the model can fail to learn a meaningful representation of the input. Indeed, the KL divergence can sometimes be too restrictive and naturally tends to make the latent code a random sample from p θ (Z). The second is that the KL divergence can make the model overfit and learn a latent code with variance tending to infinity.

So instead, we will use the Maximum Mean Discrepancy (MMD) to assess the distance between the distributions. This distance stands on the foundation that two distributions are identical if and only if their moments are the same. Let p, q be two distributions, the MMD distance is given as follows:

M M D(p(x)||q(x)) = E p(x),p(x ) (κ(x, x )) + E q(x),q(x ) (κ(x, x )) -2E p(x),q(x ) (κ(x, x )) (2)
where κ is a kernel function. We will choose a Gaussian kernel κ(x, x ) = e

-||x-x || 2 2σ 2
.

Whereas all the models mentioned previously build the latent representation in an unsupervised way, we also create latent features adapted to specific tasks like classification or survival. This idea has been used multiple times in the literature, for example in [START_REF] Zhang | Omiembed: A unified multi-task deep learning framework for multi-omics data[END_REF].

The solution relies on adding a task-related loss to the autoencoder objective function. Therefore, we denote by L task the loss such that L tot = L AE + αL task .

For the classification task, we use a categorical cross-entropy loss defined by L class = i y i log(ŷ i ), where y i is the ground truth for the i th sample, and ŷi its estimation with the downstream model.

For the survival task, we use the DeepSurv loss function. This nonlinear proportional hazard model is introduced in [START_REF] Katzman | Deepsurv: Personalized treatment recommender system using a cox proportional hazards deep neural network[END_REF]. The model is built by using the negative partial log-likelihood formula that translates in our case into:

L(θ) = - i:Ei=1 (μ(t, x i ; θ) -log j∈R(Ti) e μ(xi;θ) ) (3)
where μ(x; θ) is the risk function estimated by the output layer of the network, R(t) is the risk set, that is the set of patients still at risk of failure after time t.

Additionally, this loss will be assigned to each omic-specific network in the first training phase to create adequate sub-representations before the joint integration phase.

Interpretability

To build a more interpretable architecture, we adapt the method introduced in Withnell et al. [2021] to compute SHAP values [Lundberg and Lee, 2017] for deep variational autoencoders. Whereas this method has only been conceived for single source inputs, we expand it to the multi-source setting of CustOmics by adapting it to any deep autoencoder and applying it to each source autoencoder thanks to the dual-phase approach characterizing the mixed-integration strategy. After training the CustOmics network, we pass each intermediate autoencoder to a DeepSHAP explainer that was modified to make sure that it could take either each component of the latent vector or the prediction output for the contribution analysis. For each prediction, SHAP values corresponding to each gene or latent dimension are computed to estimate the overall contribution by averaging them over a group of samples with similar features. When combined with the mixed-integration approach through CustOmics, we can compute for each source two types of features importance.

During the first training phase, we compute the features importances of each source separately to assess the importance of genes when only considering this source. During the second training phase, the features' importance can change as it will then consider other sources, showing how some genes can have increasing or decreasing importance when introducing complementary views of the biological system.

Test Cases and Datasets

In this study, we use datasets extracted from the Genomic Data Commons (GDC) pan-cancer multi-omics study [START_REF] Grossman | Toward a shared vision for cancer genomic data[END_REF]. It is one of the most comprehensive datasets for multi-omics analysis, with high-dimensional omics data and corresponding phenotype data from The Cancer Genome Atlas (TCGA). In our experiments, we use three types of omics data, Copy Number Variations (CNV), RNA-Seq gene expressions, and DNA methylation. The CNV dataset comprises Gistic2 measurements on a total of 19,729 genes. The RNA-Seq expression dataset profile comprises around 60,484 identifiers referring to corresponding exons and measuring log2 transformed Fragments Per Kilobase of transcript per Million mapped reads (FPKM). Finally, the DNA methylation dataset was produced using the Infinium HumanMethylation450 BeadChip (450K) arrays with 485,578 probes in which beta values of probes indicate the methylation ratio of corresponding CpG sites. Moreover, we also evaluate our method on six smaller cohorts from TCGA: Bladder Urothelial Carcinoma (BLCA, n=437), Breast Invasive Carcinoma (BRCA, n=1022), Lung Adenocarcinoma (LUAD, n=498), Glioblastoma & Lower Grade Glioma (GBMLGG, n=515) and Uterine Corpus Endometrial Carcinoma (UCEC, n=538). We will perform in this study 4 different evaluations on 4 test cases for clasification and survival. The first task in our set of experiments is classifying the different tumor types in the pan-cancer study. The classification performance was measured using five metrics: Accuracy, Macro-averaged F1 score, precision, recall, and ROC-AUC [START_REF] Yang | Learning with multiclass auc: Theory and algorithms[END_REF]. We also perform a second classification task to validate our findings on a smaller dataset and test the robustness of our method. This task aims at performing a tumor subtype classification based on the PAM50 classification (IlluminlA, ILluminalB, Basal, and HER2), we use the same setup as the pan-cancer case. The third test case will be a survival study of the Pancancer dataset. Finally, the fourth test case will be about evaluating the survival performances on the 5 datasets presented earlier. For all those test cases, We compare the CustOmics model to several reference methods for multi-omics integration: first with a combination of dimensionality reduction methods, Multiple Factor Analysis (MFA), Uniform Manifold Approximation and Projection (UMAP) and non-negative matrix factorization (NMF), and also with different deep learning methods corresponding to the various strategies described in the Methods section.

Data Preprocessing

The data required some preprocessing before analysis.

• For RNA gene expression profiles, 594 exons located on the Y chromosome were removed, along with 1,904 ones with zero expression and 248 with missing values.

• For DNA methylation data, the same strategy as with gene expression profiles was used, in addition to removing probes that cannot be mapped to the human reference genome. It leaves us with 438,831 CpG sites.

Afterward, we intersected each combination of omics data in order to retrieve the maximum number of samples for each test case. We then identified and removed features with missing/consistently zero/NA values for other omic files.

Finally, for non-normalized datasets such as CNVs and RNA-Seq data, we applied a min-max normalization to ensure that each omic source was scaled identically and thus would have the same importance during integration.

Implementation Details

The CustOmics framework is based on the Pytorch deep-learning library. It can be applied to any combination of high-dimensional datasets and incorporates different integration strategies depending on the type of data and task to perform. As done in Zhang et al. [START_REF] Zhang | Omiembed: A unified multi-task deep learning framework for multi-omics data[END_REF], DNA methylation data can be divided into 23 separate blocks, each feeding a hidden layer corresponding to a chromosome to avoid overfitting and save GPU memory.

The whole architecture is built using fully-connected blocks. We use a batch normalization technique in each layer composing the neural network to address the internal covariate shift problemIoffe and Szegedy [2015]. Also, to avoid overfitting problems, we use dropout [START_REF] Srivastava | Dropout: A simple way to prevent neural networks from overfitting[END_REF]; its rate is considered a hyperparameter.

The input dataset was randomly split into training, validation, and testing set (60-20-20%) by using a stratified 5-fold cross-validation so that the proportion of samples in each tumor type between the different sets is preserved in all the folds. We perform Bayesian optimization [START_REF] Snoek | Practical bayesian optimization of machine learning algorithms[END_REF] using the validation set to find our model's best possible combination of hyperparameters.

All models were trained using an Nvidia Tesla V100S with 32 GB memory.

Results

Classification Results

We first perform the classification task on the pan-cancer dataset. Each architecture is coupled with an artificial neural network classifier composed of two hidden layers with 256 and 128 neurons. We see that the constructed layer representation succeeds at separating the data into 4 distinct clusters that we couldn't distinguish with each omic source alone. c. PAM50 gene importance: Computed SHAP values on the RNA-Seq data of the most relevant genes responsible for the discrimination between each subtype against the others using CustOmics for both integration phases.

Figure 6 and table 1 show the overall classification results (More details can be found in Table S1). Among the factorial methods, the MFA achieved the best results, so we coupled this method with the same ANN classifier used with the deep-learning representation methods as a basis for comparison. However, it does not perform as well as most deep-learning methods. It is because MFA cannot uncover nonlinear relationships between different sources, unlike deep-learning architectures. Moreover, as the MFA is an early integration, it suffers from the problems mentioned above when we introduced this strategy.

As hinted earlier, we can see that for the deep learning strategies, early integration is behind the others in terms of performance. It can be explained by the fact that RNA-Seq data hold more signals when determining tumor types or subtypes. Thus, concatenating the sources before feeding them to the VAE overshadows the other sources, and the learned representation depends mostly on RNA-Seq data without leveraging the other modalities. It is illustrated in Table S2, showing that the classification results using RNA-Seq data only are very close to those obtained with early integration, indicating that the model may overlook the interactions between sources. Late integration is not optimal either, since interactions between sources are not properly learned. Joint integration performs well in most cases, but we see that the best results displayed in Table S2 are achieved by the combination of only two sources, RNA-Seq and Methylation data, as it seems that CNV data only adds noise to the latent representation, meaning that its information is not handled well with this strategy.

Table 2: Classification performance for PAM50 classification on the TCGA-BRCA dataset, with 5 standard metrics.

We compare machine learning methods like UMAP, NMF, and MFA with deep-learning methods. We evaluate the performances on the final predicted output of the downstream classifier. These results confirm the interest of the CustOmics architecture, as it gives the best performances for all the test cases. Not only does it perform better than the others, but it also takes advantage of the complementarity and interactions between sources: As shown in Table S1, all sources bring an additional amount of information. Figure 6 gives a visualization of the different sources: Even though the initial sources are quite entangled, the CustOmics latent representation succeeds in separating the clusters using the mutual information between modalities, resulting in better classification performances. We also use the interpretability property of CustOmics introduced in the Methods section to highlight the most relevant features for the discrimination between PAM50 subtypes by computing their respective SHAP values for each source. We do it for both phases: in phase 1, we retrieve the relevant genes considered when using a single omic source, whereas, in phase 2, we investigate how the addition of other sources' signals changes the genes' importance. Figure 6 references the results of such explanations on RNA-Seq data, Figure S1 and S2 show the results for CNV and methylation data. We witness some well-referenced biomarkers for breast cancer like TFF1 [START_REF] Yi | Trefoil factor 1 (tff1) is a potential prognostic biomarker with functional significance in breast cancers[END_REF], suggesting that our method can retrieve relevant biological information.

Survival Analysis

The second task in this study is survival analysis. The goal is to predict the risk score associated with each patient from the corresponding high-dimensional omics data. Two standard metrics evaluate the performance of this downstream task: The c-index, which is a generalization of the AUC metric for censored data [START_REF] Harrell | Evaluating the Yield of Medical Tests[END_REF], and the Integrated Brier Score [START_REF] Graf | Assessment and comparison of prognostic classification schemes for survival data[END_REF].

Figure 4 and table 3 show the results for the different methods (more details can be found in Table S3) for the survival task. The same observations as for the classification task can be done regarding the differences between integration strategies. Here again, we also evaluated the performance of the CustOmics method for each combination of omics sources (Table S3).

The last task consists in evaluating the model performances for survival analysis for several specific cancer types of the TCGA datasets described in the dataset section. The objective is to evaluate the robustness of the models when dealing with smaller datasets.

Finally, we perform a more thorough analysis on the survival results that we display in Figure 4. We leave out 20% of our datasets for validation purpose, and we perform 5-fold cross-validation on the remaining 80% to compute the p-values associated to the log-rank test for different combinations of the Pan cancer test case. We show the ability of CustOmics to stratify the patients into distinct risk groups using the predicted hazard ratio. This stratification ability was also measured quantitatively using the p-value associated to the log-rank test between the different categories. Even Here again, our model outperforms the other integration strategies on both metrics. b. Log-rank test: We compute the p-value associated to the log-rank test between high and low risk groups for every integration strategy on a validation set for the pan-cancer survival test case, and we compare it to mono-omic survival predictions. c. Kaplan Meier Curves : We draw the Kaplan Meier curves and display the p-value associated to the log-rank test as computed previously for each dataset using the predicted hazard from the CustOmics model and stratify the population into high and low risk on the test set for the predicted hazard ratio. This figure shows that our method succeeds in stratifying the patients into risk subgroups. though the comparison between joint, late and mixed integration is not significant in this case, it is interesting to note that the addition of multiple omics sources has greatly affected the p-value as it was nearly multiplied by a factor 10 -5 for CNV data and 10 -2 for RNA-Seq and methylation data. Those results also show that the early integration is the strategy with minimal enrichment from other sources, which corroborates our previous intuition and is coherent with the results found in the different experiments.

Discussion

In this work, we presented a range of integration strategies for multi-source data that can handle both the high dimensionality and the heterogeneity of the data. To take the best of all those strategies, we presented the mixedintegration and the CustOmics framework to alleviate the limitations of the existing methods. This new framework can achieve better latent representations and lead to a more robust and generalizable architecture, as shown with the systematic better results than alternative strategies. Importantly, our method is capable of adapting to each omic source by handling the training independently in a first phase, which solves the issue of unbalanced signals between the sources by standardizing the representations before starting to learn cross-modality interactions. Our fusion model demonstrated better performances for both classification and survival outcome prediction across all test cases. It not only achieved great results on pan-cancer data, but it also did on smaller dataset for particular cohorts, showing the robustness of our method in situations with fewer samples. Interestingly, the performances of CustOmics in comparison to other integration strategies for different combinations of omic sources show that the mixed-integration handles better the fusion between CNV, RNA-Seq and methylation data. This aspect could probably be further enhanced by adding prior knowledge of those sources in the intermediate autoencoders, for example introducing a negative binomial prior in the RNA-Seq autoencoder. Furthermore, by adapting the SHAP method to our architecture, we were able to highlight genes of importance for specific tasks. However, there is room for improvement from a computational point of view, since the interpretability method is very expensive.

As clinical applications often suffer from a lack of samples to build efficient deep learning models and as all the modalities are not always extracted for all patients, it would be beneficial to study the benefits of per-source transfer learning during phase 1 to pick-up weaker signals and eliminate the noise. In fact, one of the advantages of phase 1 learning the sub-representations independently is that it is not necessary for all patients to have all omic sources available at all times, which helps us take the best advantage of the data at hand.

Conclusion

In conclusion, our generic and interpretable multi-source deep learning framework improves on state-of-the-art integration strategies by proposing a hybrid approach that fits well with multi-omics data. The framework is available to use on Github: https://github.com/HakimBenkirane/CustOmics. Possible steps to refine CustOmics are a more thorough study on the interpretability part of the framework and adding prior knowledge to individual autoencoders to enhance the sub-representations of omic sources during phase 1 of mixed-integration.

Supplementary Materials

Classification Results

Figure 1 :

 1 Figure 1: a. Early Integration: Sources are concatenated before being fed to a single model. b. Joint Integration:Sub-representations of each source are learned jointly before inputting into the model. c. Late Integration: Each of the sources outputs its prediction using its independent model, the predictions are then aggregated. d. Mixed Integration: Representation of the mixed integration approach. In phase 1, a specific model is trained for each source independently and embeds a sub-representation adapted to the source's specificities. In phase 2, the specific models are trained jointly to create the final output.

Figure 2 :

 2 Figure 2: a. Early Integration VAE: Variational Autoencoder architecture with early integration strategy. b. Joint Integration VAE: Variational Autoencoder architecture with joint integration strategy. c. Late Integration VAE: Variational Autoencoder architecture with late integration strategy.d. Mixed-Integration/CustOmics: This is a hierarchical architecture composed of specific per-source autoencoders that converges into a central variational autoencoder.

Figure 3 :

 3 Figure 3: a. pan-cancer and PAM50 classification results: Overall classification results for the pan-cancer tumor classification test case and the PAM50 subtype classification for breast cancer. b. T-SNE vizualization for each omic source separately, along with the latent representation constructed by CustOmics.We see that the constructed layer representation succeeds at separating the data into 4 distinct clusters that we couldn't distinguish with each omic source alone. c. PAM50 gene importance: Computed SHAP values on the RNA-Seq data of the most relevant genes responsible for the discrimination between each subtype against the others using CustOmics for both integration phases.

Figure 4 :

 4 Figure4: a. Survival Analysis Performances: We evaluate the performances of the survival model for the pan-cancer dataset using both C-index and the Integrated Brier Score (IBS). Here again, our model outperforms the other integration strategies on both metrics. b. Log-rank test: We compute the p-value associated to the log-rank test between high and low risk groups for every integration strategy on a validation set for the pan-cancer survival test case, and we compare it to mono-omic survival predictions. c. Kaplan Meier Curves : We draw the Kaplan Meier curves and display the p-value associated to the log-rank test as computed previously for each dataset using the predicted hazard from the CustOmics model and stratify the population into high and low risk on the test set for the predicted hazard ratio. This figure shows that our method succeeds in stratifying the patients into risk subgroups.

Figure 6 :

 6 Figure6: PAM50 gene importance: Computed SHAP values on the methylation data of the most relevant genes responsible for the discrimination between each subtype against the others using CustOmics for both integration phases.

Table 1 :

 1 The classification performance for the pan-cancer dataset is evaluated with 5 standard metrics for UMAP, NMF, MFA, and deep-learning methods. We evaluate the performances on the final predicted output of the downstream classifier.

	Model	Accuracy	F1-score	Precision	Recall	ROC-AUC
	UMAP	0.7598 ± 0.0036	0.7149 ± 0.0029	0.7200 ± 0.0031	0.7598 ± 0.0032	0.8740 ± 0.0012
	NMF	0.8599 ± 0.0017	0.8406 ± 0.0013	0.8460 ± 0.0018	0.8599 ± 0.0021	0.9266 ± 0.0019
	MFA	0.9167 ± 0.0012	0.9025 ± 0.0014	0.8945 ± 0.0008	0.9167 ± 0.0013	0.9565 ± 0.0003
	Early Int. VAE	0.9337 ± 0.0079	0.9314 ± 0.0086	0.9367 ± 0.0067	0.9337 ± 0.0079	0.9655 ± 0.0041
	Joint Int. VAE	0.9610 ± 0.0032	0.9600 ± 0.0032	0.9631 ± 0.0043	0.9610 ± 0.0032	0.9898 ± 0.0005
	Late Int. VAE	0.9492 ± 0.0115	0.9464 ± 0.0111	0.9498 ± 0.0079	0.9492 ± 0.0115	0.9737 ± 0.0060
	CustOmics	0.9788 ± 0.0025 0.9705 ± 0.0033 0.9728 ± 0.0041 0.9685 ± 0.0034 0.9918 ± 0.0001

Table 3 :

 3 The survival analysis performance for the pan-cancer dataset is evaluated with 2 standard metrics, C-index and IBS. We compare classical methods like UMAP, NMF, and MFA with deep-learning methods and evaluate the performances on the final predicted output of the downstream survival network.

	Model	C-index	IBS
	UMAP	0.5948 ± 0.0231	0.2486 ± 0.0327
	NMF	0.6012 ± 0.0204	0.2207 ± 0.0264
	MFA	0.6127 ± 0.0164	0.2192 ± 0.0203
	Early Int. VAE	0.6578 ± 0.0103	0.2106 ± 0.0117
	Joint Int. VAE	0.6709 ± 0.0041	0.1802 ± 0.0072
	Late Int. VAE	0.6629 ± 0.0086	0.2112 ± 0.0088
	CustOmics	0.6841 ± 0.0033 0.1745 ± 0.0052

Table 4 :

 4 Performances of state-of-the-art integration methods for survival analysis, using concordance index on 5 TCGA cohorts: Bladder Urothelial Carcinoma (BLCA), Breast Invasive Carcinoma (BRCA), Glioblastoma & Lower Grade Glioma (GBMLGG), Lung Adenocarcinoma (LUAD) and Ulterine Corpus Endometrial Carcinoma (UCEC).

	Model	BLCA	BRCA	GBMLGG	LUAD	UCEC	Overall
	UMAP	0.527 ± 0.048	0.524 ± 0.039	0.557 ± 0.028	0.530 ± 0.017	0.543 ± 0.025	0.536
	NMF	0.553 ± 0.060	0.560 ± 0.036	0.584 ± 0.067	0.589 ± 0.052	0.570 ± 0.068	0.571
	MFA	0.591 ± 0.052	0.599 ± 0.043	0.605 ± 0.036	0.597 ± 0.058	0.586 ± 0.020	0.596
	Early Int. VAE 0.603 ± 0.054	0.618 ± 0.039	0.628 ± 0.056	0.612 ± 0.041	0.609 ± 0.032	0.614
	Joint Int. VAE 0.616 ± 0.072	0.627 ± 0.030	0.635 ± 0.020	0.608 ± 0.038	0.630 ± 0.021	0.624
	Late Int. VAE	0.610 ± 0.055	0.620 ± 0.057	0.621 ± 0.067	0.595 ± 0.012	0.627 ± 0.022	0.615
	CustOmics	0.637 ± 0.050 0.633 ± 0.018 0.642 ± 0.028 0.625 ± 0.037 0.667 ± 0.022 0.640

Table 5 :

 5 Classification performances for multiple combinations of omics data using CustOMICS on the Pancancer dataset. We can see that the best performances are brought by Transcriptomics data, but the addition of other omics data increases the performances, showing that the integration is relevant. ± 0.01 0.94 ± 0.01 0.94 ± 0.01 0.94 ± 0.0132 0.97 ± 0.00 CNV + RNAseq + methyl 0.98 ± 0.01 0.97 ± 0.01 0.97 ± 0.01 0.97 ± 0.01 0.99 ± 0.00

	Omics	Accuracy	F1-score	Precision	Recall	ROC-AUC
	CNV	0.47 ± 0.03 0.47 ± 0.03 0.47 ± 0.03	0.48 ± 0.01	0.75 ± 0.02
	RNAseq	0.92 ± 0.01 0.93 ± 0.01 0.93 ± 0.01	0.92 ± 0.01	0.96 ± 0.00
	methyl	0.68 ± 0.02 0.68 ± 0.02 0.68 ± 0.02	0.68 ± 0.02	0.82 ± 0.01
	CNV + RNAseq	0.93 ± 0.01 0.92 ± 0.01 0.92 ± 0.01	0.92 ± 0.01	0.96 ± 0.00
	CNV + methyl	0.71 ± 0.02 0.69 ± 0.02 0.70 ± 0.02	0.69 ± 0.02	0.85 ± 0.01
	RNAseq + methyl	0.94

Table 6 :

 6 Classification performances for multiple combinations of omics data using joint integration on the Pancancer dataset. ± 0.01 0.96 ± 0.01 0.96 ± 0.01 0.96 ± 0.0132 0.99 ± 0.00 CNV + RNAseq + methyl 0.94 ± 0.01 0.94 ± 0.01 0.95 ± 0.01 0.94 ± 0.01 0.97 ± 0.00

	Omics	Accuracy	F1-score	Precision	Recall	ROC-AUC
	CNV	0.47 ± 0.03 0.47 ± 0.03 0.47 ± 0.03	0.48 ± 0.01	0.75 ± 0.02
	RNAseq	0.92 ± 0.01 0.93 ± 0.01 0.93 ± 0.01	0.92 ± 0.01	0.96 ± 0.00
	methyl	0.68 ± 0.02 0.68 ± 0.02 0.68 ± 0.02	0.68 ± 0.02	0.82 ± 0.01
	CNV + RNAseq	0.90 ± 0.02 0.89 ± 0.02 0.90 ± 0.02	0.88 ± 0.03	0.93 ± 0.00
	CNV + methyl	0.70 ± 0.02 0.69 ± 0.02 0.70 ± 0.02	0.70 ± 0.02	0.85 ± 0.01
	RNAseq + methyl 0.96 Survival Results

Table 7 :

 7 Survival performances for multiple combinations of omics data using CustOMICS on the Pan cancer dataset. We can see that the best performances are brought by Transcriptomics data, but the addition of other omics data increases the performances, showing that the integration is relevant.

	Omics	Accuracy	F1-score
	CNV	0.54 ± 0.05 0.25 ± 0.06
	RNASeq	0.63 ± 0.02 0.20 ± 0.02
	methyl	0.59 ± 0.02 0.23 ± 0.03
	CNV + RNAseq	0.64 ± 0.02 0.19 ± 0.02
	CNV + methyl	0.61 ± 0.03 0.21 ± 0.03
	RNAseq + methyl	0.64 ± 0.02 0.19 ± 0.02
	CNV + RNAseq + methyl 0.68 ± 0.01 0.17 ± 0.01


Aude Lagier
email: aude.lagier@univ-amu.fr

Daphné Guenoun

Thierry Legou

Robert Espesser

Antoine Giovanni

Pierre Champsaur

Control of the glottal configuration in ex vivo human models: quantitative anatomy for clinical and experimental practices

Keywords: Larynx, Glottis, Vocal folds adduction, Phonation control, Voice

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The human larynx is a complex structure of several cartilages mobilized by small muscles called intrinsic laryngeal muscles. The movements of the vocal folds mainly result of the mobilization of a set of two complex joints (i.e., the cricoarytenoid joints), under the action of the intrinsic laryngeal muscles. Prades et al. studied the biomechanics of the cricoarytenoid joint and showed that the adduction of the vocal process (to close the glottis) combines anterior tilt, superomedial shift, and a lesser degree of internal rotation of the cricoarytenoid joint [START_REF] Prades | Descriptive anatomy of the cricoarytenoid articulation: application to articular dynamics in carcinology[END_REF]. Understanding the role of the intrinsic muscles in vocal folds movements is of interest in physiological and clinical fields. Currently, three groups of intrinsic laryngeal muscles are defined: adductors for vocal fold closing, abductors for vocal fold opening, and lengtheners for increasing vocal folds length. The latter group also induces vocal fold tension, especially when they contract simultaneously with muscles that shorten the vocal folds [START_REF] Poletto | Correspondence between laryngeal vocal fold movement and muscle activity during speech and nonspeech gestures[END_REF]. Determining the role of each laryngeal muscle is possible in excised larynges [START_REF] Berry | An investigation of cricoarytenoid joint mechanics using simulated muscle forces[END_REF], and the muscular actions are consistent with the cricoarytenoid joint degrees of freedom: the anterior-posterior translation is the largest movement; medial-lateral and vertical glidings are less important [START_REF] Berry | An investigation of cricoarytenoid joint mechanics using simulated muscle forces[END_REF]. During phonation, EMG studies showed that all adductor muscles (thyroarytenoid-TA, lateral cricoarytenoid-LCA, arytenoid-IA) are active. Cricothyroid-CT and TA tend to co-contract during pitch elevation, most likely stiffening the vocal fold due to their opposite actions of lengthening and shortening the vocal fold [START_REF] Poletto | Correspondence between laryngeal vocal fold movement and muscle activity during speech and nonspeech gestures[END_REF][START_REF] Titze | Role of the thyroarytenoid muscle in regulation of fundamental frequency[END_REF]. In summary, studies on excised larynges trying to imitate the isolated action of each intrinsic muscle may not accurately predict the biomechanical effect of intact, normal human intrinsic muscles, because the cocontraction of the different muscles may lead to different resultant movements.

From the point of view of the resultant movement, the prephonatory glottal configuration is the key point. The prephonatory glottal configuration is determined by the degree of adduction of the folds and by the viscoelasticity of the vocal folds tissues. The three viscoelastic-related physical properties are mass, stiffness, and viscosity [START_REF] Jiang | Vocal fold physiology[END_REF]. Adduction and viscoelasticity of the vocal folds are very important in the quality of the voice, in the control of the fundamental frequency, and in the relationship between fundamental frequency and subglottal pressure [START_REF] Alipour | On pressure-frequency relations in the excised larynx[END_REF][START_REF] Jiang | Vocal fold physiology[END_REF]. Improving the prephonatory glottal configuration is the main objective in the surgery of unilateral laryngeal paralysis: the better the adduction, the better the voice quality [START_REF] Hoffman | Multiparameter comparison of injection laryngoplasty, medialization laryngoplasty, and arytenoid adduction in an excised larynx model[END_REF].

The objective of this paper was to identify the determining factors of the glottal prephonatory configuration from the point of view of the resulting muscular actions (i.e., arytenoids adduction, membranous vocal fold adduction, and tension) instead of the action of each intrinsic muscle.

Materials and methods

Larynges

Twenty-one human larynges (12 females and 9 males) were harvested from fresh cadavers, within 24-48 h after death, according to French law. At the time of death, the age of the persons ranged from 71 to 92 years. No history of voice disorders, thoracic or cervical cancer or surgery was known for any subject, and no abnormalities were observed on the larynges. Larynges were rapidly frozen in 0.9 % saline. They were thawed overnight in saline solution and cleaned with removal of the extraneous tissues and supraglottic structures (epiglottis and vestibular folds) before the experiment, according to the literature [START_REF] Alipour | Phonatory characteristics of the excised human larynx in comparison to other species[END_REF].

Experiments Experiment A: adduction

Eleven larynges (7 females and 4 males) were tested in the first part of the experiment, dedicated to the adduction.

Each larynx was tested in the four following conditions in a random order (Fig. 1):

• No adduction (NA) was the spontaneous position of the vocal folds • Arytenoid adduction (AA): the action of adductor muscles of the arytenoid process (i.e., lateral cricoarytenoid and inter-arytenoid muscles) was reproduced: the maximal adduction was obtained with a stitch between the vocal processes of the arytenoid cartilages. • Vocal fold medialization with Montgomery implant (MG): the action of adductor muscles of the membranous vocal fold (i.e., thyroarytenoid muscle) was reproduced according to type I thyroplasty [START_REF] Isshiki | Vocal mechanics as the basis for phonosurgery[END_REF][START_REF] Isshiki | Thyroplasty as a new phonosurgical technique[END_REF] bilaterally performed with two montgomery implants [START_REF] Montgomery | Montgomery thyroplasty implant system[END_REF] inserted through bilateral windows into the thyroid ala. • Double adduction (DA): combination of arytenoid adduction and Montgomery implants. Ten larynges (5 females and 5 males) were tested in the second part, testing the lengthening of vocal folds. In this experiment, all the larynges were set in double adduction.

The resultant action of the tensor muscles (i.e., cricothyroid muscles) was reproduced with cricothyroid approximation. A wire was attached on the midline, from the arch of the cricoid cartilage to the inferior edge of the thyroid cartilage (Fig. 2). The wire was connected to a micrometer and the approximation stepped 1 mm to 1 mm.

A wire was placed on the external aspect of the cricoid and another one on the thyroid cartilage as markers for measures and comparison of cricothyroid angle when cricoid and thyroid approximated. Due to the variability of the insertions of these two wires, only the variations of the cricothyroid angle were measured.

The length of the vocal folds was measured from the end of the vocal process to the anterior commissure for each level of cricothyroid approximation.

Measures

Two cameras were placed orthogonally around the larynx. The first one (Canon Ò 600D, lens EFS 18-135 mm) was in the axis of the larynx, at 50-cm distance to take pictures of the superior aspect of the vocal folds. The second one (Panasonic Ò Lumix F535) was on the left side of the larynx, at 50-cm distance to take pictures of the CT angle.

The open-access software Image J (http://rsb.info.nih.gov/ ij/) was used for the measurement of glottal area in experiment A, and of the vocal fold length and cricothyroid angle in experiment B was performed using. Each picture was taken with a measurer to calibrate the pixels/millimeter ratio.

In experiment A, the mean of three measures of the glottal area was calculated for each larynx in the four conditions. In experiment B, the mean of three measures of CT angle and of vocal folds length was calculated for each larynx in the 25 conditions.

Statistical analysis

The descriptive statistics are expressed as median, minimum, and maximum.

The analyses were conducted using linear mixed models to account for the grouping structure (by larynx) of the data. Package lme4 [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF] was used in R software [START_REF] Core | R: a language and environment for statistical computing[END_REF]. P values were computed using a parametric bootstrap procedure [START_REF] Halekoh | A Kenward-Roger approximation and parametric bootstrap methods for tests in linear mixed models-the R package pbkrtest[END_REF] and noted as P boot ; the significance of statistical threshold was set at 0.05.

Results

Experiment A: adduction

The mean glottal area significantly decreased from 41.2 mm 2 mean with no adduction, to 10.2 mm 2 mean with arytenoid adduction, to 9.2 mm 2 with membranous vocal fold adduction, and down to 1.1 mm 2 with the combination of arytenoid and membranous adduction.

The dependent variable was the log-transformed glottal area. The log-transformation enhanced the normality of the distribution of the residuals of the model. The predictors were the condition (NA, AA, MG, and DA), the sex and their interactions. A random intercept was added to account for the variability across the 11 larynges. Forty observations were available (the log-transformation implied that four null area were removed from the analysis). The interaction between condition and sex was not significant (P boot = 0.92), and it was removed from the final model.

The comparisons of interest needed the estimation of two variants of the model, and consequently, the significance threshold was Bonferroni-corrected at P \ 0.025. The results showed that the area for the AA condition was significantly smaller than for the NA condition (b =-1.508, t =-4.55, P boot \ 0.001); the median decreased by 22 %; the area for the AA and MG conditions could not be distinguished (b =-0.275, t =-0.83, P boot = 0.44); the area for the DA condition was significantly smaller than for the MG condition (b =-1.412, t =-3.75, P boot \ 0.001); the median decreased by 25 %; the area for the DA condition was also significantly smaller than for the AA condition (b =-1.686, t =-4.49, P boot \ 0.001); the median decreased by 19 %. As the interaction with sex was not significant, this configuration was the same for both male and female, the area for male being greater than for female by a factor 2.3 (b = 0.834, t =-3.30, P boot = 0.002). Figure 3 shows the predicted area by the model and the boxplots of the observed length on the logarithm scale.

Experiment B: lengthening

The length of vocal folds increased along with the cricothyroid approximation from 13.61 mm median (min 8.46 mm; max 16.09 mm) up to 14.48 mm median (min 9.10 mm; max 17.65 mm). The lengthening of the vocal folds was 11.6 % (min 4.9 %; max 27.6 %). The maximum vocal fold length was achieved for a median of 22 mm of CT approximation (min 16 mm, max 25 mm)

Examination of the data suggested that the evolution of the length of the vocal folds would be best modelized with a quadratic function of the CT approximation.

The dependent variable was the vocal folds length. The predictors were the CT approximation (linear and quadratic component), the sex and their interactions. A random intercept and random slopes for each component of the CT approximation were added to account for the variability across the ten larynges. All the random terms were uncorrelated. Two hundred and thirty-four observations were available.

The interaction between sex and CT approximation was not significant (Likelihood ratio test with KR approximation: F test = 0.044, df = 2, p = 0.95). This interaction was removed from the final model. The results showed that the quadratic term (b =-00033, t =-5.5, P boot \ 0.001) and the linear term (b = 0.13, t = 6.86, P boot \ 0.001) were significant. As expected the VF length was significantly longer for the males (b = 4.02, t = 4.02, P boot \ 0.001). Figure 4 shows the predicted length by the model and the boxplots of the observed lengths.

The cricothyroid angle (CT angle) decreased of 10.05°m edian (min -5.42°, max -14.25°) along with cricothyroid approximation. The maximum decrease in CT angle was achieved for a median of 22 mm of CT approximation (min 16 mm, max 25 mm).

Examination of the data suggested that the evolution of the CT angle would be best modelized with a quadratic function of the CT approximation.

The variable dependent was the CT angle. The predictors were the CT approximation (linear and quadratic component), the sex and their interactions. A random intercept and random slopes for each component of the CT approximation were added to account for the variability across the ten larynges. All the random terms were uncorrelated. Two hundred and forty-six observations were available.

The interaction between sex and CT approximation was not significant (P boot = 0.92). This interaction was removed from the final model. The results showed that the quadratic term (b = 0.0115, t = 4.63, P boot \ 0.001) and the linear term (b =-0.697, t =-10.12, P boot \ 0.001) were significant. The predictor sex was not significant (b =-0.912, t =-1.41, P boot = 0.52). Figure 5 shows the boxplots of the observed angles and the angles predicted by the model.

Discussion

The results of experiment A emphasize the subdivision of the group of adductor intrinsic muscles in arytenoids adductors (i.e., LCA and IA), and membranous vocal fold These results are interesting for the methodology of physiological studies on excised larynges as they emphasize the consequences of the different techniques of adduction: the reproduction of the actions of both arytenoid and vocal fold adductors (double adduction) leads to the narrower glottal area. This double adduction is also nearer to the reality, as shown in the in vivo literature based on electromyographic studies. Electromyographic routine uses some laryngeal gestures to test the validity of proper electrode placement: activation of TA increases during sustained vowel; LCA activation is a burst at onset and offset of phonation but not a sustained activation during prolonged phonation; the CT activation increases during the high end of a pitch glide [START_REF] Hirano | Use of hooked-wire electrodes for electromyography of the intrinsic laryngeal muscles[END_REF][START_REF] Poletto | Correspondence between laryngeal vocal fold movement and muscle activity during speech and nonspeech gestures[END_REF].

This study aimed at creating more physiological prephonatory glottal configuration, but not at achieving the very precise control used in real voice: in vivo, glottal configuration may differ with the phonetic context [START_REF] Poletto | Correspondence between laryngeal vocal fold movement and muscle activity during speech and nonspeech gestures[END_REF], but also with the expressive intention or with the different singing modes [START_REF] Freeman | A comparison of sung and spoken phonation onset gestures using high-speed digital imaging[END_REF].

The results of experiment B clearly showed direct action of cricothyroid approximation on the vocal fold length. In vivo canine experiments with stimulation of the laryngeal superior nerve found similar results, but the lengthening was much higher in dogs with a maximum range of length change 44.7 % mean [START_REF] Titze | The dynamics of length change in canine vocal folds[END_REF] vs 11.6 % mean in this study.

The choice of the range of the CT approximation (0-25 mm) was arbitrarily chosen. The parabolic aspect of the curves highlights that this range was probably too large, and the thyroid could not continue its course toward the cricoid, because they were in contact.

The two main limits of the study are the low number of larynges and the inability to reproduce the TA isometric contraction. The number of larynges studied was counterbalanced by the high number of measures and the intralarynx and inter-larynx reproducibility of responses to the tasks in the two experiments.

The design of the study, based on excised larynges, could not reproduce isometric longitudinal contraction of the TA. According to the cover-body theory of vocal folds vibration [START_REF] Hirano | Morphological structure of the vocal cord as a vibrator and its variations[END_REF], physical quantities of the vocal folds, such as length, mass, and stiffness, are the physical quantities that determine the control of fundamental frequency. These quantities depend on the complex co-activations of CT and TA muscles. CT activation is always associated with an increase in fundamental frequency. For given levels of CT activation, the effects of TA activation on F0 are different [START_REF] Lowell | Simulated effects of cricothyroid and thyroarytenoid muscle activation on adult-male vocal fold vibration[END_REF]. Simulations on 3D numerical models of vocal folds show that these co-contractions impact the vocal fold morphology and the stress distribution, and must be one of the major determinants in the effective depth of tissue vibration [START_REF] Deguchi | Cooperative regulation of vocal fold morphology and stress by the cricothyroid and thyroarytenoid muscles[END_REF]. However, this component of vocal fold physiology could not be assessed here, and the role of the TA in the control of vocal fold stiffness could not be properly studied. The results presented above are of clinical relevance: the surgical management of unilateral laryngeal paralysis (excluding re-innervation techniques) can propose injection laryngoplasty, medialization laryngoplasty (with montgomery implant or others), and finally the association of adduction arytenopexy with medialization laryngoplasty as described by Zeitels et al. It is noteworthy that adduction arytenoidopexy is always associated with membranous vocal fold medialization [START_REF] Zeitels | Adduction arytenopexy for vocal fold paralysis: indications and technique[END_REF], and is logically performed when the vocal improvement of the vocal fold medialization is not sufficient.

The excised larynx is a very convenient tool to teach the classical surgical techniques on the larynx, but also for the tests and comparisons of different techniques of laryngeal surgery [START_REF] Hoffman | Multiparameter comparison of injection laryngoplasty, medialization laryngoplasty, and arytenoid adduction in an excised larynx model[END_REF][START_REF] Isshiki | Vocal mechanics as the basis for phonosurgery[END_REF]. Hoffman et al. compared the different surgical techniques for treatment of unilateral laryngeal paralysis in canine excised larynges, and concluded that the association of medialization laryngoplasty (i.e., with montgomery implants) and arytenoid adduction leads to the greatest improvement in phonatory parameters. The medialization laryngoplasty alone did not create the arytenoid adduction. They also showed that, contrary to injection laryngoplasty, medialization laryngoplasty did not cause detrimental changes to the mucosal wave [START_REF] Hoffman | Multiparameter comparison of injection laryngoplasty, medialization laryngoplasty, and arytenoid adduction in an excised larynx model[END_REF].

In conclusion, the results of this study emphasize that the adductor group of intrinsic laryngeal muscles are functionally splitted in arytenoid and membranous vocal fold adductors. This has clinical implications in the choice of the surgical technique for unilateral laryngeal paralysis, but also for the methodology of experiments in vocal physiology on excised larynges. The action of the cricothyroid muscle on the vocal fold length was confirmed, but some significant difference in the range of the lengthening has been highlighted in comparison to animal studies. The main limit of the study was the inability to reproduce the TA longitudinal isometric contraction and, therefore, the inability to reproduce the whole control mechanisms of fundamental frequency.

Compliance with ethical standards

Conflict of interest The authors declare that there is no conflict of interest of any kind.

Fig. 1

 1 Fig. 1 Superior view of the larynx in experiment A: four conditions of adduction

Fig. 2

 2 Fig. 2 Left lateral view of the larynx set for experiment B. (1) Montgomery implant, (2) wire responsible for cricothyroid approximation, (3) cricoid wire for CT angle measurement, and (4) thyroid wire for CT angle measurement

Fig. 3 Fig. 4

 34 Fig. 3 Glottal area and the four conditions of adduction (NA no adduction, AA arytenoid adduction, MG membranous fold medialization, DA double adduction). The boxplots show the observed data, the colored point the area predicted by the model (with the fixed effects only). The colored lines indicate the 95 % confidence interval (color figure online)

Fig. 5

 5 Fig. 5 CT angle according the CT approximation and the sex. The boxplots show the observed data, the colored lines the angle predicted by the model (with the fixed effects only). The colored bands indicate the 95 % confidence interval (color figure online)


ITOLS Michael Wettstein
email: m.wettstein@itols.ch

Elyazid Mouhsine

Jean-Manuel Aubaniac

Laurent Audigé

Matthieu Ollivier

Pierre-François Leyvraz

Jean-Noël Argenson

Jean-Noël Argenson4

The torsion of the proximal femur in cementless total hip arthroplasty: a 3-dimensional evaluation

Keywords: Anatomy, arthroplasty, cementless, femur, helitorsion, hip, torsion

diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The anatomy of the proximal femur has been widely studied, aiming to standardise the conception of cementless femoral stems. Noble et al., [START_REF] Noble | The anatomic basis of femoral component design[END_REF] Rubin et al. [START_REF] Rubin | The morphology of the proximal femur. A three-dimensional radiographic analysis[END_REF] and Laine et al. [START_REF] Laine | Diversity of proximal femoral medullary canal[END_REF] showed that the anatomy of the proximal femur can be divided into intramedullary and extramedullary parts.

The extramedullary anatomy of the femur setts the position of the hip centre, depending on neck length, frontal offset, CCD angle (Centre -Colum -Diaphyseal) and neck anteversion. Despite the clinical importance of this extramedullary part, [START_REF] Ramaniraka | Noncemented total hip arthroplasty: influence of extramedullary parameters on initial implant stability and on boneimplant interface stresses[END_REF][START_REF] Sugano | Predicting the position of the femoral head center[END_REF][START_REF] Lecerf | Femoral offset: anatomical concept, definition, assessment, implications for preoperative templating and hip arthroplasty[END_REF][START_REF] Matsushita | Modular necks improve the range of hip motion in cases with excessively anteverted or retroverted femurs in THA[END_REF][START_REF] Mcmahon | The accuracy of acetabular cup positioning using patient-specific 3D orientation guidance in total hip arthroplasty[END_REF] the present study will focus mainly on the intramedullary anatomy, with consequences on neck anteversion.

In an anatomic study of 200 cadaveric femurs, Noble et al. [START_REF] Noble | The anatomic basis of femoral component design[END_REF] defined canal flare index as the width of the canal 20 mm above the lesser trochanter divided by the width at the level of the isthmus, finding 83% normal, 9% called stovepipe (large canal) and 8% champagne-fluted (narrow canal).

Rubin et al. [START_REF] Rubin | The morphology of the proximal femur. A three-dimensional radiographic analysis[END_REF] showed the same anatomic variations of the proximal femur in the frontal plane, highlighting that a standard radiological analysis was insufficient and recommending an additional CT-scan.

Husmann et al., [START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF] in a more recent CT-scan study, showed that the variations in shape described by Noble et al. [START_REF] Noble | The anatomic basis of femoral component design[END_REF] do not only occur in the frontal plane but also in the sagittal plane. When considering canal flare index in both planes they only found 30% of the femurs could be considered as normal.

In contrast with its frontal and sagittal shape, the geometry of the proximal femur in the horizontal plane is not well described, except for neck anteversion. [START_REF] Hackenbruch | Die idiopathische coxa antetorta (antetorsionssyndrom)[END_REF][START_REF] Haspl | Assessment of femoral neck-shaft and antetorsion angles[END_REF][START_REF] Starker | Validation of computerized tomography antetorsion angle measurement of the femur[END_REF][START_REF] Wissing | Determination of rotational defects of the femur by computer tomographic determination of the antetorsion angle of the femoral neck[END_REF][START_REF] Tönnis | Decreased acetabular anteversion and femur neck antetorsion cause pain and arthrosis. 1: statistics and clinical sequelae[END_REF] This anterior torsion of the femoral neck is defined as the angle between the posterior bicondylar plane and the axis of the femoral neck. According to Husmann et al., [START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF] its mean value is 24.7° ± 8.7°.

Wilkinson, [START_REF] Wilkinson | Femoral anteversion in the rabbit[END_REF] and more recently Yu et al., [START_REF] Yu | Prediction of postoperative stem anteversion in Crowe type II/III developmental dysplasia of the hip on preoperative two-dimensional computed tomography[END_REF] have shown that anterior torsion does not only depend on the position of implantation of the femoral neck on the intertrochanteric region, but also on the torsion of the proximal femoral metaphysis.

Proximal femoral torsion is defined as the angle between the posterior femoral bicondylar plane and the axis of the femoral metaphysis at the level of osteotomy of the femoral neck for total hip arthroplasty (THA). According to Husmann et al., [START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF] proximal femoral torsion, or helitorsion, has no correlation with anterior torsion of the femoral neck and shows broad individual variations between 8.3° and 44.9°, with a mean value of 24.0° ± 11.8°.

However, the choice of the level of proximal femoral torsion corresponding to the level of osteotomy of the femoral neck during hip arthroplasty does not take into account the progression of metaphyseal torsion along the proximal femur. Furthermore, this could explain why anatomic and prosthetic antetorsion may not match, as mentioned by Sendtner at al. [START_REF] Sendtner | Stem torsion in total hip replacement: CT measurements in 60 patients[END_REF] Because of the recent interest in metaphyseal fixation of cementless stems, particularly in so-called "anatomic designs", a more precise knowledge of the evolution of proximal femoral torsion seems important. This should help with the design of cementless stems in THA for various anatomies and also allow a precise adaptation of the stem to the intramedullary anatomy of each patient when needed, in order to achieve an optimal fit-and-fill. [START_REF] Callaghan | The clinical results and basic science of total hip arthroplasty with porous-coated prostheses[END_REF][START_REF] Rubin | Frictional interface micromotions and anisotropic stress distribution in a femoral total hip component[END_REF][START_REF] De Boer | Comparison of anatomic vs. straight femoral stem design in total hip replacement -femoral canal fill in vivo[END_REF] Our hypothesis was that the femoral torsion progressively changes along the proximal femoral diaphysis and metaphysis and that its progression is highly individual. The goal of this study was to determine the mean values of proximal femur torsion at different levels of the metaphysis as well as its individual differences in various types of pathologic hips found at the time of THA.

Material and methods

A total of 178 patients scheduled for cementless THA and representing various etiologies were evaluated by X-rays and CT-scans using the same protocol. The X-rays consisted of A/P pelvis, both hips A/P and profile, as well as a telemetry view. CT-scans were carried out in 5 millimetre intervals from the top of the acetabulum to the femoral isthmus, 1 slice at the knee and 1 at the foot. Because the femur was not scanned along its complete length, it was not possible to relate the level of the slices to the relative length of the femur.

Table 1 shows the number of cases analysed for each diagnosis as well as sex and age distribution.

Primary hip osteoarthritis (OA), used as reference values, is defined as osteoarthritis with coxometry values within a normal range. [START_REF] Gordenne | Results of our experience with coxometry (apropos of 104 examinations). Verification of normal angular values in the adult[END_REF] Dysplasia is defined as Wiberg 22 VCE (vertical centre-edge) angle <25°, Tschauner AC (acetabular index) angle >10° and Lequesne 23 ACE (anterior centre-edge) angle <25°. We used the expression hip dysplasia for Crowe stage I, II and III and the expression congenital dislocation of the hip for Crowe et al. [START_REF] Crowe | Total hip replacement in congenital dislocation and dysplasia of the hip[END_REF] stage I V.

A computer-assisted design program (Contouring, Symbios, Yverdon, Switzerland), was performed on CT-scan images for computer-assisted analysis. [START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF][START_REF] Sariali | Three-dimensional analysis of the proximal anterior femoral flare and torsion. Anatomic bases for metaphyseally fixed short stems design[END_REF] The first step was to determine the axis of the posterior bicondylar plane, which was automatically reported on all images of the single series. This allowed for the measurement of the angle formed with the metaphyseal axis, represented by the long axis of an oval shape matching the inner contour of the bone at each selected level of the femur and representing the proximal femoral torsion angle (Figure 1).

After determining the midpoint of the lesser trochanter (LT) as a reference point, the levels of CT-scan images to be analysed were selected at 5 mm between 20 mm below LT and 20 mm above LT, and 1 slice at 30 mm above LT. 1 slice at the basis of the femoral neck and 1 slice at the level of the centre of the femoral head were used to measure neck torsion according to the technique described by Murphy et al. [START_REF] Murphy | Femoral anteversion[END_REF] The lowest level of measure at 20 mm below LT was chosen because most cases showed an almost round shape of the femoral canal beneath, thus not allowing a precise measure of proximal femoral torsion.

25 randomly selected cases were blinded and analysed by 2 experienced observers, twice, and 2 days apart, to verify the accuracy of the measures.

All results were recorded as angular values (degrees) with mean and extreme values as well as ± 1 SD (standard deviation) indicated.

Statistical analysis was done with Stata (StataCorp LLC, College Station, TX, USA), using the 2-sample t-test with equal variances. Significance was defined as p-value <0.05.

This study protocol was accepted by the ethical committee of the Faculty of Biology and Medicine of Lausanne University (Switzerland).

Results

A group of 178 files representing 100 men and 78 women (93 right and 85 left hips) with a mean age of 59 ± 14 years (23-84 years) were analysed. We found no significant differences between the groups in means of gender, side or age.

The 25 files which were analysed twice showed a SD of torsion angle of ±2°, not considered as significant (p > 0.05), and in accordance with similar analyses. [START_REF] Sariali | Three-dimensional hip anatomy in osteoarthritis. Analysis of the femoral offset[END_REF] Intra-and inter-observer correlation ranging from 0.95 to 0.91 was considered excellent.

The mean foot progression angle showed no statistically significant difference across all groups, meaning that no torsional leg deformity was present below the knee.

Figure 2 shows the results of mean proximal femoral torsion with maxima and minima in the different groups of pathologies, including the neck angle. Mean femoral torsion in primary OA decreased from 46.1 ± 12° at -20 mm to 19 ± 7.8° at +20 mm from LT, with important individual variations at each level. Mean anterior femoral neck torsion was 10.8 ± 8.5°, varying between 3.1° and 19.5°.

This group was divided into 2 subgroups, 1 with 56 patients <65 years and the other with 44 patients >65 years. There was no significant difference in proximal femoral torsion between the 2 groups at each measured level (p > 0.05), meaning that torsion does not vary with the age of the patient.

In the femoral head osteonecrosis group, proximal femoral torsion decreased from 46.8 ± 17.5° at -20 mm to 19.1 ± 11.7° at +20 mm and in the rheumatoid arthritis group from 50.3 ± 14.4° at -20 mm to 18.6 ± 7.2° at +20 mm.

The epiphysiolysis capitis femoris group had a proximal femoral torsion decreasing from 36.6 ± 16.8° at -20 mm to 23.6 ± 9.9° at +20 mm. In contrast with all other groups, we found an important difference of anterior torsion between the level of the femoral neck basis (mean 16.3°) and the femoral head level (mean 7.9°) explained by the epiphyseal origin of this disease.

No statistically significant difference was found between the proximal femoral torsion values of the reference group of primary OA and the groups of osteonecrosis, rheumatoid arthritis and epiphysiolysis capitis femoris (p > 0.05).

When comparing the proximal femoral torsion at the level of neck osteotomy to the ideal anterior torsion of a prosthetic femoral neck of 15°, we found that 50% of the cases were within a range of ±5° and 50% were outside this range for all pathologies. Proximal femoral torsion of the dysplasia group decreased from 58.4 ± 15.8° at -20 mm to 25.2 ± 12° at +20 mm. The mean anterior torsion of the femoral neck in this group was the same as in the reference group, with a mean value of 14.4 ± 11°.

Finally, the group of congenital dislocation of the hip had a proximal femoral torsion decreasing from 63.1 ± 17.2° at -20 mm to 34.4 ± 18.3° at +20 mm. In this group, anterior neck torsion was significantly higher than in the reference group, with a mean value of 24.8 ± 20.1° (p < 0.05).

The proximal femoral torsion in the groups of dysplasia and congenital dislocation was significantly different at each level compared to the reference group especially for congenital dislocation (p = 0.02 compared to p = 0.002).

When we compared proximal femoral torsion at the level of osteotomy of the femoral neck to an ideal anterior torsion of 15° of a prosthetic femoral neck, we only found 20% of the dysplasia and 5% of the congenital dislocation cases within a range of ±5°.

Discussion

The anatomy of the proximal femur is more sophisticated than often reported and outlined by only a few previous articles. [START_REF] Noble | The anatomic basis of femoral component design[END_REF][START_REF] Rubin | The morphology of the proximal femur. A three-dimensional radiographic analysis[END_REF][START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF] The anatomic variations of the femur have been described in terms of length and anteversion of the femoral neck, as well as inclination (CCD angle). Intramedullary anatomy has shown large individual variations of shape, in the frontal and the sagittal planes. [START_REF] Noble | The anatomic basis of femoral component design[END_REF][START_REF] Rubin | The morphology of the proximal femur. A three-dimensional radiographic analysis[END_REF][START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF] Proximal femoral torsion, defined as the torsion of the upper femur useful for THA, was first reported by Husmann et al. [START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF] They only reported a value of torsion at the level of osteotomy of the femoral neck, but they did not take into account the progression of the torsion along the metaphysis. The present study has shown that torsion was located in the whole proximal femur, beginning at approximately 20 mm under the midpoint of the lesser trochanter. Our preliminary hypothesis that proximal femoral torsion changes progressively along with the proximal femoral metaphysis was verified in all studied groups of different hip diseases leading to THA.

No difference in terms of proximal femoral torsion was found in our study between younger and older patients in the group of primary OA. Noble et al. [START_REF] Noble | The effect of aging on the shape of the proximal femur[END_REF] also found no difference in male femora, but found significant differences in female femora, which may be explained by osteoporosis decreasing bone-density, possibly influencing the frontal and sagittal shape of the femur, but not affecting the horizontal plane, namely the proximal femoral torsion. In their study, Husmann et al. [START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF] found no significant difference between the mean value of proximal femoral torsion (24.0° ± 11.8°) and the mean value of anterior torsion of the neck (24.7° ± 8.7°) in patients with primary OA of the hip, despite broad variation of proximal femoral torsion between 8.3° and 44.9°. In our study, the mean value of anterior torsion of the neck in primary OA was 10.3° ± 8.5° and significantly different (p < 0.05) compared to the proximal femoral torsion at 15 millimetres above the midpoint of the lesser trochanter (approximate level of osteotomy during THA). The same finding was true for all groups of diagnoses.

No significant difference of proximal femoral torsion could be found at any level between the reference group of primary OA and the groups of aseptic osteonecrosis of the femoral head, rheumatoid arthritis and epiphysiolysis capitis femoris. This finding is explained by the epiphyseal or articular origin of the disease, which does not affect femoral metaphysis. However, there were important individual differences of torsion at each level within all these groups. This finding is important because it shows that a femur not only differs from 1 patient to the other in the frontal and the sagittal planes as described by Noble et al. [START_REF] Noble | The anatomic basis of femoral component design[END_REF] and Husmann et al. [START_REF] Husmann | Three-dimensional morphology of the proximal femur[END_REF] but also in the horizontal plane. This adds 1 component to the 3-dimensional anatomy of the femur which should be taken into account when planning or designing femoral stems for THA.

Significant differences of proximal femoral torsion were found at each level between the groups of dysplasia or congenital dislocation of the hip and the reference group of primary osteoarthritis. This difference is explained by the congenital origin of the disease as shown by Tönnis and Heinecke [START_REF] Tönnis | Decreased acetabular anteversion and femur neck antetorsion cause pain and arthrosis. 2: Etiology, diagnosis and therapy[END_REF] and Wilkinson, [START_REF] Wilkinson | Experimental bone growth[END_REF] attributed to the different rotational postures of the foetus, thus affecting the rotation of the whole proximal femur during its development.

The knowledge of the decrease of proximal femoral torsion, from metaphysis to the neck, is very important for the design of a cementless femoral stem. Indeed, the literature shows that the long-term results of cementless hip stems are related to an optimal proximal fit-and-fill to achieve primary and secondary stability with good bone ingrowth. [START_REF] Callaghan | The clinical results and basic science of total hip arthroplasty with porous-coated prostheses[END_REF][START_REF] Rubin | Frictional interface micromotions and anisotropic stress distribution in a femoral total hip component[END_REF][START_REF] Hua | Relative motion of hip stems under load. An in vitro study of symmetrical, asymmetrical, and custom asymmetrical designs[END_REF][START_REF] Robertson | Improving the fit of press-fit hip stems[END_REF][START_REF] Walker | Design and fabrication of cementless hip stems[END_REF] With a straight femoral stem design it is difficult to optimally fit-and-fill the metaphysis because of the variation of proximal femoral torsion depending on the metaphyseal level (Figure 3). Furthermore, this design may induce some peak stresses at the edges of the prosthesis where all loads are transmitted to the bone. [START_REF] Laine | The effects of cementless femoral stem shape and proximal surface texture on 'fit-and-fill' characteristics and on bone remodeling[END_REF] So-called "anatomic" designs have been developed to allow better adaptation of the prosthesis to the intramedullary anatomy. These designs are optimally adapted to cases respecting mean anatomy at each level in the frontal and sagittal planes. As shown in our study, proximal femoral torsion varies between patients within a similar diagnostic group and varies even more in groups of hip dysplasia and congenital dislocation of the hip. Therefore, a straight stem design may achieve optimal fit-and-fill in only 50% of patients, taking into account the optimal anterior torsion of the femoral neck of 15° ± 5° to achieve a stable prosthesis. [START_REF] Herrlin | Femoral anteversion and restricted range of motion in total hip prostheses[END_REF][START_REF] Herrlin | Position, orientation and component interaction in dislocation of the total hip prosthesis[END_REF] Without correction of the orientation of the prosthetic neck, a so-called anatomic stem will achieve an anterior torsion between 10° and 20° in only 50% of the cases in the groups of primary osteoarthritis, aseptic osteonecrosis, rheumatoid arthritis and epiphysiolysis capitis femoris. In the other 50%, if the femoral stem rotates during the implantation to adapt itself to the proximal femoral torsion, this can induce a malrotation of the leg leading to gait problems if the foot is turned inside or outside, or more importantly, an instability of the hip which can lead to dislocation (Figure 4). In cases Figure 3. If a rectangular stem is not adapted to proximal femoral torsion (middle green line without arrow) but implanted so as to achieve the desired anterior torsion of the femoral neck (lower purple arrowed line), the edges of the implant will induce peak loads at the femoral cortex. Furthermore, the size of the implant must be decreased to adapt to the transverse shape of the medullary canal. of dysplasia and congenital dislocation of the hip, the number of patients achieving an optimal fitting with an anterior torsion of the neck between 10° and 20° would be decreased to 20% and 5% respectively, which correlates with the findings of Sugano et al. [START_REF] Sugano | The morphology of the femur in developmental dysplasia of the hip[END_REF] Several previous articles have described techniques for correcting torsion of the femur in regards to femoral neck anterior torsion, [START_REF] Becker | Double-chevron subtrochanteric shortening derotational femoral osteotomy combined with total hip arthroplasty for the treatment of complete congenital dislocation of the hip in the adult. Preliminary report and description of a new surgical technique[END_REF][START_REF] Yasgur | Subtrochanteric femoral shortening osteotomy in total hip arthroplasty for high-riding developmental dislocation of the hip[END_REF][START_REF] Zadeh | Uncemented total hip arthroplasty with subtrochanteric derotational osteotomy for severe femoral anteversion[END_REF][START_REF] Huo | Oblique femoral osteotomy in cementless total hip arthroplasty. Prospective consecutive series with a 3-year minimum follow-up period[END_REF] such as cementing of the femoral stem, the cement acting as "fit-and-fill" between prosthesis and bone. Other techniques including prosthesis undersizing or subtrochanteric femoral osteotomy were also described with poor mid-to long-term results and high morbidity. [START_REF] Hartofilakidis | Treatment of high dislocation of the hip in adults with total hip arthroplasty. Operative technique and long-term clinical results[END_REF][START_REF] Huo | Cementless total hip replacement in patients with developmental dysplasia of the hip[END_REF] Another alternative is a conical prosthesis, influenced by neither proximal femoral torsion nor anterior torsion of the neck. However, this does not allow an optimal fit-and-fill of the metaphysis, even if medium term results in dysplastic patients are satisfactory. [START_REF] Cherubino | Conus" uncemented stem in developmental hip dysplasia[END_REF] Therefore, if a cementless hip prosthesis is to be used, we believe that the better way to simultaneously obtain an optimal proximal fit-and-fill adapted to each proximal femoral torsion and to achieve an ideal anterior torsion of 15° ± 5° of the femoral neck is to plan every single patient THA in 3 dimensions based on CT scans. A precise analysis of the torsional anatomy of each individual femur allows the determination of the single best femoral stem for each case, both intramedullary and extramedullary. [START_REF] Sariali | Accuracy of reconstruction of the hip using computerised three-dimensional pre-operative planning and a cementless modular neck[END_REF][START_REF] Wegrzyn | The tridimensional geometry of the proximal femur should determine the design of cementless femoral stem in total hip arthroplasty[END_REF] It will also indicate those cases where the torsion of the proximal femoral metaphysis does not allow us to simultaneously achieve both a perfect fit-and-fill and anterior torsion (Figure 5). This is the case in dysplasia (80%) and congenital dislocation of the hip (95%), but also in approximately 50% of the patients with other diagnoses. A modular stem could be an alternative, adapting the orientation of the stem to the metaphyseal torsion, and correcting the anterior torsion to its optimal angulation. However, this adaptation is often limited to a few degrees and the modularity may also introduce potential metal ion problems. In cases where no optimal adaptation between proximal femoral torsion and anterior torsion can be found, a custom femoral stem may be an attractive alternative in order to offer the patient a femoral stem with the most optimal correction and bony adaptation. [START_REF] Argenson | Anatomy of the dysplastic hip and consequences for total hip arthroplasty[END_REF][START_REF] Argenson | Threedimensional anatomy of the hip in osteoarthritis after developmental dysplasia[END_REF][START_REF] Wettstein | Threedimensional computed cementless custom femoral stems in young patients: midterm followup[END_REF][START_REF] Tsiampas | Custom-made femoral implants in total hip arthroplasty due to congenital disease of the hip: a review[END_REF] The limitations of this study include the small number of patients analysed and the absence of measurements on the whole femur. However, the main portion affected by the torsion of the femur at the time of THA is located on the proximal metaphysis, which is covered in our study. We believe this series represents a significant contribution to anatomic knowledge of the proximal femoral metaphyseal torsion in various aetiologies leading to THA.

The knowledge of proximal femoral torsion, its values and variations, are of tremendous importance to the surgeon who believes in the principles of fit-and-fill. [START_REF] Hua | Relative motion of hip stems under load. An in vitro study of symmetrical, asymmetrical, and custom asymmetrical designs[END_REF][START_REF] Robertson | Improving the fit of press-fit hip stems[END_REF][START_REF] Walker | Design and fabrication of cementless hip stems[END_REF] As proximal femoral torsion can only be assessed by CT scan, we suggest that, ideally, each patient with a high life expectancy or suspected dysplasia and who is scheduled for a cementless THA, may benefit from 3-dimensionnal computer assisted planning with the goal of a perfect intramedullary adaptation and anatomic extramedullary reconstruction. We also believe the results of the present study could help in the design of cementless stems in THA for various anatomies. This can be considered as one further step in optimising the success rate of cementless total hip arthroplasty. In patients with an increased proximal femoral torsion (middle green line without arrow), if the goal is to achieve optimal adaptation of the stem to the intramedullary shape of the bone without compromising femoral neck anterior torsion (lower purple arrowed line), a modular or custom stem should be considered.

Figure 1 .

 1 Figure 1. Measurement of anterior torsion: (a) at the level of the femoral neck (N) and proximal femoral torsion and (b) at the level of the proximal femoral metaphysis (F). The oval shape helps measuring the axis of the metaphysis.

Figure 2 .

 2 Figure 2. Mean, minimal and maximal femoral torsion at several levels in all diagnostic groups.

Figure 4 .

 4 Figure 4. A stem which is implanted according to an increased proximal femoral torsion (middle green line without arrow) will induce an increased femoral neck anterior torsion (lower purple arrowed line).

Figure 5 .

 5 Figure 5.In patients with an increased proximal femoral torsion (middle green line without arrow), if the goal is to achieve optimal adaptation of the stem to the intramedullary shape of the bone without compromising femoral neck anterior torsion (lower purple arrowed line), a modular or custom stem should be considered.

Table 1 .

 1 Distribution of gender, side and age according to diagnosis.

	Diagnosis	Gender	Side	Mean age (years)
	Primary osteoarthritis	55 men, 45 women	52 right, 48 left	60.5 ± 12.5 (30-84)
	Dysplasia	6 men, 14 women	12 right, 8 left	53.9 ± 17.2 (23-78)
	Congenital dislocation of the hip	6 men, 14 women	11 right, 9 left	57.8 ± 13.9 (27-78)
	Aseptic osteonecrosis of femoral head	15 men, 5 women	9 right, 11 left	55.4 ± 16.1 (32-80)
	Rheumatoid polyarthritis

men, 4 women

right, 3 left 63.5 ± 11.3 (47-77) Epiphysiolysis capitis femoris 6 men, 3 women 5 right, 4 left 55.3 ± 12.7 (35-70)

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.