Boris P Hejblum
email: boris.hejblum@u-bordeaux.fr

Chariff Alkhassim

Raphael Gottardo

Rodolphe Thiébaut

SEQUENTIAL DIRICHLET PROCESS MIXTURES OF MULTIVARIATE SKEW T -DISTRIBUTIONS FOR MODEL-BASED CLUSTERING OF FLOW CYTOMETRY DATA

Keywords: Primary 62H30, 62P10; secondary 62L12 Automatic gating, Bayesian Nonparametrics, Dirichlet process, Flow cytometry, HIV, Mixture model, Skew t-distribution

proposed nonparametric Bayesian mixture models of Gaussian distributions, that directly estimate this number of cell populations. All these methods, except those of Lin et al. (2013), of Cron et al. (2013) and of Dundar et al. (2014), were evaluated by Aghaeepour et al. (2013).

) who proposed to use hierarchical Dirichlet process mixture (DPM) of Gaussian distribution models to analyze multiple samples simultaneously. Yet in the case of repeated measurements of flow cytometry data, it can be useful to perform a sequential analysis as the samples are acquired (samples are often

Flow cytometry is a high-throughput technology used to quantify multiple surface and intracellular markers at the level of a single cell. This enables us to identify cell sub-types, and to determine their relative proportions. Improvements of this technology allow us to describe millions of individual cells from a blood sample using multiple markers. This results in high-dimensional datasets, whose manual analysis is highly time-consuming and poorly reproducible. While several methods have been developed to perform automatic recognition of cell populations, most of them treat and analyze each sample independently. However, in practice, individual samples are rarely independent, especially in longitudinal studies. Here we analyze new longitudinal flow-cytometry data from the DALIA-1 trial which evaluates a therapeutic vaccine against HIV, by proposing a new Bayesian nonparametric approach with Dirichlet process mixture (DPM) of multivariate skew t-distributions to perform model based clustering of flow-cytometry data. DPM models directly estimate the number of cell populations from the data, avoiding model selection issues, and skew t-distributions provides robustness to outliers and non-elliptical shape of cell populations. To accommodate repeated measurements, we propose a sequential strategy relying on a parametric approximation of the posterior. We illustrate the good performance of our method on simulated data and on an experimental benchmark dataset. This sequential strategy outperforms all other methods evaluated on the benchmark dataset, and leads to improved performance on the DALIA-1 data.

1. Introduction. Flow cytometry is a high-throughput technology used to quantify multiple surface and intracellular markers at the level of single cell. More specifically, cells are stained with multiple fluorescently-conjugated monoclonal antibodies directed to cell surface receptors (such as CD4) or intracellular markers (such as cytokines) to determine the type of cell, their differentiation, and their functionality. With the improvement of this technology leading currently to the measurement of up to 18 markers at the same time (using 18 colors for flow cytometry), multi-parametric description of millions of individual cells can be generated.

Analysis of such data is generally performed manually. This results in analyses that are: i) poorly reproducible [START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF], ii) expensive (highly time-consuming) and consequently iii) only focused on specific cell populations (i.e. specific combination of markers), ignoring other cell populations. There has been an effort in the recent years to offer automated solutions to overcome these limitations [START_REF] Lo | Automated gating of flow cytometry data via robust model-based clustering[END_REF][START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF][START_REF] Gondois-Rey | Multi-parametric cytometry from a complex cellular sample: Improvements and limits of manual versus computational-based interactive analyses[END_REF]. Quite a lot of different methodological approaches have been proposed to perform automatic recognition of cell populations from flow cytometry data. Clustering meth-collected across several time points in a population of patients). In such a case, one would want to use previously acquired samples as prior information in the analysis of a new sample. In this paper, the proposed approach includes a strategy of sequential approximations of the posterior distribution for multiple data samples, presented in Section 3.2. Our approach offers three advantages: i) it quantifies the uncertainty of the posterior clustering, ii) it can make use of prior knowledge to inform on the structure of the data, potentially building up on previous analyses, and iii) it allows the analysis of multiple samples without requiring to process all the data at once, alleviating both the computational burden and the necessity for all data to be readily available before any analysis can be performed.

The automatic recognition of cell populations from flow cytometry data is a difficult task which can be seen as an unsupervised clustering problem [START_REF] Lo | Automated gating of flow cytometry data via robust model-based clustering[END_REF]. It is characterized by two big challenges. First, the total number of cell populations to identify is unknown. Second, the empirical distributions of the populations are heavily skewed, even when optimal transformation of the data is applied [START_REF] Lo | Automated gating of flow cytometry data via robust model-based clustering[END_REF][START_REF] Pyne | Automated high-dimensional flow cytometric data analysis[END_REF][START_REF] Lo | Flexible mixture modeling via the multivariate t distribution with the Box-Cox transformation: An alternative to the skew-t distribution[END_REF], and the data generally present many outliers. To address all these points together, our approach considers a Bayesian nonparametric model-based approach, where the flow cytometry data are assumed to be drawn from a DPM of multivariate skewt distributions. First, this approach enables the number of cell populations to be inferred from the data and avoids the challenging problem of model selection. Second, it has been demonstrated that the Gaussian assumption for the parametric shape of a cell population fits poorly flow cytometry data [START_REF] Mosmann | SWIFT-scalable clustering for automated identification of rare cell populations in large, high-dimensional flow cytometry datasets, Part 2: Biological evaluation[END_REF]. Indeed, even after state-of-the-art transformation of raw cytometry data, such as the biexponential transformation [START_REF] Finak | Optimizing transformations for automated, high throughput analysis of flow cytometry data[END_REF], cell population distributions are typically skewed. [START_REF] Pyne | Automated high-dimensional flow cytometric data analysis[END_REF] have showed the advantages of the skew t-distribution [START_REF] Azzalini | Distributions generated by perturbation of symmetry with emphasis on a multivariate skew t-distribution[END_REF] for modeling cell populations in flow cytometry data. Numerous parameterizations have been proposed for the multivariate skew t-distribution [START_REF] Lee | On mixtures of skew normal and skew tdistributions[END_REF][START_REF] Murray | Mixtures of skew-t factor analyzer[END_REF][START_REF] Azzalini | On nomenclature for, and the relative merits of, two formulations of skew distributions[END_REF][START_REF] Mclachlan | Comment on On nomenclature, and the relative merits of two formulations of skew distributions by A[END_REF], most notably the restricted and the unrestricted multivariate skew t-distributions (denoted rMST and uMST respectively) which are generalizations of the skew normal distribution [START_REF] Azzalini | The multivariate skew-normal distribution[END_REF] with a heavier tail (making it more robust to outliers). [START_REF] Lee | Finite mixtures of canonical fundamental skew t-distributions: The unification of the restricted and unrestricted skew t-mixture models[END_REF] recently proposed the canonical fundamental skew t-distribution (CFUST) as a generalization that encompasses both the rMST and the uMST. To avoid identifiability issues associated with the uMST and the CFUST [START_REF] Lee | Finite mixtures of canonical fundamental skew t-distributions: The unification of the restricted and unrestricted skew t-mixture models[END_REF], in this work we will adopt the rMST formulation, and in the remainder of this article we will be referring to the rMST formulation when mentioning the skew t-distribution. Frühwirth-Schnatter and [START_REF] Frühwirth-Schnatter | Bayesian inference for finite mixtures of univariate and multivariate skew-normal and skew-t distributions[END_REF] proposed a finite mixture model of rMST. We extend this model to the infinite mixture case in a Bayesian nonparametric framework. Of interest, quantifying the uncertainty around the estimated partition is straightforward in this Bayesian paradigm, from the posterior distribution of the partition. While a skewed distribution could be fitted either by a skew t or by a mixture of Gaussians, using the latter requires to estimate separately the overall number of clusters and the skewness. On the contrary, our proposed approach jointly estimates the two thus taking into account the uncertainty associated with both. Furthermore, the use of a Bayesian framework enables the use of informative priors. In the case of repeated measurements for instance, we propose to sequentially estimate the posterior partition of flow cytometry using posterior information from time point t as prior information for time point t `1.

The proposed method is evaluated on simulated data and on a benchmark clinical dataset from [START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF], and is applied to analyze an original experimental longitudinal dataset from a phase I HIV clinical trial DALIA-1 (featuring depending time-course data where the sequential approach is of particular interest). The method is implemented in the R package NPflow, available on the CRAN at https://CRAN.R-project.org/ package=NPflow.

Statistical Model.

2.1. Motivation and problem set-up.

2.1.1. Motivating example. Our motivating example for developing a sequential model-based clustering approach for longitudinal flow cytometry data comes from the DALIA-1 trial. DALIA-1 is a phase I trial for a therapeutic vaccine candidate against HIV [START_REF] Lévy | Dendritic cell-based therapeutic vaccine elicits polyfunctional HIV-specific T-cell immunity associated with control of viral load[END_REF]. This vaccine candidate was based on ex-vivo generated interferon-α dendritic cells loaded with HIV-1 lipo-peptides, and activated with lipopolysaccharide. The primary objectives of this trial were to evaluate the safety of the vaccine strategy and to evaluate the immune response. As part of this trial, 12 HIV positive patients had their cellular populations quantified repeatedly by flow-cytometry, generating an important amount of data, for which comprehensive manual gating would take several months. Hence we aimed at developing an automatic gating approach suitable for longitudinal measurements.

2.1.2. Problem set-up and notations. We first consider a single sample per subject, i.e. one data matrix where each row represents a cell and each column contains the fluorescence intensities for one marker measured by the flow cytometer. The case of the sequential estimation of multiple datasets will be addressed in Section 3.2. We let y c P R d denote the data, c " 1, . . . , C corresponding to the vector of fluorescence intensities measured for the cell c. Typically, the observations y c have been transformed (to help visualization and gating) from the raw measurements of fluorescence through a biexponential or Box-Cox transformation [START_REF] Finak | Optimizing transformations for automated, high throughput analysis of flow cytometry data[END_REF]. We assume that these observations are independent and identically distributed (i.i.d.) from some unknown distribution F :

(2.1)

y c |G i.i.d. " F for c " 1 . . . , C
where F is a mixture of distributions:

(2.2)

F pyq " ż Θ f θ pyqGpdθq
with f θ pyq is a known probability density function, parameterized by θ P Θ a set of parameters, and defining the shape of a cluster. G is the unknown mixing distribution, which carries the weights and locations of the mixture components. In a parametric approach, G " ř K k"1 π k δ θ k where π k is the weight of the k th mixture component. Maximum likelihood or Bayesian estimates of F can be derived for such models [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. In a nonparametric perspective (where the number of clusters is unknown) G is written as an infinite sum of atoms: G " ř `8 k"1 π k δ θ k . The Dirichlet process is a conjugate prior for the infinite atomic discrete distribution, which makes it very useful for unsupervised clustering approaches.

2.2. Dirichlet process mixture of skew t-distributions.

2.2.1. Dirichlet process mixture. We assume that the random mixing distribution G is drawn from a Dirichlet process [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF]:

(2.3) G " DPpα, G 0 q
where DPpα, G 0 q denotes the Dirichlet process of concentration parameter α ą 0 and base probability distribution G 0 . A draw G " DPpα, G 0 q is almost surely discrete [START_REF] Sethuraman | A constructive definition of Dirichlet priors[END_REF] and gives a nonparametric mixing distribution G "

ř `8 k"1 π k δ θ k with θ k i.i.d.
" G 0 , and π " pπ k q k"1,2,... drawn from a so-called "stick-breaking" distribution, written as the Griffiths-Engen-McCloskey (GEM) distribution [START_REF] Pitman | Combinatorial Stochastic Processes[END_REF]. The model defined by Equations (2.1), (2.2) and (2.3) yields the following hierarchical model known as a Dirichlet process mixture model [START_REF] Lo | On a class of Bayesian nonparametric estimates: I. Density estimates[END_REF][START_REF] Escobar | Bayesian Density Estimation and Inference Using Mixtures[END_REF][START_REF] Teh | Dirichlet Process[END_REF] with a Gamma hyperprior on α: α|a, b " Gammapa, bq (2.4a)

π ˇˇα " GEMpαq (2.4b) θ k ˇˇG 0 " G 0 for k " 1, 2, . . . (2.4c) c ˇˇπ " Multpπq for c " 1, 2, . . . , C (2.4d) y c ˇˇ c , pθ k q " f θ c for c " 1, 2, . . . , C (2.4e)
where c is the latent cluster-allocation for cell c. G 0 tunes the prior information about the cluster locations while α tunes the prior distribution on the overall number of clusters K that will be uncovered within C cells. In particular we have ErK|Cs "

ř C´1 c"0 α α`c .
2.2.2. The multivariate skew t-distribution. Frühwirth-Schnatter and Pyne (2010) rely on [START_REF] Azzalini | The multivariate skew-normal distribution[END_REF]'s parametrization of the multivariate skew normal (SN ) to propose a truncated normal random-effects model representation of this distribution: Y " ξ `ψZ `ε with Z " N r0,`8r p0, 1q and ε " N p0, Σq. If X " SN p0, Ω, ηq and W " Gammap ν 2 , ν 2 q, [START_REF] Azzalini | Distributions generated by perturbation of symmetry with emphasis on a multivariate skew t-distribution[END_REF] show that Y " ξ `1 ? W X then follows a multivariate skew t-distribution: Y " ST pξ, Ω, η, νq. Following Frühwirth-Schnatter and Pyne (2010), we write the density of a multivariate skew t-distribution as:

f ST py; ξ, Ω, η, νq "2f T py; ξ, Ω, νqT ν`d ˜η1 ω ´1py ´ξq d ν `d ν `Qy (2.5)
with ω " a DiagpΩq, Q y " py ´ξq 1 Ω ´1py ´ξq, f T the multivariate Student t-distribution density, and T ν`d the cumulative distribution function of the standard univariate Student's t-distribution with ν `d degrees of freedom. This parametrization of the skew t is referred as the restricted multivariate skew t distribution by [START_REF] Lee | On mixtures of skew normal and skew tdistributions[END_REF], and it admits the following random-effect model representation:

(2.6) Y " ξ `ψ Z ? W ` ? W 2.2.3.
Dirichlet process mixture of multivariate skew t-distributions. Combining model (2.4) with a random-effects model representation (2.6) of the skew t-distribution, we propose the following model:

α|a, b " Gammapa, bq (2.7a) π ˇˇα " GEMpαq (2.7b) for k " 1, 2, . . . ξ k , ψ k , Σ k , ν k " G 0 (2.7c) for c " 1, 2, . . . , C c ˇˇπ " Multpπq (2.7d) γ c ˇˇ c , tν k u " Gamma ´ν c 2 , ν c 2 (2.7e) s c ˇˇγ c " N r0,`8r ˆ0, 1 γ c
(2.7f)

y c ˇˇ c , γ c , s c , pξ k , ψ k , Σ k q " N ˆξ c `ψ c s c , 1 γ c Σ c (2.7g)
where G 0 is the product of a structured normal-inverse-Wishart (sN iW ) and of a prior on ν: G 0 " sN iW pξ 0 , ψ 0 , B 0 , Λ 0 , λ 0 qP 0,ν .

2.3. Discussion on the model assumptions. In model (2.7), the base distribution parameter G 0 conveys the prior information on the cluster parametric shape. For the parameters ξ k , ψ k and Σ k , we have conditional conjugacy with the random-effects model representation using joint priors taking the form of a structured normal-inverse-Wishart distribution. See Online Supplement A for details [START_REF] Hejblum | Supplement to "Sequential Dirichlet process mixtures of multivariate skew tdistributions for model-based clustering of flow cytometry data[END_REF]. Frühwirth-Schnatter and [START_REF] Frühwirth-Schnatter | Bayesian inference for finite mixtures of univariate and multivariate skew-normal and skew-t distributions[END_REF] pointed out that the prior on Σ k can have a big impact on the posterior number of clusters. Indeed, setting the scale of the prior on Σ k too small will result in an inflated number of clusters in the posterior, whereas too large values tend to cluster all the observations together. Adding a Wishart hyperprior on Σ k , that carries on conjugacy with the inverse-Wishart, enables us to reduce this impact of the prior (Frühwirth-Schnatter and [START_REF] Frühwirth-Schnatter | Bayesian inference for finite mixtures of univariate and multivariate skew-normal and skew-t distributions[END_REF][START_REF] Huang | Simple Marginally Noninformative Prior Distributions for Covariance Matrices[END_REF]. Assuming prior independence between each ν k and also from the three parameters mentioned above, we can use any of the three priors proposed in [START_REF] Juárez | Model-Based Clustering of Non-Gaussian Panel Data Based on Skew-t Distributions[END_REF] for instance (such as an objective Jeffrey's prior, see Online Supplement A [START_REF] Hejblum | Supplement to "Sequential Dirichlet process mixtures of multivariate skew tdistributions for model-based clustering of flow cytometry data[END_REF]).

Estimation.

3.1. Posterior Estimation via Gibbs sampling. For making inference on the model (2.7), MCMC methods can be used to sample the partition t 1:C u and the corresponding cluster parameters tθ k u " ttξ k u, tψ k u, tΣ k u, tν k uu from the marginal posterior distribution. Extending results from Frühwirth-Schnatter and Pyne (2010) and [START_REF] Caron | Bayesian nonparametric Placket-tLuce models for the analysis of preferences for college degree programmes[END_REF], it is possible to implement an efficient and valid partially collapsed Gibbs sampler with a Metropolis-Hastings step [START_REF] Van Dyk | Partially Collapsed Gibbs Samplers[END_REF][START_REF] Van Dyk | Metropolis-Hastings within Partially Collapsed Gibbs Samplers[END_REF]. The use of slice sampling [START_REF] Neal | Slice sampling[END_REF][START_REF] Kalli | Slice sampling mixture models[END_REF] enables the straightforward parallelization of the latent allocation sampling (thanks to conditional conjugacy) in such an MCMC algorithm (even in the skew normal and skew t cases), which can lead to substantial computation speed up when the number of observations C (cells) per sample increases. Each iteration of our Gibbs sampler proceeds in the following order (details are provided in Online Supplement A [START_REF] Hejblum | Supplement to "Sequential Dirichlet process mixtures of multivariate skew tdistributions for model-based clustering of flow cytometry data[END_REF]):

1. Update the concentration parameter α given the previous partition t 1:C u using the data augmentation technique from [START_REF] Escobar | Bayesian Density Estimation and Inference Using Mixtures[END_REF]. 3.2. Sequential Posterior Approximation. In flow cytometry experiments, it is common to actually have multiple datasets y piq (with i " 1, . . . , I) corresponding to multiple individuals or repeated measurements of the same individual. In such cases, it is of interest to use previous time points or previous samples results as prior information, in order to leverage all the information available to estimate the mixture. However, incorporating prior information into Dirichlet process mixture models is not straightforward [START_REF] Kessler | Marginally specified priors for non-parametric bayesian estimation[END_REF]. Here we propose to use the posterior MCMC draws obtained from previous dataset y piq as prior information to analyze the next dataset y pi`1q . To do so, first we consider the hierarchical model using all observations from both y piq and y pi`1q at once :

α " Gammapa, bq (3.1a) G|α " DP pα, G 0 q (3.1b) y piq , y pi`1q |G i.i.d. " ż Θ f θ p¨qdGpθq (3.1c)
We are interested in estimating ppG|y piq , y pi`1q q9ppG|y piq qppy pi`1q |Gq. The idea is to first approximate ppG|y piq q by a Dirichlet process through MCMC draws from the model described in 2.1:

(3.2) ppG|y piq q » ż DP pG; α, G 1 qGammapα; a 1 , b 1 qdα
where G 1 , a 1 , b 1 are parameters to be estimated from the MCMC approximation of the true posterior: i) p a 1 and p b 1 can be taken as MLE estimates from the MCMC samples α pjq ; ii) x G 1 is a parametric approximation of the posterior mixing distribution G 1 (the true posterior is not suitable for being directly plugged in as a base distribution parameter of another DP as it is nonparametric). In the case of a skew t-distribution mixture model, we approximate G 1 with the following joint distribution: G 1 » psN iW, P 0,ν q where P 0,ν is the chosen prior for the skew t-distribution degrees of freedom. To estimate G 1 , we estimate the maximum a posteriori (MAP) from the posterior MCMC samples (see Online Supplement B [START_REF] Hejblum | Supplement to "Sequential Dirichlet process mixtures of multivariate skew tdistributions for model-based clustering of flow cytometry data[END_REF]). Now using this posterior parametric approximation, we have the same hierarchical model as before but conditional on y piq :

α|y piq " Gammap p a 1 , p b 1 q (3.3a) G|α, y piq " DP pα, x G 1 q (3.3b) y pi`1q |G, y piq i.i.d. " ż Θ f θ p¨qdGpθq (3.3c)
Note that under this approximate posterior model, the cluster parameters θ k are i.i.d. according from G 1 . Such an approach can be iterated a number of times, if for instance several time points are observed, iteratively approximating the successive posteriors. This approach allows us to finally account for all the previous information in the mixture model estimation. This model hypothesizes that all the data are originating from the same mixture model and as more data are acquired, the successive posteriors will concentrate: since the overall posterior is our target, this concentration effect is desired. When the data are anticipated to be non stationary, which might be the case in real life applications, we propose a forgetful factor by adding a vague component in the base distribution of the prior.

3.3. Point estimate of the clustering. Getting a representation of the partition posterior distribution is difficult [START_REF] Medvedovic | Bayesian infinite mixture model based clustering of gene expression profiles[END_REF]. One can use the maximum a posteriori, i.e. using the point estimation from the MCMC sample that maximizes the posterior density. However this ignores all the information about the uncertainty around the partition provided by the Bayesian approach.

Another way is to rather consider a co-clustering posterior probability (or similarity) matrix ζ on each pair pc, dq of observations. Such a matrix can be estimated by averaging the co-clustering matrices from all the explored partitions in the posterior MCMC draws:

(3.4) p ζ cd " 1 N N ÿ i"1 δ piq c piq d
where N is the number of MCMC draws from the posterior, i the MCMC iteration, and δ kl " 1 if k " l, 0 otherwise. An optimal partition point estimate t p 1:C u can then be derived in regard of this similarity matrix through stochastic search with the explored partitions in the posterior MCMC draws [START_REF] Dahl | Model-Based Clustering for Expression Data via a Dirichlet Process Mixture Model[END_REF], by using a pairwise coincidence loss function [START_REF] Lau | Bayesian Model-Based Clustering Procedures[END_REF]) such as the one proposed by [START_REF] Binder | Bayesian Cluster Analysis[END_REF][START_REF] Binder | Approximations to Bayesian Clustering Rules[END_REF] which optimizes the Rand index [START_REF] Fritsch | Improved criteria for clustering based on the posterior similarity matrix[END_REF]):

(3.5) t p 1:C u " arg min

t piq 1:C uP ! t p1q 1:C u,...,t pN q 1:C u ) C´1 ÿ c"1 C ÿ d"c`1 2 ´δ piq c piq d ´p ζ cd

¯2

The computational complexity of this approach however is of the order OpN C 2 q due to the necessity of computing all the similarity matrices.

A different optimal partition point estimate t ˜ 1:C u can also be derived using the F-measure as our loss function. The F-measure is widely used as a way to summarize the accordance between 2 methods, one being considered as a reference (gold-standard). It is the harmonic mean of the precision and recall:

(3.6)

F " 2P rRe P r `Re

In order to use the F-measure to evaluate our clustering method, we rely on the definition proposed in the online methods from [START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF].

In this unsupervised clustering setting, the precision P r is the number of cells correctly assigned to a given cluster divided by the total number of cells assigned to that cluster (also called Positive Predictive Value). The recall Re is the number of cells correctly assigned to a given cluster divided by the number of cells that should be assigned to this cluster according to the gold-standard. Since in our problem the labels of the different clusters are exchangeable, the F-measure is computed for each combination of the reference clusters and the predicted clusters. Let G " tg 1 , . . . , g m u be a set of m reference clusters and H " th 1 , . . . , h n u be set of n predicted clusters.

For each combination pair of a reference cluster g q and a predicted cluster h r , the F-measure is computed as follows:

(3.7) P rph r , g q q " |g q X h r | |h r | and Reph r , g q q " |g q X h r | |g q |

(3.8) Fph r , g q q " 2P rpg q , h r qRepg q , h r q P rpg q , h r q `Repg q , h r q

This F-measure is comprised in r0, 1s and the closer it is to 1 the better the agreement is between the predicted cluster and the reference cluster. The total F-measure for a predicted partition H given a gold-standard G is then defined as the weighted sum of the best matched F-measure:

(3.9)

F tot pH, Gq " 1 ř m q"1 |g q | m ÿ q"1 |g q | max rPt1,...,nu
Fph r , g q q This total F-measure is again between 0 and 1, and the closer it is to 1 the better the predicted partition agrees with the gold-standard. The optimal partition point estimate in respects of this F-measure is then obtained with the partition that maximizes its average F-measure over all the other explored partitions in the posterior MCMC draws:

(3.10) t ˜ 1:C u " arg max t piq 1:C uP ! t p1q 1:C u,...,t pN q 1:C u ) 1 N N ÿ j"1 j‰i F tot ´t piq 1:C u, t pjq 1:C u

Note

the F-measure is computed here only between sampled partitions, and a gold-standard partition is unnecessary.

Simulation Study.

4.1. Weakly informative prior. First, to assess the performance of the Dirichlet process mixture of skew t distributions model in a simple clustering case, 100 simulations in 2-dimensions were performed. In each simulation, 10,000 observations were drawn from 5 distinct clusters representing respectively 50%, 29.9%, 15%, 5%, and 0.1% of the data. After 20,000 MCMC iterations (18,000 iterations burnt and a thining of 20 gave 100 partitions sampled from the posterior; the chain was initialized with 30 clusters), we compared the partition point estimate obtained from our approach with the true clustering of the simulated data through the resulting mean F-measure. When the clustering problem is well characterized, i.e. when the true clusters are well separated [START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF], the median F-measure was 0.998. When considering overlapping clusters, i.e. when the true clustering is not entirely recoverable, our approach was nevertheless able to maintain good performance with a median F-measure of 0.895. Figure 1 shows an example of the partition point estimate obtained for one of those simulation runs in both cases, where one can see that NPflow is able to correctly recover the 5 clusters, including the extremely small one of 0.1% (green diamond). Of course, if this extremely small cluster was not well separated from larger clusters in Figure 1B, the data would hardly contain any evidence of its presence and the model would likely not recover it. As a comparison, k-means only reached a F-measure of 0.920 and 0.823 in the well-separated and overlapping scenarios respectively, in spite of having the correct number (5) of true clusters specified.

Comparison between skew-t and Gaussian kernels. Figure 2 illustrates the improvement due to the use of a skew t kernel over a Gaussian kernel in a nonparametric mixture model. Figures 2A and2C both show that the F-measure is significantly better with a skew t kernel than with a Gaussian kernel for well-separated and overlapping settings respectively, while Figures 2B and 2D both show that skew t kernels allow to accurately recover the true number of clusters in a majority of cases for well-separated and overlapping settings respectively. Both criteria also highlight the increased difficulty in the overlapping setting as the clustering problem becomes less well characterized.

Sequential posterior approximation plugged-in as informative prior.

To illustrate how the sequential posterior approximation strategy compares to the standard weakly informative prior setting, we ran simulations where we considered two samples derived from the same infinite mixture model. The first sample is simulated for a time t, and the second sample at t `1. As all observations originate from the exact same distribution, regardless of the Partition point estimate from one of the 100 2-dimensional simulations with weakly informative prior (10,000 observations with 5 clusters representing respectively 50%, 29.9%, 15%, 5% and 0.1% of the data). A: well separated clusters ; B: overlapping clusters sample, the hypothesis of the sequential posterior approximation strategy is satisfied. One of the major gain observed is the time to convergence for the partition. Using an informative prior derived from the sample at time t to estimate the partition of the sample from t `1 makes it more than three time faster to converge according to the Gelman-Rubin statistics.

In further simulations, we also investigated the performance of this sequential posterior approximation strategy. opposed to using the standard weakly informative prior strategy, it shows substantial gains when the amount of information brought by the prior is substantial compared to the amount available from the data at time t `1 alone. As the amount of information available at time t `1 increases, the gain from using this strategy can become less noticeable, as shown using the F-measure in Figure 3. But even when the number of observations available at time t `1 is the same as at time t, the accuracy for rare cell populations is still improved by using an informative prior. This is not necessary visible at the scale of the total F-measure, because it is masked by the larger clusters. However, when computing a limited F-measure, that only takes into account smaller clusters, the use of an informative prior in this sequential strategy seems to always im- prove the clustering accuracy for smaller clusters (see Supplementary Figure S1 in Online Supplement C [START_REF] Hejblum | Supplement to "Sequential Dirichlet process mixtures of multivariate skew tdistributions for model-based clustering of flow cytometry data[END_REF]).

Application to real datasets.

5.1. Benchmark dataset. The Graft versus Host Disease (GvHD) dataset is a public dataset that was first analyzed (manually gated) in [START_REF] Brinkman | High-content flow cytometry and temporal data analysis for defining a cellular signature of graft-versus-host disease. Biology of blood and marrow transplantation : journal of the American Society for Blood and[END_REF], with the objective of identifying a cellular signature that correlates or predicts the Graft versus Host disease. These GvHD data were used as benchmark data in the FlowCAP challenge [START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF]. Flow cytometry data was collected for 12 samples, and original manual gates are being regarded as the true cell clustering (actually a consensus over eight different manual operators). In an attempt to mitigate further the well known reproducibility issues with manual gating [START_REF] Ge | flowPeaks: a fast unsupervised clustering for flow cytometry data via K-means and density peak finding[END_REF][START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF] 1 shows the performance of our proposed approach NPflow on this dataset, compared to the other approaches reviewed by [START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF]. The F-measure is computed for all samples available for a given dataset and the mean over all samples is reported, as well as a bootstrap 95% confidence interval. No algorithm is performing significantly better than NPflow thus placing NPflow among the top methods for automatic gating. Additionally, we compared the use of a skew t kernel by NPflow with the use of a Gaussian kernel (denoted NPflowG in Table 1), the latter reaching a mere 0.61 F-measure on average thus demonstrating the benefit of the skew t distributions for modeling real flow-cytometry data. Thanks to its use of sequential information, the sequential posterior model would ideally improve results by analyzing each individual sample sequentially. The GvHD benchmark data are not longitudinal but, as long as the different samples are similar enough, one can expect an improvement. On this benchmark, our sequential approach NPflow-seq reaches a mean Fmeasure of 0.89 (0.85, 0.94) compared to a value of 0.85 (0.80, 0.90) with the standard NPflow model (Table 1). Of note, the only other approach not analyzing the samples independently, and that relies on a hierachical Dirichlet process Gaussian mixture model (HDPGMM), only reaches a value of 0.35 (0.30, 0.39) [START_REF] Cron | Hierarchical modeling for rare event detection and cell subset alignment across flow cytometry samples[END_REF]; [START_REF] Johnsson | BayesFlow: latent modeling of flow cytometry cell populations[END_REF]. This illustrates that integrating all samples in a simultaneous model does not necessarily yield better results (e.g. if the global model across samples is misspecified or not flexible enough). For the GvHD benchmark dataset, our sequential strategy thus exhibits the highest F-measure compared to standard NPflow and to HDPGMM, and also to competing unsupervised automatic gating methods evaluated in [START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF] that however analyze each sample independently. It is worth noting that since our sequential strategy performs sequential approximations of the posterior and provides intermediate results for each sample, the order in which observations are included can have an impact (especially for the first sample). Here we analyzed the GvHD samples in their original order as provided by their identifiers in FlowCAP-I (from 001 to 012).

5.2.

Original DALIA-1 data: a longitudinal real data study. We applied our method to analyze an original dataset from DALIA-1, a phase I trial evaluating a therapeutic vaccine against HIV [START_REF] Lévy | Dendritic cell-based therapeutic vaccine elicits polyfunctional HIV-specific T-cell immunity associated with control of viral load[END_REF]. For our purpose here, we are interested in the 12 HIV positive patients who had their cellular populations quantified at 18 time points during the trial. More specifically, we focused on two time points (at week 24 and week 26 of the trial) immediately following antiretroviral treatment (HAART) interruption which took place at week 24. Following this interruption, the increase of viral replication is associated with changes in cell populations [START_REF] Thiébaut | Immunological markers after long-term treatment interruption in chronically HIV-1 infected patients with CD4 cell count above 400 x 10(6) cells/l[END_REF][START_REF] Lévy | Effect of intermittent interleukin-2 therapy on CD4+ T-cell counts following antiretroviral cessation in patients with HIV[END_REF]. Here we especially looked at the CD4+ effector T-cells, defined as CD45RA+CD27-among the CD3+CD4+ cells [START_REF] Larbi | From "truly naïve" to "exhausted senescent" T cells: When markers predict functionality[END_REF], that are one of the first cell populations to be affected during the viral rebound [START_REF] Lévy | Effect of intermittent interleukin-2 therapy on CD4+ T-cell counts following antiretroviral cessation in patients with HIV[END_REF]. Since flow-cytometry measurements were repeated at each time point for each patient, we used the sequential strategy at week 26, in the hope to use the information from week 24 to better identify the CD4+ effector T-cell population at the next time point. Figure 4 illustrates the overall efficiency gain at week 26 from using the sequential strategy. The average limited F-measure (considering available manual gating as gold-standard) on those 12 samples is 0.58 for NPflow with a non-informative prior, and increases to 0.63 with the sequential strategy. By comparison, flowMeans (the second best method on the benchmark GvHD dataset) gives an average limited F-measure of 0.49 (see Online Supplement D for details). We also compared our approach to the HDPGMM proposed by [START_REF] Cron | Hierarchical modeling for rare event detection and cell subset alignment across flow cytometry samples[END_REF] that is specifically focusing on small cell populations (even if it had the lowest F-measure on the benchmark dataset -see Table 1). In spite of this example representing its ideal use-case, it performed slightly worse than our approach giving an average limited F-measure of 0.62 (see Online Supplement D for more details [START_REF] Hejblum | Supplement to "Sequential Dirichlet process mixtures of multivariate skew tdistributions for model-based clustering of flow cytometry data[END_REF]). Figure 5 gives an example of a patient for who the sequential strategy was especially improving the identification of the CD4+ effector T-cells. In this case, the percentage of CD4+ effector T-cells was estimated at 31.7 by the manual gating, at 7.6 by NPflow, and at 38.1 by the sequential strategy. Figure 6 shows the slight increase (of about 2%) of CD4+ effector T-cell proportions after treatment interruption (see Online Supplement D for more details [START_REF] Hejblum | Supplement to "Sequential Dirichlet process mixtures of multivariate skew tdistributions for model-based clustering of flow cytometry data[END_REF]).

In addition to providing a point estimate of the partition, our method also quantifies the uncertainty around the posterior clustering through posterior co-clustering probabilities. Figure 7 displays such a co-clustering posterior probability matrix where we can clearly identify four core clusters, with some uncertainty between two overlapping clusters that both belong to the CD4 effector T-cell population (with moderate and high expression of the CD45RA marker, respectively). proaches, our sequential strategy using Dirichlet process mixtures of multivariate skew t-distributions allowed a better recovery of the effector T-cell population after a meaningful perturbation of this targeted population following HAART interruption, highlighting their expected increase.

Our proposed method extends the classical multivariate Dirichlet process Gaussian mixture model to multivariate skew t-distribution mixtures, based on Frühwirth-Schnatter and Pyne (2010) parametrization of the restricted multivariate skew t-distribution [START_REF] Lee | On mixtures of skew normal and skew tdistributions[END_REF]. Automatic gating of cell populations from flow cytometry data is an open research problem and the proposed approach features two important characteristics for this task: i) it avoids the difficult issue of model selection by estimating directly the number of components in the mixture ; ii) it uses skew and heavy tailed distributions in the form of skew t-distributions, of which the skew normal and the normal are particular cases. Further domain-knowledge can be incorporated in the proposed model by specifying more informative priors on the Dirichlet process parameters for instance. Thanks to the use of the rMST formulation of the skew t, we avoid the identifiability issues mentioned by [START_REF] Lee | Finite mixtures of canonical fundamental skew t-distributions: The unification of the restricted and unrestricted skew t-mixture models[END_REF]. Estimation of the pairwise posterior co-clustering probabilities allows to quantify the uncertainty about the posterior partition, and an optimal point estimate of the clustering is provided by minimizing a cost function in regards to the average posterior co-clustering matrix. We have developed and implemented an efficient collapsed Metropolis within Gibbs sampler for estimating such models. One of the advantages of our proposed sampler is the absence of label switching issue, as it uses directly the partition of the data without having to deal with labels [START_REF] Jasra | Markov Chain Monte Carlo Methods and the Label Switching Problem in Bayesian Mixture Modeling[END_REF]. The computational cost of fitting our model is linear in the number of observations as well as in the number of clusters, whilst the computational cost of the partition point estimate depends of the optimal criterion chosen. The use of state-of-the-art MCMC techniques along with inner parallelization allow us to mitigate the computational cost that comes with such approaches on large data. As an indication of runtime, around 3,000 MCMC iterations can be run on average for a real dataset of around 30,000 observations over 6 dimensions, using one Intel R Xeon R x5675 processor for one hour. Besides, instead of using a partially collapse Gibbs sampler algorithm, it could be of interest to also investigate the use of sequential Monte-Carlo algorithms, especially for the sequential modeling strategy or other possible dynamic extensions of the model proposed here [START_REF] Caron | Bayesian Inference for Linear Dynamic Models With Dirichlet Process Mixtures[END_REF][START_REF] Caron | Generalized Pólya Urn for Time-Varying Pitman-Yor Processes[END_REF].

In case of repeated measurements of flow cytometry data, we propose to use sequential parametric approximations of the posterior as refined infor-mative priors. The proposed sequential analysis strategy enables to analyze each sample sequentially, as the data are acquired. It does not require to wait for the last sample to perform the automatic gating nor to analyze all data at once, but it still uses available prior knowledge. This contrasts with hierarchical extensions of the Dirichlet Process Mixture Model such as those proposed by [START_REF] Cron | Hierarchical modeling for rare event detection and cell subset alignment across flow cytometry samples[END_REF] or [START_REF] Dundar | A non-parametric Bayesian model for joint cell clustering and cluster matching: identification of anomalous sample phenotypes with random effects[END_REF], where the complete dataset must be analyzed at once. This sequential strategy allows one to analyze the samples as they are acquired, which can be useful in clinical trials where there are often intermediate analyses for instance. Moreover in large studies the size of the data can make it challenging to analyze all samples at once, and such a sequential approach then makes practical sense [START_REF] Huang | Sampling for Bayesian Computation with Large Datasets[END_REF]. Futhermore, this use of sequentially informed priors does not face the usual complications of cluster matching arising when an algorithm is run on each sample separately [START_REF] Cron | Hierarchical modeling for rare event detection and cell subset alignment across flow cytometry samples[END_REF]. In our simulation study, this sequential posterior approximation strategy improves the fit of the model. In addition, such a strategy exhibits accelerated convergence and greater accuracy for small clusters, as long as the different samples are similar enough. Besides, the parametric prior can also be specified to inform the model with expert knowledge, e.g. to favor a range for the expected number of clusters. On real flow-cytometry, data we showed that the sequential strategy also improves the clustering performances. On the benchmark dataset, it outperformed all other methods investigated in by [START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF], and in the DALIA-1 trial, the sequential strategy also improved the automatic gating results. It is worth noting however that in other cases, for instance if the data distributions were too different between samples, the sequential posterior model would not necessarily improve the clustering results, and could even gave a diminished F -measure compared to the non sequential strategy.

Manual gating is still considered the gold-standard when evaluating an automatic gating strategy on real flow cytometry data. Yet one should keep in mind that manual gating has reproducibility issues, often resulting in a partial and subjective clustering [START_REF] Ge | flowPeaks: a fast unsupervised clustering for flow cytometry data via K-means and density peak finding[END_REF][START_REF] Welters | Harmonization of the intracellular cytokine staining assay[END_REF][START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF][START_REF] Gondois-Rey | Multi-parametric cytometry from a complex cellular sample: Improvements and limits of manual versus computational-based interactive analyses[END_REF]. Therefore using manual gating as the gold-standard might not be actually the best way to assess the performance of automatic gating algorithms on real data, because of its inherent flaws.

Mass cytometry is a technology very similar to flow cytometry. Using ions in place of colors, CyTOF is able to measure up to 40 cell markers at once, generating even more data than flow cytometry. Efficient automated gating method are therefore all the more needed in the context of CyTOF [START_REF] Melchiotti | Cluster stability in the analysis of mass cytometry data[END_REF]. The approach proposed here could be directly applied to such data. More generally, we propose here a framework for Dirichlet process mixtures of multivariate skew t-distributions modeling that is suitable for any kind of data modeled as such a mixture, especially when the number of mixture components is unknown. We provide an efficient implementation of our method within the R package NPflow that is available on CRAN at https://CRAN.R-project.org/package=NPflow.

  2. Update the mixing distribution G given α, tξ k u, tψ k u, tΣ k u and the base distribution G 0 via slice sampling. 3. For c " 1, . . . , C update the individual skew parameter s c given tξ k u, tψ k u, tΣ k u and the new c . 4. Update tξ k u, tψ k u, tΣ k u given the base distribution G 0 , the updated partition t 1:C u and the updated individual skew parameters ts 1:C u. 5. Finally jointly update the degrees of freedom and the individual scale factors ptν k u, tγ 1:C uq in an Metropolis-Hastings (M-H) within Gibbs step. First an M-H step is performed to update the tν k u where the tγ 1:C u are integrated out, immediately followed by a Gibbs step to sample the tγ 1:C u from their full conditional distribution. This ensures that the reduced conditioning performed in the M-H step does not change the stationary distribution of the Markov chain (van Dyk and Jiao, 2015) -see Online Supplement A (Hejblum et al., 2018).

Fig

  Fig 1. Partition point estimate from one of the 100 2-dimensional simulations with weakly informative prior (10,000 observations with 5 clusters representing respectively 50%, 29.9%, 15%, 5% and 0.1% of the data). A: well separated clusters ; B: overlapping clusters

Fig 2 .

 2 Fig 2. Comparison of a Gaussian kernel versus a skew t kernel. A, B: well separated clusters ; C,D: overlapping clusters

Fig 4 .

 4 Fig 4. Limited F-measures for the CD4+ effector T-cell population from the DALIA-1 trial two weeks after HAART interruption for NPflow with or without the sequential strategy, compared to manual gating.

Fig 5 .

 5 Fig 5. CD3+CD4+ cells of patient 3 from the DALIA-1 trial two weeks after HAART interruption (at week 26).

Fig 6 .

 6 Fig 6. Proportion of CD4+ effector T-cells in the DALIA-1 trial following HAART interruption (from manual gating).

Fig 7 .

 7 Fig 7. Heatmap of the posterior co-clustering probabilities for the CD3+CD4+ cells of patient 3 at week 26 from DALIA-1 with NPflow seq.

  , only the most concordant clusters Fig 3. Mean F-measure according to the number of observations available at time t `1, while 1,000 observations are available at time t, over 300 simulations

		0.94
	Mean F-measure	0.90 0.92								All data (t, t+1) Data at t+1 only (prior informed on t) Data at t+1 only (non informative prior)
		0.88
		100	200	300	400	500	600	700	800	900 1000
		Number of additionnal observations available at t+1

between the eight gatings (i.e. with an F-measure above 0.8) were used for comparison with the automated results, as was done in

[START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF]

. The data were downloaded from the FlowCAP project website [http://flowcap.flowsite.org/] as part of the FlowCAP-I challenge [http: //flowcap.flowsite.org/codeanddata/FlowCAP-I.zip]. Table

Table 1

 1 Mean F-measures across all the 12 samples from the GvHD benchmark dataset

	Method	F-measure NPflow
		0.85 (0.80, 0.90)
	NPflow-seq M	0.89 (0.85, 0.94)
	NPflowG	0.61 (0.57, 0.66)
	ADICyt :	0.81 (0.72, 0.88)
	CDP :	0.52 (0.46, 0.58)
	FLAME :	0.85 (0.77, 0.91)
	FLOCK :	0.84 (0.76, 0.90)
	flowClust/Merge : 0.69 (0.55, 0.79)
	flowMeans :	0.88 (0.82, 0.93)
	FlowVB :	0.85 (0.79, 0.91)
	L2kmeans :	0.64 (0.57, 0.72)
	MM :	0.83 (0.74, 0.91)
	MMPCA :	0.84 (0.74, 0.93)
	SamSPECTRAL : 0.87 (0.81, 0.93)
	SWIFT :	0.63 (0.56, 0.70)
	HDPGMM ;M	0.35 (0.30, 0.39)

˚95% Confidence Intervals are calculated on 10,000 bootstrap samples of the F -measures.

M methods that do not analyze the 12 samples independently. : estimates from

[START_REF] Aghaeepour | Critical assessment of automated flow cytometry data analysis techniques[END_REF]

.

; estimates are from

[START_REF] Johnsson | BayesFlow: latent modeling of flow cytometry cell populations[END_REF]

.

Acknowledgements. The authors are extremely grateful to Jean-Louis Palgen for his time and efforts to manually gate the effector T-cells in the DALIA-1 trial at weeks 24 and 26, as well as to the DALIA-1 study group. The authors also thank Cedric Lachat for his precious help in getting the dpmix software to run despite its lack of maintenance. The authors also thank Nima Aghaeepour for his help in using supplementary data provided with his publication (Aghaeepour et al., 2013). Boris P. Hejblum was a recipient of a Ph.D. fellowship from the École des Hautes Études en Santé Publique (EHESP) Doctoral Network. Part of this work has been supported by the BNPSI ANR project n ˝ANR-13-BS-03-0006-01. Part of this work has been supported by the IMI2 grant EBOVAC2. Computer time for this study was partly provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour.

SUPPLEMENTARY MATERIAL

Online Supplement to "Sequential Dirichlet Process Mixtures of Multivariate Skew t-distributions for Model-based Clustering of Flow Cytometry Data" (doi: COMPLETED BY THE TYPESETTER; .pdf). We provide additional mathematical details for the proposed Gibbs samplers and the parameter estimations, as well as additional plots showing the good performance of the sequential strategy.


Boris P Hejblum

Jing Cui

Lauren J Lahey

Andrew Cagan

Jeffrey A Sparks

Jeremy Sokolove

Tianxi Cai

MD, MPH Katherine P Liao
email: kliao@partners.org

Affiliations

Association between anti-citrullinated fibrinogen antibodies and coronary artery disease in rheumatoid arthritis

Conclusion

Anti-cit-fibrinogen antibodies as a group were associated with CAD outcomes in our RA cohort, with the strongest signal for association arising from a subset of the autoantibodies.

Significance and Innovations

• Previous studies found an association between anti-citrullinated protein antibodies (ACPA) directed against citrullinated fibrinogen (cit-fibrinogen) and surrogate markers of cardiovascular (CV) risk; this study examined the association between ACPAs and coronary artery disease (CAD) outcomes

• Our study demonstrated a specific association between a subset of anti-cit-fibrinogen antibodies and CAD in RA, present even after adjusting for anti-CCP status • We provide a roadmap for a data driven approach that can be applied to study association between a group of biomarkers with an outcome using state-of-the-art biostatistics methods; the traditional approach of testing individual biomarkers with an outcome and adjusting for multiple comparisons are challenging due to the high degree of correlation between ACPAs

INTRODUCTION

Cardiovascular disease (CVD) is a leading cause of death for patients with rheumatoid arthritis (RA) [START_REF] Avina-Zubieta | Risk of incident cardiovascular events in patients with rheumatoid arthritis: a meta-analysis of observational studies[END_REF].

Several studies have demonstrated that traditional CV risk calculators significantly underestimate CV risk in RA [START_REF] Crowson | Usefulness of risk scores to estimate the risk of cardiovascular disease in patients with rheumatoid arthritis[END_REF]. These findings have motivated studies to examine whether biomarkers can assist rheumatologists in identifying RA patients at elevated CV risk. Anti-citrullinated protein antibodies (ACPAs) are used in the diagnosis of RA, and recent studies have identified epitope targets several epitope targets implicated in the disease [START_REF] Sokolove | Autoantibody epitope spreading in the pre-clinical phase predicts progression to rheumatoid arthritis[END_REF][START_REF] Brink | Rheumatoid factor isotypes in relation to antibodies against citrullinated peptides and carbamylated proteins before the onset of rheumatoid arthritis[END_REF][START_REF] Raijmakers | Elevated levels of fibrinogen-derived endogenous citrullinated peptides in synovial fluid of rheumatoid arthritis patients[END_REF]. A specific class of ACPAs, autoantibodies targeting citrullinated fibrinogen (anti-cit-fibrinogen) have been studied as potential markers for cardiovascular disease (CVD) in RA [START_REF] Sokolove | Brief report: citrullination within the atherosclerotic plaque: a potential target for the anti-citrullinated protein antibody response in rheumatoid arthritis[END_REF].

In a previously published study, Sokolove an colleagues observed that protein lysates of atherosclerotic plaques from aortic specimens were enriched for citrullinated fibrinogen (cit-fibrinogen) [START_REF] Sokolove | Brief report: citrullination within the atherosclerotic plaque: a potential target for the anti-citrullinated protein antibody response in rheumatoid arthritis[END_REF]. Further, they found that RA patients with higher titers of anti-cit-fibrinogens also had higher atherosclerotic burden as measured by coronary artery calcium (CAC) scores. These data suggested a mechanistic pathway whereby inflammation promoted cit-fibrinogen in atherosclerotic plaques, and thus the development of anti-cit-fibrinogen may be a useful marker of inflammation, atherosclerotic burden and CV risk in RA. Interestingly, a subsequent study from an independent cohort using a different platform to measure anti-cit-fibrinogen found no association between carotid intima media thickness or CAC, and titers of anti-cit-fibrinogen [START_REF] Montes | Brief report: lack of replication of an association between anti-citrullinated fibrinogen and subclinical atherosclerosis in patients with rheumatoid arthritis[END_REF].

Since previous studies have only examined surrogate CV outcomes, the objective of this study was to examine the association between anti-cit-fibrinogen and the outcome of CAD. Additionally, fibrinogen is a relatively large molecule [START_REF] Kollman | Crystal structure of human fibrinogen[END_REF] with autoantibodies targeting different portions of the citrullinated fibrinogen. Thus, a secondary objective of this study was to determine whether there was specificity to the type(s) of anti-cit-fibrinogen associated with CAD.

METHODS

Study population

We studied a validated RA cohort identified in the electronic medical records (EMR) of two large tertiary care hospitals in Boston, MA, Brigham and Women's Hospital (BWH) and Massachusetts General Hospital (MGH). Briefly, the RA cases were identified using a validated RA phenotype algorithm utilizing a combination of structured data, e.g. International Classification of Diseases, 9 th Revision (ICD9), and concepts extracted from text notes, e.g. smoking status, using natural language processing (NLP). The positive predictive value (PPV) for RA using the algorithm was 94% [START_REF] Liao | Electronic medical records for discovery research in rheumatoid arthritis[END_REF]. Discarded blood samples were collected using a biospecimen repository linked to anonymized EMR data and all samples were measured for antibodies to cyclic citrullinated peptide (anti-CCP) using a commercial kit as described in a previous publication [START_REF] Kurreeman | Genetic basis of autoantibody positive and negative rheumatoid arthritis risk in a multi-ethnic cohort derived from electronic health records[END_REF]. We studied all RA subjects with clinical data and available blood samples.

CAD Outcome definition

Subjects were classified as having CAD if they had diagnosis of CAD in the medical record with supporting evidence of disease through documentation of CABG, percutaneous coronary intervention (PCI) with stent or balloon angioplasty, positive stress test, or EKG changes consistent ischemia [START_REF] Liao | Methods to Develop an Electronic Medical Record Phenotype Algorithm to Compare the Risk of Coronary Artery Disease across 3 Chronic Disease Cohorts[END_REF].

Subjects with a diagnosis of CAD without supporting documentation or subjects with no mention of CAD were considered not to have the disease.

Covariates

Information on covariates and potential confounders were extracted from the structured EMR data including age, gender, self-reported race classified as white or non-white, hypertension (ICD9 401.x, 402.x. 403.x. 404.x, 405.x) , diabetes (249.xx or 250.xx), hyperlipidemia (272.xx), and RA treatments (ever/never use): disease modifying anti-rheumatic drugs (DMARDs), categorized into non-biologic and biologic DMARDS, and prednisone. Smoking status was extracted using NLP [START_REF] Liao | Electronic medical records for discovery research in rheumatoid arthritis[END_REF].

Laboratory testing

All plasma samples were measured for anti-cit-fibrinogen using a bead-based antigen microarray as described previously [START_REF] Sokolove | Brief report: citrullination within the atherosclerotic plaque: a potential target for the anti-citrullinated protein antibody response in rheumatoid arthritis[END_REF]. In total, we measured autoantibodies directed against 10 epitopes of citfibrinogen (Table 2). Anti-CCP positivity was determined using a commercially available assay (CCP3 IgG ELISA, INOVA Diagnostics, San Diego, CA).

Statistical analysis

To test the association between autoantibodies as a group with the primary outcome CAD, we used the score test [START_REF] Vasan | Biomarkers of cardiovascular disease: molecular basis and practical considerations[END_REF][START_REF] Liao | Phenome-Wide Association Study of Autoantibodies to Citrullinated and Noncitrullinated Epitopes in Rheumatoid Arthritis[END_REF]. The use of the score test allowed us to integrate information from multiple correlated ACPAs at once (which individually may not be statistically significant) to determine the overall association of the group with the outcome. The p-value from the score test provides a summary statistic for whether the biomarkers as a group are associated with the outcome. Prior studies on ACPAs tested the association between individual ACPAs with an outcome. Due to the high degree of correlation between the ACPAs, this approach cannot account for collinearity between ACPAs, and can further result in a loss of power to detect associations. As well, each individual ACPA may have a modest effect size, however, as a group, the ACPAs may have a significant association with an outcome that may be missed with individual testing. The score test provides information about whether there is a significant association between a group and an outcome. Notably, since it is testing a group, a single effect size cannot be provided for the group.

Since the anti-cit-fibrinogens are known to be highly correlated, we performed a Principal Component Analysis (PCA) to determine whether anti-cit-fibrinogen data cluster into groups based on the data.

Briefly, PCA is a factor analysis method that highlights linear relationship between correlated variables while extracting a low number of key Principal Components (PC). The first PC explained more than 80% of the anti-cit-fibrinogen antibody variability. PC2 explained an additional 7.7% of variability and PC3 an additional 3.5% of variability. Together the 3 PCs explained over 90% of the variability in anti-citfibrinogen levels. From the PCA, anti-cit-fibrinogen autoantibodies clustered into three groups (Figure 1). In addition, the same clustering was found when performing a hierarchical clustering of the data (data not shown). We then tested the association between each group of anti-cit fibrinogen with CAD outcomes using a score test in a logistic regression model, adjusting for traditional CV risk factors.

We performed 3 sensitivity analyses. First, we performed an association study using a published group of ICD9 codes for ischemic heart disease (IHD, ICD9 411.x or 414.x) [START_REF] Denny | Systematic comparison of phenome-wide association study of electronic medical record data and genome-wide association study data[END_REF], where subjects with ³1 ICD9 code for any of those codes was defined as having IHD. Second, we tested for association between the anti-cit-fibrinogen group with CAD outcome adjusting for the same covariates in the primary analysis with the addition of high sensitivity CRP (hsCRP) levels. Third, we further added anti-CCP to the adjusted model.

Statistical analyses were performed in R 3.2.0 (R Foundation for Statistical Computing, Vienna, Austria).

All aspects of this study were approved by the Partners Institutional Review Board (IRB).

RESULTS

We studied a total of 1006 RA patients with mean age 61.0 (SD) years, 79.0% females, and 72.2% anti-CCP positive (Table 1).

From the PCA analysis PC1 accounts for the majority of variability, 81.4%, among all anti-cit fibrinogen antibodies(Figure 1). PC2 explains 7.7% and PC3 3.5% of the variability across all anti-cit fibrinogens respectively. From this PCA analysis, we observed that the anti-cit-fibrinogen antibodies clustered into three main groups (Table 2). The anti-cit-fibrinogens (10 autoantibodies total) as a group were associated with a higher risk of CAD, adjusted for age, gender, race, diabetes, hypertension, hyperlipidemia, smoking, non-biologic DMARD, biologic DMARD and prednisone use, p=1.1x10 -4 (Table 2). When the associations were examined by subgroups, the strongest association was observed in Group 2 (p=0.015).

The sensitivity analysis using ICD9 codes showed a similar association between anti-cit-fibrinogen as a group with IHD, as well as between Group 2 and IHD (Table 2). We observed no association between anti-CCP positivity (commercial assay) with CAD, OR 1.3 (95% CI 0.7, 2.3), using the same covariates for adjustment. Further we studied the subjects with hsCRP data (n=897). The association between anti-cit fibrinogens and CAD remained significant adjusting for CV risk factors and hsCRP (p=2.2x10 -4 , for Group 2, p=7.1x10 -3 ). As well, the association remained significant after adjusting the model including hsCRP and anti-CCP status (p=5.7x10 -8 ; for group 2 alone, p=2.8x10 -4 ).

Discussion

In this study, we observed a significant association between anti-cit-fibrinogens with CAD outcomes in a cohort of RA patients. Using a data-driven approach, we observed that the association was driven by a subset of autoantibodies in Group 2 (Table 2). These findings demonstrated an association between anti-cit-fibrinogen with the clinical outcome of CAD, adding to previous studies examining surrogate measures of CV risk. Moreover, this study may explain why differences have been observed with regards to association between anti-cit-fibrinogen and surrogate measures of CV risk [START_REF] Sokolove | Brief report: citrullination within the atherosclerotic plaque: a potential target for the anti-citrullinated protein antibody response in rheumatoid arthritis[END_REF][START_REF] Montes | Brief report: lack of replication of an association between anti-citrullinated fibrinogen and subclinical atherosclerosis in patients with rheumatoid arthritis[END_REF]. While we detected a significant association between all anti-cit-fibrinogens with CAD, the signal appeared to be driven by a group of targets within the large fibrinogen molecule. Thus, if a particular ACPA assay does not measure antibodies targeting specific regions of fibrinogen, an association with CAD, may not be detected.

Additionally, an association may not be detected using individual ACPA tests.

In a study by Montes et al [START_REF] Montes | Brief report: lack of replication of an association between anti-citrullinated fibrinogen and subclinical atherosclerosis in patients with rheumatoid arthritis[END_REF] no relationship was observed between anti-cit-fibrinogen and cIMT and CAC scores. Their examined antibodies directed against citrullinated fibrinogen beta and not alpha. In the present study, we found the strongest signal to be comprised of a group of autoantibodies directed at citrullinated chains for both fibrinogen alpha and beta. Using the same biomarker measurement platform, the results of this study corroborates the previous associations observed between anti-cit fibrinogen and CAC scores [START_REF] Sokolove | Brief report: citrullination within the atherosclerotic plaque: a potential target for the anti-citrullinated protein antibody response in rheumatoid arthritis[END_REF]. This association remained significant after adjusting for potential confounders including age, gender, race, and traditional CV risk factors. The association remained significant after further adjustment for CRP and anti-CCP status. Thus, the concentrations of anti-cit fibrinogen provided additional information in identifying which RA patients had CAD in addition to anti-CCP. These data also support the notion that specificity for epitope target may be important for identifying CAD in RA; the commercial anti-CCP assay is a mixture of cyclic antibodies against a proprietary set of epitope targets, whereas the anti-cit fibrinogen assays were developed using tissue from RA patients with characterized epitope targets. As found in a previous study, cit-fibrinogen is present in atherosclerotic plaques and are recognized by ACPAs from RA patients [START_REF] Sokolove | Brief report: citrullination within the atherosclerotic plaque: a potential target for the anti-citrullinated protein antibody response in rheumatoid arthritis[END_REF].

We designed this study as a "Phase 3" biomarker validation study. Phase 3 biomarker studies focus on the association between previously identified candidate biomarkers, anti-cit fibrinogens with CAD outcome, rather than novel biomarker identification [START_REF] Vasan | Biomarkers of cardiovascular disease: molecular basis and practical considerations[END_REF]. Therefore, we did not perform screening for other ACPAs with the CAD outcome. However, our recently published results from a Phenome-Wide Association Study (PheWAS) corroborate findings from previous studies as well as our own. This hypothesis generating screen performed for all ACPAs and all potential phenotypes (defined using billing code data), showed a signal for association between anti-cit fibrinogens with billing codes for ischemic heart disease [START_REF] Liao | Phenome-Wide Association Study of Autoantibodies to Citrullinated and Noncitrullinated Epitopes in Rheumatoid Arthritis[END_REF]. We did not observe association of other ACPAs with CV-related phenotypes.

The potential mechanism or biology behind the association between cit-fibrinogen, anti-cit-fibrinogen and CV risk is beyond the scope of this study. However, fibrinogen is a large multi-domain glycoprotein and is one of the best characterized citrullinated protein targets of ACPAs in RA [START_REF] Kollman | Crystal structure of human fibrinogen[END_REF][START_REF] Fernandes-Cerqueira | Targeting of anti-citrullinated protein/peptide antibodies in rheumatoid arthritis using peptides mimicking endogenously citrullinated fibrinogen antigens[END_REF]. Fibrinogen is a known constituent of the atherosclerotic plaque for which at least a portion has been found to be citrullinated [START_REF] Smith | Fibrinogen, fibrin and fibrin degradation products in relation to atherosclerosis[END_REF]. These data suggest a potential pathway connecting citrullination of fibrinogen with CAD. Whether pathology occurs as a result of citrullination of fibrinogen or antibody binding of citrullinated fibrinogen is not known.

A limitation to our study is that the samples were collected in patients who have prevalent CAD. Thus, this study cannot answer the question of whether anti-cit-fibrinogen is predictive of future risk.

Notably, the association between the anti-cit fibrinogen as a group remained significant after adjusting for anti-CCP status. Misclassification of RA and CAD is another potential limitation. However, we believe that misclassification, particularly with CAD is more likely to attenuate the association.

Additionally, while we adjusted for CRP, we did not have detailed information on disease activity. In summary, previous studies observed an association between anti-cit-fibrinogen antibodies with preclinical CAD, and in this study with actual CAD outcomes supporting its utility as a potential marker for assessing CV risk in RA. While our study adjusted for traditional risk CV risk factors, future Phase 4 biomarker development studies investigating the predictive value of these biomarkers in an independent cohort are needed. Further studies can assess whether these markers as a group improve

Figure 1 .

 1 Figure 1. Principal component analysis of anti-cit-fibrinogen antibodies demonstrating 3 main groupings of anti-cit-fibrinogen antibodies. (FA=fibrinogen alpha, FB= fibrinogen beta)

Table 1 .

 1 Clinical characteristics of cohort (n=1006).

	Clinical characteristics	CAD yes, n=93 (9.2%)	CAD no, n=913 (90.8%)
	Mean age, years (SD)	71.0 (10.0)	60.0 (12.8)
	Female, n (%)	56 (60.2)	739 (80.9)
	Anti-CCP positive, n (%)	69 (74.2)	657 (72.0)
	HTN	74 (79.6)	321 (35.2)
	Hyperlipidemia	64 (68.8)	246 (26.9)
	Diabetes mellitus	29 (31.2)	84 (9.2)
	Mean hsCRP*, mg/L (SD)	13.1 (17.3)	11.4 (18.0)
	DMARD, ever use, n (%)	87 (93.5)	865 (94.7)
	Biologic DMARD, ever use, n (%)	43 (46.2)	437 (47.9)
	Non-biologic DMARD, ever use, n (%)	62 (66.7)	564 (61.8)
	Methotrexate	50 (53.8)	589 (64.5)
	TNFi	41 (44.1)	417 (45.7)
	Smoking, ever, n (%)	73 (78.5)	496 (54.3)
	Prednisone, ever use, n (%)	82 (88.2)	682 (74.7)
	* hsCRP value available for n=897 subjects

Table 2 .

 2 The association between anti-cit-fibrinogens and CAD outcomes*. Adjusted by age, gender, self-reported race, hypertension, diabetes mellitus, hyperlipidemia, nbDMARD ever use, bDMARD ever use, smoking status and prednisone ever use.

	Group	Anti-cit-fibrinogen	Association with CAD	Association with IHD
			p-value	p-value
	Group 1	1: FA (616635) cit3 sm-cyclic	0.33	0.65
	(red)	2: FA (556-575) cit sm cyclic
	Group 2	3: FA (211-230) cit sm cyclic	0.015	0.0040
	(orange)	4: FB (54-72) cit1
		5: FA (41-60) cit3 cyclic
		6: FB (36-52) cit
		7: FA (27-43) cit
		8: FB (246-267) cit
		9: FA (582-599) cit
	Group 3	10: Fibrinogen cit	0.34	0.29
	(blue)
	All anti-cit-	All 10 autoantibodies above	0.00011	0.00029
	fibrinogen

*

Note: FA= Fibrinogen alpha, FB=Fibrinogen beta, sm=small

Rheumatology Research Foundation Disease Targeted Research Pilot Grant, NIH R01 HL127118, NIH U54 HG007963

CV risk stratification for RA patients in the context of validated general population based CV risk estimators.


Z V Khaidukov

M A Zubkov

Chiral torsional effect

Keywords: numbers: 71, 10, -w, 67, 30, H-

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The non -dissipative transport effects have been widely discussed recent years [START_REF] Landsteiner | Anomalous Transport from Kubo Formulae[END_REF][START_REF] Chernodub | A condensed matter realization of the axial magnetic effect[END_REF][START_REF] Gorbar | Chiral separation and chiral magnetic effects in a slab: The role of boundaries[END_REF][START_REF] Miransky | Quantum field theory in a magnetic field: From quantum chromodynamics to graphene and Dirac semimetals[END_REF][START_REF] Valgushev | Chiral Magnetic Effect in finite-size samples of parity-breaking Weyl semimetals[END_REF][START_REF] Buividovich | Chiral magnetic conductivity in an interacting lattice model of parity-breaking Weyl semimetal[END_REF][START_REF] Buividovich | Spontaneous chiral symmetry breaking and the Chiral Magnetic Effect for interacting Dirac fermions with chiral imbalance[END_REF][START_REF] Buividovich | Anomalous transport with overlap fermions[END_REF]. These effects are to be observed in the non -central heavy ion collisions [9]. They have also been considered for the Dirac and Weyl semimetals [START_REF] Parameswaran | Probing the chiral anomaly with nonlocal transport in Weyl semimetals[END_REF][START_REF] Vazifeh | Electromagnetic response of weyl semimetals[END_REF][START_REF] Chen | Axion response in Weyl semimetals[END_REF][START_REF] Chen | Weyl fermions and the anomalous Hall effect in metallic ferromagnets[END_REF][START_REF] Ramamurthy | Patterns of electro-magnetic response in topological semi-metals[END_REF][START_REF] Zyuzin | Topological response in Weyl semimetals and the chiral anomaly[END_REF][START_REF] Goswami | Axionic field theory of (3+1)-dimensional Weyl semi-metals[END_REF][START_REF] Liu | Chiral gauge field and axial anomaly in a Weyl semimetal[END_REF] and in 3 He-A [START_REF] Volovik | The Universe in a Helium Droplet[END_REF]. It was expected that among the other effects their family includes the chiral separation effect (CSE) [START_REF] Metlitski | Anomalous Axion Interactions and Topological Currents in Dense Matter[END_REF], the chiral magnetic effect (CME) [20][START_REF] Kharzeev | The Chiral Magnetic Effect and Anomaly-Induced Transport[END_REF][START_REF] Kharzeev | Chiral Magnetic Effect in High-Energy Nuclear Collisions -A Status Report[END_REF][START_REF] Kharzeev | Chern-Simons current and local parity violation in hot QCD matter[END_REF], the chiral vortical effect (CVE) [24], the anomalous quantum Hall effect (AQHE) [START_REF] Zyuzin | Topological response in Weyl semimetals and the chiral anomaly[END_REF]25,26]. All those phenomena have the same origin -the chiral anomaly. The naive "derivation" of the CME from the standard expression for the chiral anomaly has been presented, for example, in [START_REF] Zyuzin | Topological response in Weyl semimetals and the chiral anomaly[END_REF]20,[START_REF] Son | Berry curvature, triangle anomalies, and chiral magnetic effect in Fermi liquids[END_REF].

However, the more accurate consideration demonstrates that the original equilibrium version of the CME does not exist. In [START_REF] Valgushev | Chiral Magnetic Effect in finite-size samples of parity-breaking Weyl semimetals[END_REF][START_REF] Buividovich | Chiral magnetic conductivity in an interacting lattice model of parity-breaking Weyl semimetal[END_REF][START_REF] Buividovich | Spontaneous chiral symmetry breaking and the Chiral Magnetic Effect for interacting Dirac fermions with chiral imbalance[END_REF][START_REF] Buividovich | Anomalous transport with overlap fermions[END_REF] this subject has been elaborated via the numerical lattice methods. In the condensed matter theory the absence of CME was reported in [31]. The same conclusion has been obtained on the basis of the no -go Bloch theorem [32]. The analytical proof of the absence of the equilibrium CME was presented in [START_REF] Zubkov | Wigner transformation, momentum space topology, and anomalous transport[END_REF][START_REF] Zubkov | Absence of equilibrium chiral magnetic effect[END_REF] using the technique of Wigner transformation [START_REF] Zachos | Quantum Mechanics in Phase Space[END_REF][START_REF] Robert G Littlejohn | The semiclassical evolution of wave packets[END_REF][38][START_REF]Wigner transformation, momentum space topology, and anomalous transport M.A Zubkov[END_REF] applied to the lattice models 1 . The same technique allows to reproduce the known results on the AQHE [START_REF] Zubkov | Wigner transformation, momentum space topology, and anomalous transport[END_REF]. It also confirms the existence of the equilibrium CSE [START_REF] Khaidukov | Chiral Separation Effect in lattice regularization[END_REF]. Notice, that in the framework of the naive nonregularized quantum field theory the CSE was discussed, for example, in [START_REF] Gorbar | Chiral separation and chiral magnetic effects in a slab: The role of boundaries[END_REF] in the technique similar to the one that was used for the consideration of the CME [START_REF] Zyuzin | Topological response in Weyl semimetals and the chiral anomaly[END_REF]20].

The investigation of momentum space topology was developed previously mainly within the condensed matter physics. It allows to relate the gapless boundary fermions to the bulk topology for the topological insulators and to describe the stability of the Fermi surfaces of various kinds. Besides, this technique allows to calculate the number of the fermion zero modes on vortices (for the review see [START_REF] Volovik | The Universe in a Helium Droplet[END_REF][START_REF] Volovik | Topology of quantum vacuum[END_REF]). Momentum space topology of QCD has been discussed recently in [START_REF] Zubkov | Momentum space topology of QCD[END_REF]. The whole Standard Model of fundamental interaction was discussed in this framework in [START_REF] Volovik | Standard Model as the topological material[END_REF].

Recently the new non -dissipative transport effects were proposed: the rotational Hall effect [START_REF] Zubkov | Hall effect in the presence of rotation[END_REF] and the scale magnetic effect [START_REF] Chernodub | Anomalous Transport Due to the Conformal Anomaly[END_REF]. In the present paper we will propose one more non -dissipative transport effect -the chiral torsional effect. Namely, we will discuss the emergence of axial current in the presence of torsion. It will be shown that this effect is intimately related to the chiral vortical effect [24]. In conventional general relativity [START_REF] Diakonov | Lowenergy General Relativity with torsion: a systematic derivative expansion[END_REF] torsion vanishes identically, it appears only in its various extensions. The background (non -dynamical) gravity with torsion emerges, in particular, in certain condensed matter systems. For example, elastic deformations in graphene and in Weyl semimetals induce the effective torsion experienced by the quasiparticles [START_REF] Volovik | [END_REF][START_REF] De Juan | Dislocations and torsion in graphene and related systems[END_REF]. In 3 He-A torsion appears dynamically when motion of the superfluid component is non -homogeneous. Below we will consider mostly the massless free Dirac fermions in the presence of the background gravity with torsion. We assume that this model is the relevant ingredient of the low energy effective field theory for the appropriate condensed matter system (either the Weyl/Dirac semimetal or 3 He-A) or of the theory of elementary particles in the presence of real background gravity with torsion. The other ingredients of the mentioned theories are interactions between the fermions due to the exchange by various interaction carriers including the gauge bosons. On the present level of consideration we omit the interactions completely.

Besides, we disregard the natural anisotropy existing without emergent gravity in 3 He-A and Weyl semimetals. Therefore, while applying expression for the chiral torsional effect obtained below to these materials one should take care of this anisotropy. It actually may be restored in the final expression in the straightforward way.

RESPONSE OF THE AXIAL CURRENT TO TORSION

We are considering the model of massless Dirac fermions, which serves as the low energy approximation for the electronic excitations in Weyl semimetals and also for the fermionic excitations in 3 He-A. For simplicity we will restrict ourselves by one flavor of Dirac fermions, and do not take into account the anisotropy of the Fermi velocity. In order to take it into account we should simply rescale the coordinates in an appropriate way.

Following [START_REF]Wigner transformation, momentum space topology, and anomalous transport M.A Zubkov[END_REF][START_REF] Khaidukov | Chiral Separation Effect in lattice regularization[END_REF] we will calculate here the response of the axial current j 5µ (x) = ψ(x)γ µ γ 5 ψ(x) to torsion in terms of the Wigner transformation of the Green function. We consider the case of vanishing spin connection and the nontrivial torsion encoded in the vielbein e a µ (x) = δ a µ + δe a µ (x). The inverse vielbein is denoted here by

E µ a = δ µ a + δE µ a , it obeys equation e a µ (x)E µ b (x) = δ a b , which gives δe a b (x) ≈ -δE a b . Momentum operator is to be substituted by p a → (δ b a +δE b a (x))p b ≈ (p a -δe b a (x)p b )
. Product of the coordinate -depending vielbein and momentum operator should be symmetrized in order to lead to the Hermitian Hamiltonian of the quasipaticles 2 . The term δe b i p b = A i to a certain extent is similar to the gauge field (see also [51,[START_REF] Sumiyoshi | Torsional Chiral Magnetic Effect in a Weyl Semimetal with a Topological Defect[END_REF]). Suppose, that G(p) is the two point fermion Green function in momentum space for the system of Dirac spinors. In order to obtain expressions for the axial current we deform it as follows:

G(p) → G(p -B(x)γ 5 ) (1)
Next, the axial current is to be expressed as the functional derivative of the partition function with respect to B. The partition function is built using the Dirac operator defined in momentum space as G -1 (p -B(x)γ 5 ). In this consideration x is to be treated as parameter. Following the same steps as those of [START_REF] Khaidukov | Chiral Separation Effect in lattice regularization[END_REF] in the leading order in the derivative of the vielbein we obtain the following expression for the axial current:

j 5k = i 2 T a ij ∫ dωd 3 p (2π) 4 p a Tr G∂ pi G -1 ∂ pj Gγ 5 ∂ p k G -1 (2) T a ij = ∂ i e a j -∂ j e a i
Let us consider the situation when the only nonzero components of torsion are given by:

T 0 ij = ϵ 0 .f ij S f (3)
2 In this form gravity emerges in graphene and Weyl semimetals.

where S is a certain three -vector. In case of the noninteracting theory of massless Dirac fermions in continuous space in the presence of finite temperature T we obtain the following expression:

j 5k = 4S k T ∑ ωn ∫ d 3 p (2π) 3 ω 2 n (ω 2 n + p 2 ) 2 (4)
Here the sum is over the discrete Matsubara frequencies ω n = 2π(n + 1/2)/T , n ∈ Z. In the expression of Eq. ( 2) we do not see the signature of the topological protection from various corrections. Therefore, we do not exclude that the small deformations of the Green's function may lead to the variation of the axial current. Thus, there will, possibly, appear corrections to this expression due to interactions.

SUM OVER THE MATSUBARA FREQUENCIES

The obtained above expression ( 4) is divergent, and requires regularization. Below we present the results of our calculations made using the two complementary methods: via the zeta-regularization and via the sum over the Matsubara frequencies with the contribution of vacuum subtracted. In both cases the results coincide. Let us consider first the direct summation over the Matsubara frequencies:

j 5z = 4S 3 T ∑ ωn ∫ d 3 p (2π) 3 ω 2 n (ω 2 n + p 2 ) 2 = j 5z vac + j 5z medium (5)
where ω n = 2πT (n + 1/2). Let us perform summation over n:

4T ∑ ωn ω 2 n (ω 2 n + p 2 ) 2 = -2 ∫ dz 2πi z 2 (z 2 -p 2 ) 2 th ( z 2T ) = 2z 2 (z + p) 2 d dz th ( z 2T ) z=p + 2z 2 (z -p) 2 d dz th ( z 2T ) z=-p + 1 p th ( p 2T ) = 1 2T ch 2 ( p 2T ) + 1 p th ( p 2T ) (6)
In the first row the integral is over the contour surrounding all poles of the hyperbolic tangent. This contour is then deformed to surround points ±p. In this expression we separate the contribution of thermal quasiparticles subtracting the contribution of vacuum (the expression of Eq. ( 5) for T = 0). We denote x = p 2T and obtain:

j 5i medium = 2T 2 S i ∫ ∞ 0 x 2 dx π 2 ( 1 ch 2 x + 1 x thx - 1 x ) ( 7
)
Vacuum part is formally divergent, and we do not consider it here. After integration over x we obtain the following expression for the separate contribution of medium (i.e. of the thermal quasiparticles):

j 5k medium = T 2 12 S k = - T 2 24 ϵ 0kij T 0 ij (8)
It is well -known that ζ -regularization subtracts systematically all divergencies. It allows to include all logarithmic divergencies to the renormalization of coupling constants automatically. At the same time the powerlike divergencies are subtracted from all expressions completely. Formally integration over momenta in Eq. ( 5) gives

⃗ j 5 = 1 2π ⃗ ST ∑ ∞ n=-∞ |ω n |.
In zeta -regularization instead of this expression we obtain

⃗ j 5 = 1 π ⃗ ST ∞ ∑ n=0 1 ω s n s=-1
Here the sum over n determines the function of s > 1.

The final answer is given by its analytical continuation to s = -1. This results in the following resummation

⃗ j 5 = 2 ⃗ ST 2 ( ∞ ∑ n=0 1 (n + 1/2) s ) s=-1 = 2 ⃗ ST 2 ζ H (-1, 1/2) = T 2 12 ⃗ S ( 9
)
In the second row ζ H (s, ν) is the Hurwitz zeta -function.

One can see that the expressions for the current produced by the thermal quasiparticles of ( 5) and ( 9) coincide.

CVE EFFECT AS THE PARTICULAR CASE OF THE CHIRAL TORSIONAL EFFECT

Let us consider the particular non -relativistic system (having the emergent relativistic symmetry at low energies) in the presence of macroscopic motion with velocity v. When v is constant or slowly varying, the one -particle Hamiltonian is given by Ĥ = Ĥ0 + pv. Here Ĥ0 is the Hamiltonian in the reference frame accompanying the motion. For the case of emergent massless Dirac fermions Ĥ0 = γ 0 ∑ k=1,2,3 γ k p k . Then the combination p 0 -Ĥ may be written as γ 0 E µ a γ a p µ with the inverse vierbein

E 0 0 = 1, E k a = δ k a , E k 0 = -v k , E 0 a = 0 for a, k = 1, 2, 3
. Therefore, we are able to treat velocity of the macroscopic motion as the components of the vielbein:

e 0 k = v k , k = 1, 2, 3. Correspondingly, torsion appears with the nonzero components T 0 ij = ∂ i v j -∂ j v i = 2ϵ 0ijk ω k
, where ⃗ ω is vector of angular velocity. Eq. ( 8) gives

j 5k medium = - T 2 24 ϵ 0kij T 0 ij = T 2 6 ω k (10)
which is nothing but the conventional expression for the chiral vortical effect (see also [START_REF] Landsteiner | Gravitational Anomaly and Transport[END_REF]). Recall that the latter effect for vanishing chemical potential is the appearance of the axial current in the system of rotating thermal quasiparticles. Typically in the consideration of this effect the vacuum as a whole is assumed to be at rest, only the excitations over vacuum are rotated. Therefore, indeed, the subtraction performed above (of the vacuum contribution to the chiral torsional effect) effectively results in the description of the chiral vortical effect for vanishing µ.

CONCLUSIONS AND DISCUSSIONS

To conclude, we have calculated the axial current produced by thermal quasiparticles in the presence of torsion for the system of the noninteracting massless Dirac fermions. Eq. ( 8) is obtained in two ways: via zeta regularization and via the summation over Matsubara frequencies with the vacuum contribution subtracted. Although we did not consider the effects of interactions directly, there are certain indications, that the coefficient in Eq. ( 8) is subject to the renormalization due to interactions. Notice, that in Eq. ( 8) there is no general coordinate invariance. First of all, this invariance is absent because we consider the system with vanishing spin connection, which corresponds to the case of emergent gravity in real Weyl materials [START_REF] Cortijo | Emergent gravity in the cubic tight-binding model of Weyl semimetal in the presence of elastic deformations[END_REF]. This corresponds to the partial gauge fixing, which eliminates the general coordinate invariance. Moreover, we considered the system with finite temperature. Therefore, even the Lorentz invariance is broken and Eq. ( 8) actually possesses the O(3) symmetry of spacial rotations.

It is also worth mentioning that we considered only the components of torsion T 0 ij with i, j = 1, 2, 3. The possible appearance of the axial current in the presence of the other components of torsion remains open and requires an additional consideration, which is out of the scope of the present paper. The resulting current of Eq. ( 8) represents the new anomalous transport effect. The corresponding current is non -dissipative, at least in the approximation of the slowly varying vielbein/torsion. The explanation for this is the mentioned above analogy between the composition δe 0 i ω and the electromagnetic vector -potential. According to this analogy H k = T 0 ij ϵ 0ijk ω plays the role of the effective magnetic field. The latter does not perform the work and thus cannot lead to dissipation.

The divergency of the vacuum axial current in the presence of background torsion prompts that the quantum field theory in the presence of background gravity with torsion is ill -defined. Presumably, there are the aspects of the theory, in which it depends on regularization. In the other words, suppose, that the microscopically different real physical systems have naively the same continuum low energy effective field theory. In the presence of emergent torsion these effective theories become dif-ferent. In particular, we expect the appearance of such differences between the low energy description of electrons in Weyl semimetals (in the presence of elastic deformations) and the description of the quasiparticles in the 3 He-A superfluid (in the presence of the in -homogeneous motion of the superfluid component). This difference may, possibly, be related to the appearance of the Nieh -Yan -like term [54] in the chiral anomaly for the solids and its absence for 3 He-A. Notice, that experiments with vortices in 3 He-A prompt that the chiral anomaly in this superfluid does not contain the Nieh -Yan term (see [START_REF] Volovik | The Universe in a Helium Droplet[END_REF] and references therein). At the same time the recent theoretical consideration of Weyl semimetals indicates the appearance of this term in the chiral anomaly for these materials [START_REF] Parrikar | Torsion, parity-odd response, and anomalies in topological states[END_REF].

It would be interesting to apply the lattice regularization in the spirit of [START_REF] Khaidukov | Chiral Separation Effect in lattice regularization[END_REF] for the calculation of axial current of massless Dirac fermions in the presence of nontrivial vielbein. Such a calculation may shed led on various questions related to the vacuum contribution to the chiral torsional effect and to the possible appearance of the Nieh -Yan term in the chiral anomaly. Z.V.K. is grateful to the physics department of the Ariel University for kind hospitality. The work of Z.V.K. was supported by Russian Science Foundation Grant No 16-12-10059. Both authors kindly acknowledge useful discussions with M.Suleymanov.