Ana Arribas-Gil

Catherine Matias
email: catherine.matias@math.cnrs.fr

A time warping approach to multiple sequence alignment

Keywords: Alignment, Dynamic time warping, Multiple sequence alignment, Warping

We propose an approach for multiple sequence alignment (MSA) derived from the dynamic time warping viewpoint and recent techniques of curve synchronization developed in the context of functional data analysis. Starting from pairwise alignments of all the sequences (viewed as paths in a certain space), we construct a median path that represents the MSA we are looking for. We establish a proof of concept that our method could be an interesting ingredient to include into refined MSA techniques. We present a simple synthetic experiment as well as the study of a benchmark dataset, together with comparisons with 2 widely used MSA softwares.

Introduction

Multiple sequence alignment (MSA) is one of the most fundamental tasks in bioinformatics. While there are many attempts to handle comparative sequence analyses without relying on MSA, it still represents a starting point for most evolutionary biology methods. Pairwise sequence alignment has been conceptualized as early as the 1970's, starting with global alignments that attempt to align entire sequences [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF] and then introducing a decade later local alignments that focus on the identification of subsequences sharing high similarity [START_REF] Smith | Comparative biosequence metrics[END_REF]. The standard computational formulation of both tasks is to maximize a scoring function obtained as the sum of the score for each aligned pair of residues (nucleotides or amino acids, the highest scores being attributed to pairs of residues with highest similarity), minus some gaps penalties. Since these seminal works, an abundant literature has flourished exploring this topic in many different directions, from the pairwise problem to the more complex task of aligning more than 3 sequences [one of the very first attempts appearing in 15], from exact solutions that scale exponentially with sequence lengths to faster heuristic approaches used in the most common tools, and from the scoring formulation of the 1 alignment problem that requires to choose the scoring parameters to probabilistic formulations in which those parameters are estimated [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF][START_REF] Durbin | Biological sequence analysis: probabilistic models of proteins and nucleic acids[END_REF]. However, manually refined alignments continue to be superior to purely automated methods and there is a continuous effort to improve the accuracy of MSA tools [START_REF] Edgar | Multiple sequence alignment[END_REF]. We refer the reader to the reviews [START_REF] Edgar | Multiple sequence alignment[END_REF][START_REF] Kumar | Multiple sequence alignment: In pursuit of homologous DNA positions[END_REF][START_REF] Notredame | Recent evolutions of multiple sequence alignment algorithms[END_REF][START_REF] Wallace | Multiple sequence alignments[END_REF] for more details on MSA.

Dynamic time warping (DTW) is a general version of the dynamic programing algorithm that solves exactly the pairwise biological sequence alignment problem. It is a well-known and general technique to find an optimal alignment between two given (time-dependent) sequences. In time series analysis, DTW is used for constructing an optimal alignment of two sequences with possible different lengths by stretching or contracting time intervals [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF]. In functional data analysis, the time warping approach consists in modeling a set of curves exhibiting time and amplitude variation with respect to a common continuous process [START_REF] Liu | Functional convex averaging and synchronization for timewarped random curves[END_REF]. Thus, time warping techniques are used in many different areas concerned by sequence or curve comparisons, one of its most famous successes being on human-speech recognition [START_REF]An overview of sequence comparison: Time warps, string edits, and macromolecules[END_REF].

Here, we propose a simple and fast procedure for MSA, inspired from recent techniques of curve synchronization developed in the context of functional data analysis [START_REF] Arribas-Gil | Pairwise dynamic time warping for event data[END_REF][START_REF] Tang | Pairwise curve synchronization for functional data[END_REF]. In this setup, one often observes a set of curves which are modeled as the composition of an 'amplitude process' governing their common behavior, and a 'warping process' inducing time distortion among the individuals. Specifically, y i (t) = x i •h i (t), t ∈ [a, b], i = 1, . . . , K, are observed, with x i being i.i.d. realisations of the amplitude process X, and h i strictly monotone functions such that h i (a) = a and h i (b) = b being i.i.d. realisations of the warping process H. Aligning pairs of curves (that is eliminating time variation, which comes to estimating the warping functions h i ) is a first step before estimating the common amplitude process. These authors proposed to first estimate pairwise warping functions between all possible trajectories pairs which are then used to create estimators of the underlying individual warping functions in a second step. Sample means or more robust median-based estimates come into play to solve this second step. This procedure is an alternative to the widely used approach of template registration, that consists in aligning every observed curve to some given template, which should be a good estimate of the mean amplitude process. The drawback of this methodology is that it heavily relies on the choice of the template, which is not straightforward. Now, in the MSA context, the warping process is the insertion-deletion (or indel) process that stretches or contracts the initial sequence, while the amplitude process is the substitution process that modifies the value of the sequence base. The equivalent of template registration in the MSA context would be the alignment of every sequence to some estimate of the ancestral sequence, which is, of course, not available. However, exploiting the previous ideas, we show how pairwise alignments can be combined with a simple median-based approach to obtain an estimate of the multiple alignment of the sequences.

Our aim is to establish a proof of concept that our new method could be an interesting ingredient to include into refined MSA techniques. Indeed, the method is able to align a large number K of sequences (that are assumed to share a common ancestor) in a quite simple and fast manner, although a bit rough w.r.t. accuracy. We would like to stress that we do not claim to be competitive with actual aligners. Let us recall that there already exist many competitors to solve the MSA problem whose respective performances have been largely explored [see for e.g 21, and the references therein]. Here, we would rather like to point out to recent developments from curve synchronization that could open the way to new improvements in MSA. While we do not pretend to propose a most accurate method, it is important to note that our approach could be used as a starting point in an iterative refinement strategy [START_REF] Gotoh | Significant improvement in accuracy of multiple protein sequence alignments by iterative refinement as assessed by reference to structural alignments[END_REF]. Those strategies, included in many widely used tools such as ProbCons [START_REF] Do | ProbCons: Probabilistic consistency-based multiple sequence alignment[END_REF], MUSCLE [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF], MAFFT [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: Improvements in performance and usability[END_REF] or MUMMALS [START_REF] Pei | MUMMALS: multiple sequence alignment improved by using hidden markov models with local structural information[END_REF], mostly consist in repeatedly dividing the set of aligned sequences into two random groups and realign those groups by optimizing an objective function. Thus, our simple and fast procedure could be combined with similar refinement strategies that would improve its efficiency. An advantage of our method is that it only uses pairwise comparison but (unlike progressive aligners) is not sensitive to any order on which we consider the sequence pairs. Indeed, progressive aligners rely on a guiding tree and progressively build K -1 pairwise alignments from the union of the initial set of K sequences with consensus sequences built at the internal nodes of the tree. As such, the consensus sequences heavily depend on the guiding tree and the sequences ordering which impacts the resulting MSA. Moreover, progressive aligners are known to tend to propagate errors appearing at the early stages of the method. For instance, ClustalW [START_REF] Thompson | CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, positionspecific gap penalties and weight matrix choice[END_REF] is criticized for being responsible of the heuristic rule 'once a gap, always a gap' that makes errors in pairwise alignments to propagate through MSA. The method proposed here is not sensitive to that issue because a gap in an alignment pair will not necessarily be selected for the MSA.

The manuscript is organized as follows. Section 2 first recalls concepts coming from the pairwise alignment and time warping problems (Section 2.1), sets our framework (Section 2.2) and describes our procedure (Section 2.3). Then Section 3 presents our experiments, starting with synthetic datasets (Section 3.1), then relying on a benchmark dataset (Section 3.2), namely the Balibase dataset [START_REF] Thompson | BAliBASE 3.0: Latest developments of the multiple sequence alignment benchmark[END_REF] and concludes with possible extensions of our approach (Section 3.3).

Framework and procedure

In what follows, we first introduce some ideas around pairwise alignment and the time warping problem. We then explain how to use these concepts in the context of multiple sequence alignment.

Pairwise alignment and time warping

For any two sequences

X = X 1 • • • X n and Y = Y 1 • • • Y m with
values in a finite alphabet A, any (global) pairwise alignment of X and Y corresponds to an increasing path in the grid [0, n] × [0, m] composed by three elementary steps {(1, 0), (0, 1), (1, 1)}, as shown in Figure 1. Note that for biological reasons, such path is often restricted to never contain two consecutive steps in {(1, 0), (0, 1)} (a gap in one sequence may not be followed by a gap in the other sequence). We do not use that constraint in what follows but a postprocessing step could be applied to any alignment in order to satisfy it (the simplest way is then to replace those 2 successive gaps with a (mis)-match move; from a graphical point of view a horizontal+vertical or vertical+horizontal move is replaced by a diagonal one). For notational convenience, we extend the path to [-, n] × [-, m] for some small > 0. Now, any such path may be viewed as an increasing function φ X,Y : [-, n] → [-, m] with the convention that it is càd-làg (continuous on the right, limit on the left). By letting φ(u-) denote the left limit of φ at u, we moreover impose that φ X,Y (0-) = 0 and φ X,Y (n) = m. Note that we can define a unique generalized inverse function φ Y,X :

= φ -1 X,Y : [-, m] → [-, n] constrained to be càd-làg and such that φ X,Y •φ -1 X,Y = φ -1 X,Y
•φ X,Y = Id (the identity function). From a graphical point of view, the path corresponding to φ -1 X,Y is obtained as the symmetric of the path corresponding to φ X,Y with respect to the diagonal line y = x (see Figure 1).

Functions φ X,Y , φ -1 X,Y may be viewed as time warping functions that describe the homology (and thus also the indels) between sequences X, Y obtained from their alignment. Let us first explain this idea on a simple example. Consider an alignment between X = ACAGT AGT and Y = CT T AAG given as follows

A C A G T A -G T -C T -T A A G -
(This alignment corresponds exactly to the thick and blue path depicted in Figure 1). In this alignment, any character from one sequence may be associated to a unique character from the other sequence as follows: when the position corresponds to a match or a mismatch, the character is associated to the (mis)-matching character (in the other sequence) while at positions corresponding to indels, we (arbitrarily) associate it to the character in the previous (mis)-matching position in the other sequence. In order to deal with the case where no such previous (mis)-matching position exists (namely when the alignment begins with an indel), we introduce extra artificial characters (X 0 , Y 0 ) at the beginning of each sequence, these two artificial characters being aligned with each other. In other words, we consider the alignment

X 0 A C A G T A -G T Y 0 -C T -T A A G -
The bare alignment (the alignment without specification of nucleotides) of these two sequences is as follows

• B B B B B B -B B • -B B -B B B B -
Here B stands for base and we do not specify to which character it corresponds. The bare alignment naturally appears to separate the indel process from the substitution process [see for e.g. [START_REF] Arribas-Gil | Statistical alignment with a sequence evolution model allowing rate heterogeneity along the sequence[END_REF][START_REF] Thorne | An evolutionary model for maximum likelihood alignment of DNA sequences[END_REF]. Given the sequences, the knowledge of the bare alignment is sufficient to recover the corresponding full alignment. In this context, the time warping process is exactly the indel process that stretches or contracts an initial sequence, while the amplitude process is the substitution process that modifies the value of the sequence. We focus on the time warping process encoded in the bare alignment.

Let → denote this (asymmetric) relation induced by the bare alignment. In our example, we have that X 1 is associated to the artificial character Y 0 and we denote X 1 → Y 0 . In the same way,

X 2 → Y 1 , X 3 → Y 2 , X 4 → Y 2 , X 5 → Y 3 , X 6 → Y 4 , X 7 → Y 6 and X 8 → Y 6 .
Note that no position of X is associated to Y 5 (this is because of a discontinuity of the function φ X,Y at u = 6 where φ X,Y (6-) = 4 and φ X,Y (6) = 5). More importantly, these associations are exactly encoded in the mapping φ X,Y as we have

X u → Y v if and only if φ X,Y (u-) = v.
More generally, for any real value u ∈ [0, n], letting u denote the smallest integer larger or equal to u, we set X(u) := X u (and similarly for sequence Y ). From the (bare) pairwise alignment of X, Y , we obtain a time warping function φ X,Y such that this alignment is entirely described by the association

X(u) → Y • φ X,Y (u).
In what follows, we denote this functional association by X → Y • φ X,Y . Note that we equivalently have

Y → X • φ -1 X,Y = X • φ Y,X .
We also mention that when the pairwise alignment between X, Y is extracted from a multiple sequence alignment containing at least 3 sequences, say X, Y, Z, the above relation → should be associative in the sense that whenever X → Y • φ and Y → Z • φ , we should also have X → Z • φ • φ. This property will be used in the next section.

Multiple sequence alignment

In this section, we consider a set of sequences S 1 , . . . , S K with values in a finite alphabet A and respective lengths n 1 , . . . , n K and we assume that they all share some latent ancestor sequence A with values in A and length N .

A multiple sequence alignment of the set of sequences S 1 , . . . , S K is given by the knowledge of homologous positions as well as inserted positions. To fix the ideas, let us consider for example the following multiple alignment of 3 sequences

(H) (H) (H) (H) X 0 A C A G T A - - T Y 0 - C T - C C A G - Z 0 T C A C -- - C T
where the first line indicates homologous positions (H) while the other are inserted positions. Homologous positions describe the characters that derived from the ancestor sequence, so that there are exactly N homologous positions in the alignment (whenever an ancestral position was deleted in all the sequences, this can not be reconstructed and such position would not appear in our reconstructed ancestral sequence). For each homologous position, there is at most one character in sequence S i that is associated to it. This means that homologous columns in the multiple sequence alignment may contain matches, mismatches and even gaps (when the corresponding ancestral position has been deleted in that particular sequence). Between two consecutive homologous positions, each sequence might have a certain number of characters that are inserted. These characters do not derive from an ancestral position. Note that these insert runs are not aligned in the sense that the choice of how to put the letters in the insert regions is arbitrary and most MSA implementations simply left-justify insert regions. Now, given the set of sequences, this multiple alignment may be completely encoded through the knowledge of the homologous positions in the sequences (see Section 2.4).

Our goal is to estimate the alignment of each S i to the ancestor A and thus the global alignment of the set of sequences S 1 , . . . , S K , by relying on the set of pairwise alignments of each S i to all the other sequences S 1 , . . . , S K . To do this, we will implicitly assume that a) the multiple sequence alignment of all the sequences S i is well approximated by the alignment we would obtain from all the sequences {A, S i ; i = 1, . . . , K}; b) all the pairwise alignments of S i , S j are good approximations to the extracted pair alignments from the multiple alignment of all sequences {A, S i ; i = 1, . . . , K}.

First, any sequence S i is derived from the common ancestor sequence A through an evolutionary process that can be encoded in the alignment of these two sequences. This alignment induces a time warping process φ i := φ S i ,A : [-, n i ] → [-, N ] such that we have S i → A • φ i . For the moment, we assume that the warping functions φ i 's are i.i.d. realisations of the same random process Φ. For identifiability reasons, we constrain Φ such that E(Φ) = Id (see for example [START_REF] Arribas-Gil | Pairwise dynamic time warping for event data[END_REF]). Moreover, note that assuming the φ i 's are i.i.d. realisations of the same process boils down to considering a star tree leading from the ancestor sequence to any observed sequence S i . In the case of an ultra-metric tree, we would obtain non independent but identically distributed realisations. We will discuss later how to handle more realistic setups.

Since the ancestral sequence A is not observed, the idea of [START_REF] Arribas-Gil | Pairwise dynamic time warping for event data[END_REF] is to estimate the warping function φ i from all the pairwise alignments of S i with the other sequences. More precisely, when considering the pairwise alignment of S i with S j , we obtain a warping function

g ji = φ -1 j • φ i satisfying S i → A • φ i → S j • φ -1 j • φ i = S j • g ji .
As a consequence, relying on the desired associative property of →, we should have S i → S j • g ji so that we can estimate the warping function g ji from a pairwise alignment of {S i , S j }. Note that since pairwise alignment of {S i , S j } does not depend on the order of the sequences, the proposed estimators satisfy ĝji = ĝ-1 ij . Then, we want to combine these estimates {ĝ ji ; j = 1, . . . , K} to obtain an estimate of φ i . In [START_REF] Tang | Pairwise curve synchronization for functional data[END_REF], the authors use the empirical mean and notice that this makes sense since E(g ji |φ i ) = φ i and under the assumption E(Φ) = Id [see also 4]. Let us recall that these works come from functional data analysis where the functions φ i 's are real-valued functions. Here, we are dealing with discrete sequences and the empirical mean of a set of paths g ji would not give a proper path and thus would not correspond to an alignment of S i to the ancestor A. That is why instead of using empirical mean, we will rather rely on median values. Note that -in general -medians have the advantage of being more robust than mean values.

From an algorithmic perspective, a global alignment of

X = X 1 • • • X n and Y = Y 1 • • • Y m
may be encoded trough the sequence of coordinates Z 1 , . . . , Z T , where T is the length of the alignment, Z 1 = (0, 0), Z T = (n, m) and for any t = 1, . . . , T -1, we have Z t+1 -Z t ∈ {(1, 0), (0, 1), (1, 1)}. We let

Z n,m = {(Z 1 , . . . , Z T ); T ≥ 1, Z 1 = (0, 0), Z T = (n, m), such that Z t+1 -Z t ∈ {(1, 0), (0, 1), (1, 1)} for all t = 1, . . . , T -1}, (1)
be the set of possible alignments between X and Y . The link between the path sequence Z and the warping function φ X,Y is as follows

∃t ≥ 1 such that Z t = (k, l) ⇐⇒ φ X,Y (k-) = l.
Now, considering a median path boils down to computing median coordinates of points Z t as explained in the next section.

The median procedure

In this section, we explain how for each sequence S i , we combine the set of pairwise alignments of S i with S j (j = 1, . . . , K) to obtain an estimate of the warping function φ i . We thus fix a sequence S i of length n i and consider its pairwise alignments with all the other sequences S j (j = 1, . . . , K), including itself (the alignment path is then the identical function from [0, n i ] to [0, n i ]). We first define an estimate of the true number N of homologous positions in our MSA (that is positions that were present in the ancestral sequence A, which are the only positions of A that may be recovered), as a median value (constrained to be an integer number) of the set of lengths {n j ; j = 1, . . . , K}, namely N = median{n j ; j = 1, . . . , K}.

More precisely, as a convention and for the rest of this work, we set the median of a sequence of integers (n 1 , . . . , n K ) with ordered values (n (1) , . . . , n (K) ) as follows. Whenever K is even and the mean value m between the ordered statistics n ([K/2]) and n ([K/2]+1) is not an integer, we (arbitrarily) set the median as m . Now, we aim at defining a path from [0, n i ] to [0, N ] by using all the estimated paths ĝji obtained by aligning {S i , S j }. For any position u = 1, . . . , n i , we consider the character S i u throughout all the pairwise alignments. The character S i u may either be aligned with another character S j v or to a gap after position v in sequence S j , as may be viewed in Figure 2 (left part). As already explained, each pairwise alignment of S i with S j corresponds to a sequence of coordinates Z j = (Z j 1 , . . . , Z j T j ) ∈ Z n i ,n j as defined above, see Equation ( 1). Note that we dropped the index i on which the sequence Z j also depends. Now, for each value u = 1, . . . , n i , there is a unique time point t j u and two integers v j u,1 ≤ v j u,2 such that

Z j t j u = (u -1, v j u,1
) and Z j t j u +1

= (u, v j u,2 ).

In practice, either v j u,2 = v j u,1 + 1 when S i u is aligned to a character or v j u,2 = v j u,1 when S i u is aligned to a gap. We then compute two median points

Z1 (u) = u -1, median j=1,...,K {v j u,1 } and Z2 (u) = u, median j=1,...,K {v j u,2 } , (2)
which gives the part of the median path that concerns character S i u (Figure 2, left part). Note that at this step, we did not take into account the possible vertical steps in the alignments ĝji . We consider this now. Indeed, the concatenation of the positions

( Z1 (1), Z2 (1), Z1 (2), Z2 (2), . . . , Z1 (n i ), Z2 (n i ))
is almost what we wanted, namely the median path that gives our estimate of the alignment between S i and ancestral sequence A. However one should notice that for each value 1 ≤ u ≤ n i -1, we either have Z2 (u) = Z1 (u + 1) or Z2 (u) < Z1 (u + 1).

In the first case ( Z2 (u) = Z1 (u + 1)), one of the two points Z2 (u) or Z1 (u + 1) is redundant and should be removed to obtain a proper path in Z n i , N . Otherwise, vertical movements are included in the median path and we keep both values Z2 (u) and Z1 (u + 1).

Let us note that at the last position u = n i , we have Z2 (n i ) = (n i , median j=1,...,K n j ) = (n i , N ) so that the median path correctly ends at (n i , N ).

MSA from median warping

We now combine the different steps to obtain our global alignment procedure as explained in pseudocode in Algorithm 1. This algorithm outputs a set of homologous positions (Hom) as well as a set of inserted positions (Ins). Let us recall that in each sequence, homologous positions correspond to positions that are aligned with an ancestor position A u (u = 1, . . . , N ). Then two homologous positions may be separated by insert runs which are not aligned. We use a table Ins with K rows and N + 1 columns. Each entry (i, u) with u = 2, . . . , N of this table contains the number of insertions in S i between two homologous positions with respective ancestor positions A u-1 and A u , while entry (i, 1) (resp. entry (i, N + 1)) is the number of insertions before first (resp. after last) homologous position. We also use a table Hom with K rows and N columns. Each entry (i, u) of this table contains the position in sequence S i that is homologous to ancestor position A u . When there is no such position (a deletion in sequence S i ), the entry is set to 0.

For instance, going back to the example of a multiple sequence alignment of Section 2.2, we may describe this alignment with the following tables Note that the information contained in table Ins is redundant with that in table Hom but we use it for convenience. Indeed, it is helpful to output the final alignment whose size (i.e. total number of columns) depends on the number of homologous positions plus the number of inserted positions. Now for each sequence S i with i = 1, . . . , K we reconstruct its alignment with ancestral sequence A as follows. For each position u = 1, . . . , n i , we compute the coordinates Z1 (u), Z2 (u) from Equation [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF]. We update the quantity Path(i) that stores the alignment of sequence S i to ancestor A up to position u. Then we either update the table Ins if an inserted position has been obtained (case Z1 (u) [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF] = Z2 (u) [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF] where X [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF] is the second coordinate of vector X) or the table Hom otherwise.

With respect to the algorithmic complexity of the proposed procedure, we can assume that the K sequences have average length similar to N . Then, the time complexity of the MSA from median warping (once the pairwise alignments are given) is O(N K 2 ) since we need K • N iterations in which we compute (twice) a median over K values (the complexity of the median calculation being O(K) relying for e.g. on the Quickselect algorithm [START_REF] Floyd | Algorithm 489: The algorithm select—for finding the ith smallest of n elements [m1][END_REF]). We provide running time comparisons with MAFFT on the BAliBASE dataset at the end of Section 3.2.

Results

Synthetic experiments

In this section, we propose a simple synthetic experiment in order to assess the performances of our approach. We start with an ancestral sequence A on the set of nucleotides A = {A, C, G, T } with length N = 100. From this ancestral sequence, we simulate i.i.d. sequences S 1 , . . . , S K as follows. We rely on the simple Thorne-Kishino-Felsenstein [28, hereafter TKF] model that includes an insertion-deletion process and that can be combined with any substitution process. Here we set parameters for TKF λ = µ = 0.03 and use the Jukes-Cantor substitution model [see for e.g. 30] with substitution rate α = 0.1 (all nucleotide frequencies being set to 1/4). That is, we are simulating K nucleotide sequences related by a star tree and whose branch lengths are set to be equal, as done for instance in [START_REF] Arribas-Gil | Parameter estimation in multiple hidden i.i.d. models from biological multiple alignment[END_REF]. We repeat this experiment M = 100 times for each K, and we let the number of sequences vary in {10, 20, 30, 40, 50}.

To obtain multiple alignments with our procedure, we first conduct all the pairwise alignments between pairs of sequences through the Needleman-Wunsch algorithm for global alignment [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF] as implemented in the pairwiseAlignment function of the R Biostrings library [START_REF] Pages | Biostrings: String objects representing biological sequences, and matching algorithms[END_REF]. We set the parameters for the pairwise alignment as: gap opening penalty equal to -10, gap extension penalty equal to -0.5, nucleotide substitution matrix with diagonal values equal to 5 and non-diagonal values equal to -4. These are commonly used as default parameter values for the global alignment of DNA sequences. That is, we do not look at optimising the alignment parameters, but just at showing that our procedure can produce reasonable results under general conditions.

In order to assess the performance of the method, we compute for each alignment two scores that measure its overall concordance with the simulated one, SP (sum of pairs) and TC (total column) scores (see next section for details), and we compare them to the scores obtained by two well established multiple alignment softwares, namely ClustalW [START_REF] Thompson | CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, positionspecific gap penalties and weight matrix choice[END_REF] in its 2.1 version and T-coffee [START_REF] Notredame | T-coffee: a novel method for fast and accurate multiple sequence alignment[END_REF] in its 10.00.r1613 version. Both softwares are run under default parameter values. Finally, we also run our procedure on the reference pairwise alignments, that is, those extracted from the simulated multiple alignment for each pair of sequences. This is done to assess the performance of the method in the best case scenario, as a way to validate the proposed algorithm for combining pairwise alignments.

Results are presented in Figure 3. The first thing to notice is that our method, when used with the reference pairwise alignments, reaches the maximum score possible in almost every case. Although these are not realistic results since reference pairwise alignments are unknown in practice, they confirm the soundness of our method that combines pairwise alignments from median warping. With respect to the other three methods, in all the scenarios T-coffee always provides the best performance. In general, TC scores tend to get worse as the number of sequences increases whereas SP scores are more stable. However, when comparing our procedure (with estimated pairwise alignments) and ClustalW (which is a widely used software), we see that ClustalW is better for a small number of sequences, but from K = 20, the performance of our method is superior.

Benchmark results

MSA methods are typically benchmarked on sets of reference alignments, the most widely used being the BAliBASE [START_REF] Thompson | BAliBASE 3.0: Latest developments of the multiple sequence alignment benchmark[END_REF]. It is a large scale benchmark specifically designed for multiple sequence alignment. It contains test cases based on 3D structural superpositions of protein sequences that are manually refined to ensure the correct alignment of conserved residues. It is organized into several reference sets, designed to represent real multiple alignment problems. Reference 1 contains alignments of equidistant sequences with 2 different levels of conservation (RV11 and RV12); Reference 2 (RV20) contains families aligned with one or more highly divergent 'orphan' sequences; Reference 3 (RV 30) contains divergent subfamilies; Reference 4 (RV 40) contains sequences with large N/C-terminal extensions; and finally Reference 5 (RV 50) contains sequences with large internal insertions. In addition, three separate Reference Sets [START_REF] Durbin | Biological sequence analysis: probabilistic models of proteins and nucleic acids[END_REF][START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF][START_REF] Edgar | Multiple sequence alignment[END_REF], are devoted to the particular problems posed by sequences with transmembrane regions, repeats, and inverted domains [START_REF] Thompson | BAliBASE 3.0: Latest developments of the multiple sequence alignment benchmark[END_REF]. These last 3 sets are rarely included in benchmark analyses and we exclude them from ours. For each reference set among the first five ones, except for RV 40, two versions of the same multiple alignments are provided: one with the original sequences (which is noted BB), and one with a shorter version of the sequences contained in the alignment (which is noted BBS). In general, BBS alignments should be easier to recover than the original ones, since the less alignable parts of the original sequences have been removed.

The BAliBASE dataset comes with a function called baliscore used to assess the quality of a MSA. Two different criteria are used to score an alignment: SP (sum of pairs) and TC (total column) scores [START_REF] Thompson | A comprehensive comparison of multiple sequence alignment programs[END_REF]. These criteria only use core blocks of the alignment. SP score is the percentage of the correctly aligned residue pairs in the alignment. It is defined as the number of correctly aligned residue pairs found in the test alignment divided by the total number of aligned residue pairs in core blocks of the reference alignment. TC score is the percentage of the correct columns in the alignment. It is defined as the number of correctly aligned columns found in the test alignment divided by the total number of aligned columns in core blocks of the reference alignment.

In order to assess the performance of the method on the Balibase data sets, we proceed as for the simulated alignments of the preceding section. The only differences being that we now consider BLOSUM62 as the default substitution matrix for pairwise alignments (since we now have protein sequences) and that for the last two reference sets (RV40 and RV50) we consider 'overlap' pairwise alignments instead of global ones. Indeed, these two sets are characterized by large differences in sequences lengths so it is convenient to allow for gaps at the beginning and the end of pairwise alignments. In this section we also report the results of ClustalW and T-coffee for reference. We refer to [START_REF] Pais | Assessing the efficiency of multiple sequence alignment programs[END_REF] for a full comparison of available methods on Balibase.

Figure 4 shows that while our method is less performant than ClustalW and T-coffee, it still provides reasonable results. These results should be put in perspective with the level of complexity of the method that is exceptionally low (only pairwise alignments are needed as input and a simple median path is computed). Moreover, as it was the case for the simulated alignments, the scores based on the reference pairwise alignments are almost always equal to 1 in all sets for both SP and TC scores. Again, these are not realistic results since the pairwise alignments extracted from the reference multiple alignment might be far from the optimal pairwise alignment between the two sequences for which no extra information on the remaining sequences is available. Nevertheless, these results are encouraging since once again they serve to validate our method from a theoretical point of view in the sense that the median warping approach provides sound results.

Finally, let us give an idea of the running times. We ran both MAFFT [in its version 7 with default values, 11], which is a very fast and widely used tool, and our procedure on the whole BAliBASE dataset. On a computer with Intel Core i7 at frequency 3.5 GHz and 16 Go RAM (on only one run), we obtained the alignments in 1m23.759s with MAFFT and 1m32.932s with our procedure. Though our algorithm runs on R with no optimization of any sort and MAFFT was run through a command line (outside R environment), these running times (which are not exactly comparable) are roughly the same, which confirms that our procedure is indeed very fast.

Conclusions

In this work, we propose a proof of concept that a simple method derived from recent techniques of curve synchronization in the context of functional data analysis could be of potential interest to MSA. Our method is able to align a large number of sequences in a quite simple and fast manner, although a bit rough w.r.t. accuracy. While we do not claim to be competitive with actual aligners, we believe that our procedure could be successfully included (for e.g. as a starting point) in more refined MSA techniques. While it is out of the scope of the current work to provide such a refined tool, our simulations as well as the use of the Balibase dataset establish that the method has good potential, particularly when looking at the almost perfect results obtained by relying on the (unknown) reference pairwise alignments (that is, the pairwise alignments extracted from the MSA). In particular, while our experiments rely here on exact scoring alignment with default parameter values, the method could be combined with more refined pairwise alignments, such as probabilistic methods that automatically select optimal scoring parameters for the sequences at stake.

While the method implicitly assumes that a star ultra-metric tree describes the evolution of the sequences from their common ancestor, it could be improved to gain in robustness w.r.t. this assumption. In particular, let us assume that additional to the sequences one has access to a non ultra-metric guide tree describing the common evolution of these sequences. Then, we propose to weight the sequences in a way inversely proportional to their distance to the root. Let d i denote the evolutionary distance (i.e. branch length) from sequence S i to root and d max = max i=1,...,K d i . We fix > 0 and set the weight values to Then our method could be generalized to the computation of a weighted median path, namely relying on weighted median values instead of simple ones. In such a way, sequences farther from the ancestor will have a lower weight in the MSA. We leave this for later investigation.

(n, m) • The median value is shown in dashed and red. On the right: a set of K = 3 paths (in blue, including diagonal path) standing for the alignment of S i with all other S j (j = 1, . . . , K) and the median path (in dashed and red). q q 0.5 0.6 0.7 0.8 0.9 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 q q q q q q q q MSA from PW Global ClustlalW T-Coffee MSA from PW Ref.

X 1 X n Y 1 Y n
q q q q q q q q q q MSA from PW Global ClustlalW T-Coffee MSA from PW Ref.

q q q q q q q q q q q q q q q q 0.0

w i = 1 -

 1 d i /(d max + ) k=1,...,K 1 -d k /(d max + ) ∈ (0, 1).

Figure 1 :Figure 2 :

 12 Figure 1: Graphical representation of an alignment between two sequences X 1 • • • X n and Y 1 • • • Y m . The blue and thick line corresponds to a càd-làg function φ X,Y whose generalized inverse φ Y,X = φ -1 X,Y is shown in dashed and red.

Figure 3 :

 3 Figure 3: SP and TC scores on the synthetic data sets. Distribution over M =100 simulated multiple alignments for each value of

Algorithm 1: Pseudocode for MSA from median warping // Pairwise alignments Align(i, j); For any pair of sequences {S i , S j }, compute their pairwise alignment Align(i, j).

// Reconstructing alignment Path(i) of each S i with ancestor A; Z2 (0) ← ((0, 0)); Ins ← table of 0's with K rows and N + 1 columns; Hom ← table of 0's with K rows and N columns;

// Update insertion and homologous positions tables; if Z1 (u) [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF] = Z2 (u) [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF] then

Hom[i, Z2 (u) [START_REF] Arribas-Gil | Parameter estimation in pair-hidden Markov models[END_REF]] ← u // Compute the maximal number of insertions before first, between 2 and after last homologous positions, respectively;

// Compute size of multiple alignment; q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q 0.0 MSA from PW Ref. BBS q q q q q q q q MSA from PW Global BB

MSA from PW Ref. BBS q q q q q q q q q q q q q q q q q q q q q MSA from PW Global BB

MSA from PW Ref. BBS q q q q q q q q q q q q q q q q q MSA from PW Global BB MSA from PW Ref. BBS q q q q q q q q MSA from PW Overlap BB ClustlalW BB T-Coffee BB MSA from PW Ref. BB q q q q q q q q 0.0


C Xi

Latnie

Tan Jl

Wall St

Lee Lc
email: lclee@egr.msu.edu

Patient-specific computational modeling of ventricular mechanics in pulmonary arterial hypertension

Patient-specific biventricular computational models associated with a normal subject and a patient diagnosed with pulmonary arterial hypertension (PAH) were developed to investigate the effects of this disease on ventricular mechanics. These models were developed using geometry reconstructed from magnetic resonance (MR) images, and constitutive descriptors of passive and active mechanics. Model parameter values associated with ventricular mechanical properties and myofiber architecture were obtained by fitting the models with the corresponding measured pressure-volume loops and the circumferential strain calculated from the MR images using a hyperelastic warping method. Our results shows that the peak right ventricle (RV) pressure was substantially higher in the PAH patient when compared to the normal (65 mmHg vs. 20 mmHg). Circumferential strain (E cc ) and ejection fraction (EF) were comparatively lower in both the left ventricle (LV) and RV of the PAH patient (LV EF: 39% vs. 66% and RV EF: 18% vs. 64%; LV E cc : -2.1% vs -9.4% and RV E cc -6.8% vs. -13.2%). On the other hand, passive stiffness, contractility and myofiber stress were all found to be substantially increased in the PAH patient in both the RV and the left ventricle (LV). Septum in the PAH patient was also found to possess a smaller curvature than the LV free wall. Simulations using the PAH model with varying RV preload and afterload revealed an approximately linear relationship between the septum curvature and the transseptal pressure gradient at early-diastole and end-systole, respectively. These findings suggest that PAH can induce LV remodeling and measurements of septum curvature may be useful in quantifying the transseptal pressure gradient in PAH patients.

A 3T Philips scanner with ECG gating was used to acquire cine magnetic resonance (MR) images from one PAH patient as well as one gender-and age-matched normal human subject. Left and right ventricular, atria and arterial pressure were also acquired from the PAH patient by right heart catheterization as part of their clinical care. All data were acquired in the National Heart Center of Singapore and all enrolled participants gave written informed consent.

Introduction

Pulmonary arterial hypertension (PAH) is a disease associated with vasoconstriction of the pulmonary arterial vessels that results in an elevated pulmonary arterial pressure. The clinical diagnosis for PAH is an elevated mean pulmonary artery (PA) pressure greater than 25 mmHg with normal wedge pressure [START_REF] Mclaughlin | ACCF/AHA 2009 Expert Consensus Document on Pulmonary Hypertension. A Report of the American College of Cardiology Foundation Task Force on Expert Consensus Documents and the American Heart Association Developed in Collaboration With the American College o[END_REF] . Without treatment, PAH can quickly lead to decompensated right heart failure and death. The current prognosis of PAH remains poor with about 15% mortality within 1 year on modern therapy 2 , and a low 3-year survival rate of 67% [START_REF] Humbert | Survival in incident and prevalent cohorts of patients with pulmonary arterial hypertension[END_REF] .

The elevated afterload in the right ventricle (RV) caused by PAH can lead to changes in ventricular mechanics [START_REF] Vonk Noordegraaf | The role of the right ventricle in pulmonary arterial hypertension[END_REF] . For example, elevated RV pressure (RVP) can lead to abnormal ventricular deformation in the form of a "left ventricular septal bow" (LVSB), which is a leftward motion of the septum into the left ventricle (LV) cavity. Moreover, PAH is also associated with long-term functional, structural and geometrical changes in the RV [START_REF] Jessup | The effect of chronic pulmonary hypertension on left ventricular size, function, and interventricular septal motion[END_REF] . These changes and their impacts on ventricular mechanics are currently not well-understood [START_REF] Borgdorff | Right ventricular failure due to chronic pressure load: What have we learned in animal models since the NIH working group statement?[END_REF] . Although animal models have been used to understand ventricular mechanics associated with RV remodeling in PAH [START_REF] Hill | Structural and Mechanical Adaptations of Right Ventricle Free Wall Myocardium to Pressure Overload[END_REF][START_REF] Wang | Progressive right ventricular functional and structural changes in a mouse model of pulmonary arterial hypertension[END_REF] , similar studies in humans are lacking and most clinical investigations have been confined to quantifying global ventricular mechanics through pressure-volume (PV) loop measurements [START_REF] Mccabe | Right ventricular dysfunction in chronic thromboembolic obstruction of the pulmonary artery: a pressure-volume study using the conductance catheter[END_REF][START_REF] Trip | Clinical relevance of right ventricular diastolic stiffness in pulmonary hypertension[END_REF] .

Image-based computational models are increasingly used to analyze ventricular mechanics in heart diseases (e.g., myocardial infarction [START_REF] Lee | A Novel Method for Quantifying In-Vivo Regional Left Ventricular Myocardial Contractility in the Border Zone of a Myocardial Infarction[END_REF] , left branch bundle block [START_REF] Aguado-Sierra | Patient-specific modeling of dyssynchronous heart failure: A case study[END_REF] , mitral valve regurgitation [START_REF] Baillargeon | Human Cardiac Function Simulator for the Optimal Design of a Novel Annuloplasty Ring with a Sub-valvular Element for Correction of Ischemic Mitral Regurgitation[END_REF][START_REF] Wenk | First finite element model of the left ventricle with mitral valve: insights into ischemic mitral regurgitation[END_REF] ) as well as heart failure treatments [START_REF] Lee | Bioinjection treatment: Effects of post-injection residual stress on left ventricular wall stress[END_REF] . The focus of these computational analyses is, however, on the left ventricle (LV) that has historically received more attention than the RV. On the other hand, computational modeling analyses of PAH have mostly been performed on pulmonary arteries [START_REF] Hunter | Computational simulation of the pulmonary arteries and its role in the study of pediatric pulmonary hypertension[END_REF][START_REF] Kheyfets | Patient-specific computational modeling of blood flow in the pulmonary arterial circulation[END_REF] while the RV, whose function is a major determinant of prognosis in PAH [START_REF] Voelkel | Right ventricular function and failure: Report of a National Heart, Lung, and Blood Institute working group on cellular and molecular mechanisms of right heart failure[END_REF] , has largely been overlooked.

To fill the void in the current research of PAH, we describe in here an image-based computational analysis of the ventricular mechanics in a PAH patient. Specifically, we seek to quantify regional ventricular myofiber stresses, myofiber strain, contractility and passive tissue stiffness in a PAH patient, and compare them to those found in an age-and gendermatched healthy human subject. We also use the PAH computational model to investigate and quantify the effects of interventricular pressure gradients on the septum curvature.

Method

Acquisition

Biventricular geometry

Three-dimensional biventricular geometries of the PAH patient and normal subject were reconstructed by manually delineating the contours of the LV endocardium, RV endocardium and epicardium in different short and long-axis views in the cine MR images. These contours were then fitted to obtain 3D surfaces of the LV endocardial, RV endocardial and epicardial surfaces using the medical images analysis software MeVisLab (http://www.mevislab.de). Finite element (FE) meshes for the two cases were generated in the volume enclosed by these surfaces using GMSH [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] . The meshes consist of approximately 4000 quadratic tetrahedral elements and 7000 nodes.

Pressure volume loops

Cavity volumes of the LV and RV were measured in the PAH patient and normal subject at different time frames within a cardiac cycle using the MR-reconstructed biventricular geometries described in the previous section. These measured volumes were used together with the acquired pressure data to reconstruct the LV and RV PV loops of the PAH patient. For the normal subject, normal right ventricular pressure waveform from a previous study [START_REF] Redington | Characterisation of the normal right ventricular pressure-volume relation by biplane angiography and simultaneous micromanometer pressure measurements[END_REF] and a scaled normal left ventricular pressure waveform (with end-systolic pressure equal to 0.9 of the measured cuff pressure 21 ) were used together with the cavity volume measurements to reconstruct the PV loops.

Strain Analysis of MR images

Regional circumferential strain was estimated from the PAH patient and the normal subject using hyperelastic warping in FEBio (http://www.febio.org) with the FEWarp plugin. Briefly, in hyperelastic warping, a FE model of the biventricular geometry from the template image was deformed into alignment with the corresponding object in the target image via a computed image-based local body force term that depends on (1) the difference in image intensity between the template and target images, (2) the target image intensity gradient and (3) a prescribed penalty factor [START_REF] Veress | Measurement of strain in the left ventricle during diastole with cine-MRI and deformable image registration[END_REF] . In our analysis, the penalty factor was increased incrementally until the strain field converged. Spatially-averaged circumferential strain was computed in the LV and the RVFW at different cardiac time points using end-diastole as the reference configuration. Circumferential direction was prescribed using the Laplace-Dirichlet Rule-Based algorithm [START_REF] Bayer | A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models[END_REF] with myofiber angle set to be zero.

Biventricular Finite Element Model

Biventricular geometries from the corresponding MR images acquired at early diastole, in which the LV and RV pressures were at their lowest, were used as the reference geometries for the biventricular FE models of the PAH patient and normal subject (Figure 1a). Myocardial fiber orientation was prescribed in the models using a Laplace-Dirichlet Rule-Based algorithm [START_REF] Bayer | A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models[END_REF] . Based previous experimental measurements [START_REF] Streeter | Fiber orientation in the canine left ventricle during diastole and systole[END_REF] , myofiber helix angle was prescribed in the LV with a linear transmural variation from 60°at the endocardium to -60°at the epicardium (i.e., 60°/-60°) (Figure 1b). Given the lack of human RV myofiber data, myofiber orientation in the RVFW was described similarly using a linear transmural variation β° -β°, in which β is a model parameter found from fitting the clinical measurements.

Cardiac Mechanics

A cellular-based cardiac electromechanics model [START_REF] Sundnes | Improved discretisation and linearisation of active tension in strongly coupled cardiac electro-mechanics simulations[END_REF][START_REF] Lee | An Integrated Electromechanical-Reversible Growth Heart Model for Simulating Cardiac Therapies[END_REF] was used to model the mechanical behavior of the biventricular unit over a cardiac cycle. Because the focus here is on biventricular mechanics, we have simplified the model by assuming a homogeneous activation in the biventricular unit. The resultant model can therefore be described by the following system of ordinary differential equations (ODEs) and partial differential equations (PDEs):

∂s ∂t = f (v, s, λ) (1a)
∂v ∂t + I ion (v, s, λ) = I s (1b) ∇ ⋅ σ = 0 (1c)
Briefly, Eq. ( 1a) and (1b) consist of a system of ODEs that describe the local coupling between cellular electrophysiology 27 and cross-bridge cycling process 28 . Here, s, v, I ion and λ denote a vector of state variables consisting of various membrane channels and intracellular ionic concentration, the total ionic current that is scaled with the membrane capacitance, the transmembrane potential and the myofiber stretch, respectively. A stimulus current I s was prescribed over a short time interval to excite the biventricular unit.

Equation (1c) describes the mechanical equlibrium of the cardiac tissue with σ denoting the Cauchy stress tensor. An active stress formulation was used to model the contraction of the tissue. In this formulation, σ was decomposed into a passive component σ p and an active component σ a that accounts for the stiffening of the tissue during cross-bridge cycling, i.e.,

σ = σ p + σ a (s, λ, λ; T ref ) . (2)
The active stress tensor σ a depends on the time evolution of the state variables s, myofiber stretch λ, rate of elastic myofiber stretch λ. The contractile stress was scaled by a reference tension T ref , which a model parameter that is associated with the tissue contractility. The passive stress σ p was described by a Fung-type transversely isotropic hyperelastic constitutive model with the following strain energy function:

Ψ = 1 2 C(e Q -1) + C compr (J ln J -J + 1) (3a)
where

Q = b f f E 2 f f + b xx (E 2 ss + E 2 nn + E 2 sn + E 2 ns ) + b f x (E 2 f n + E 2 nf + E 2 f s + E 2 sf ) (3b)
In the above equation, E = 1 2 (F T F -I) is the Green-Lagrange strain tensor, J is the Jacobian of the deformation gradient tensor F, and E ij with i, j ∈ (f, s, n) are the strain components with f , s, n respectively denoting the fiber, sheet and sheet normal directions. The passive material parameters are C compr , C, b f f , b xx b f x . The passive stress σ p is related to the strain energy function as follows:

σ p = 1 J F ∂Ψ ∂E F T (4)

Simulation of a cardiac cycle

A full cardiac cycle was simulated based on an open-loop circulatory model. Specifically, passive filling was simulated by incrementally applying pressure to the LV and RV endocardial surfaces until the prescribed end-diastolic pressures (EDP) were reached. The biventricular model was then stimulated, and the LV and RV cavity volumes were constrained to remain constant during the isovolumic contraction phase. The ejection phase was simulated by coupling the LV and RV cavity volumes each to a 3-element Windkessel model when the LV and RV pressures exceeded the aortic and pulmonary valve opening pressures, respectively (Figure 1c). The isovolumic relaxation phase was simulated by constraining the cavity volumes and begins when outflows from the LV and RV became negative.

Model parameterization

Biventricular FE models of the PAH patient and normal subject were each divided into two material regions, namely, the LV and the RVFW with different passive stiffness C and contractility T ref (Figure 1a). Together with the Windkessel model peripheral resistance and the RVFW myofiber helix angle β, these material parameters were manually adjusted to fit the (1) experimentally measured PV loops in the LV and RV, and (2) the regional circumferential strain measurements from hyperelastic warping in each case. All other material parameters were prescribed with the same values in both cases. Parameter values in the active contraction model were similar to those in Rice et al. [START_REF] Rice | Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations[END_REF] whereas those in the passive strain energy function (i.e., C compr , b f f , b f x , and b xx ) were similar to those in Genet et al. [START_REF] Genet | Distribution of normal human left ventricular myofiber stress at end-diastole and end-systole-a target for in silico studies of cardiac procedures[END_REF] .

Septum curvature

The fitted PAH model was used to quantify the effects of varying RV afterload and preload on the regional curvature in the LV endocardial surface. Regional curvature κ was calculated at the mid-ventricular level (half-way between the apex to base) from the biventricular FE model at different time point within a cardiac cycle. Specifically, contour of LV endocardial slice at the mid-ventricular level was parameterized as r(θ) in a polar coordinate system (Figure 1d) and κ was computed at each point in the contour using:

κ(θ) = r 2 + 2r ′2 -rr ′′ (r 2 + r ′2 ) 3 2 , ( 5
)
where r ′ ≡ dr dθ. Septum and LV free wall (LVFW) curvature were computed by averaging κ for all the points that are located in their respective region. 3 Results

Pressure-volume loops

The model-predicted PV loops of the RV and LV were in close agreements with the measurements in both PAH and normal cases (Figure 2). The differences between measurements and the model predictions were within 1%. Peak RV pressure was about 3.5 times higher in the PAH patient than that found in the normal subject (65 mmHg vs. 20 mmHg), whereas peak LV pressure was only slightly higher in the PAH patient (110 mmHg vs. 100 mmHg). End-diastolic volume (EDV) and end-systolic volume (ESV) of the LV and RV were also larger in the PAH patient than in the normal subject with the largest difference found the RV (RV EDV: 275 ml vs. 100 ml; RV ESV: 225 ml vs. 36ml). Correspondingly, the ejection fraction (EF) of the PAH patient was substantially lower when compared to that of the normal subject (LV EF: 39% vs. 66% and RV EF: 18% vs. 64%).

Regional circumferential strain

Circumferential strain measured using hyperelastic warping was substantially depressed in both LV and RVFW of the PAH case when compared to the normal case (Figure 3a,b). Specifically, the spatially-averaged peak circumferential strain in a cardiac cycle was -2.1% and -9.4% in the RVFW of the PAH and the normal cases, respectively. Left ventricular peak circumferential strain was also depressed in the PAH patient when compared to the normal subject (-6.8% vs. -13.2%). In both cases, peak circumferential strain in the LV was found to be larger than that in the RVFW. The model predictions agreed reasonably well with these strain measurements, with maximum differences between the model-predicted peak strain and the measurements at about 5% in the LV of the normal subject . We also found that the circumferential strain is very sensitive to the fiber orientation, and in order to capture the differences in strain between the LV and RVFW, transmural variation of myofiber angle in the RVFW has to be larger than that in the LV.

Myofiber stress

Myofiber stress in the PAH patient was substantially higher than that of the normal subject in the entire biventricular region, with regions of high myofiber stress found in the LVFW and septum (Figure 4). Specifically, peak average myofiber stress in the RVFW and LV of the PAH patient were almost twice the values found in the normal subject (LV: 99.7kPa vs. 51.2kPa; RVFW: 73.2kPa vs. 43kPa). For both of the PAH patient and the normal subject, the stress in the LV is higher than that in the RVFW. Table I shows the model parameters for the PAH patient and the normal subject that were calibrated to against the PV loops and the regional circumferential strains. The model parameters C and T ref reflects the tissue's passive stiffness and contractility, respectively. We found that the tissue's passive stiffness and contractility were both higher in the LV and RV of the PAH patient than those in the normal subject. Peripheral resistance in the RV was also substantially higher in the PAH patient than in the normal subject, but was similar in the LV of both. We also found that the transmural variation of myofiber angle was steeper in the RV with the fibers oriented more longitudinally than in the LV in both cases. This regional variation was found necessary to produce a smaller circumferential strain in the RV (compared to the LV) as found in the hyperelastic warping measurements (Figure 3). Two-dimensional curvature κ calculated from Eq. ( 5) on the LV endocardial surface at the mid-ventricular short-axis slice shows a substantially lower κ in the septum of the PAH patient (Figure 5). Compared to the normal subject model, which has a homogeneous curvature in the LV, the septum curvature in the PAH patient was substantially lower than that found in the LVFW. To eliminate the effects of LV diameter on κ, the septum curvature was normalized against the average LV curvature to give a normalized curvature κ n . We found that the normalized curvature κ n was lower in the septum of the PAH patient than that in the normal subject.

Model parameters

Septum curvature

Effects of afterload

The septum curvature κ decreased with increasing RV afterload (Figure 6a), which was adjusted in the PAH model by simultaneously increasing the pulmonary peripheral resistance R per,p and decreasing the RVFW contractiliy T ref so that the stroke volume was relatively unchanged (Figure 6c). The largest variation of septum κ with RV afterload occurs primarily during the systole phase, with the maximum curvature κ varying almost linearly(r = 0.984) with the end-systolic transseptal pressure gradient between the LV and RV (i.e., LVP -RVP) (Figure 6b). On the other hand, we found that the curvature κ was less than the normal case (Figure 5b) even after reducing the afterload in the PAH model to a value close to the normal physiological RV peak pressure (30 mmHg). The septum curvature also decreased with increasing RV preload (Figure 7a), which was simulated by increasing the prescribed end-diastolic pressure. Due to length-dependent effects of the contraction model, RVFW contractiliy T ref was simultaneously adjusted to maintain a constant stroke volume in each of the different preload cases with respect to the baseline PAH case (Figure 7c). In contrast to the effects of afterload, the largest variation of septum curvature κ occurs primarily during early diastole, whereas the curvature was little changed in systole. Interestingly, the maximum diastolic curvature κ max seems to vary almost linearly (r = 0.997) with the transseptal pressure gradient of between LV and RV at early diastole. Specifically, we found that an increase in transseptal gradient by 1 mm Hg produce an average decrease in the curvature by 0.0156 cm -1 .

Effects of preload

Discussion

We have described the development of patient-specific biventricular computational models using pressure measurements and MR data from a PAH patient and a normal human subject. In each case, we have reconstructed the LV and RV PV loops using the pressure data and cavity volumes quantified from the MR images. A hyperelastic warping method was also used to compute the LV and RVFW circumferential strains from cine MR images. These strain measurements and PV loops were used to fit the biventricular computational models.

Using these fitted models, we have compared differences in the in vivo regional ventricular mechanics between the PAH patient and normal subject, which includes regional myofiber stress distribution, mechanical properties and regional curvature. Additionally, we have also quantify the effects of RV preload and afterload on the septal curvature using the PAH model. As most previous computational studies of PAH have focused on the arteries [START_REF] Hunter | Computational simulation of the pulmonary arteries and its role in the study of pediatric pulmonary hypertension[END_REF][START_REF] Kheyfets | Patient-specific computational modeling of blood flow in the pulmonary arterial circulation[END_REF] , this study, to the best of our knowledge, is the first patient-specific computational analysis of the ventricular mechanics associated with PAH.

Circumferential Strain

Compared to the LV, quantitificaion of RV strain is significantly more challenging because of its complex geometry and thin wall. Nevertheless, RV strain analysis on PAH patients have been performed using speckle tracking on echocardiographic images [START_REF] Smith | Three-dimensional speckle tracking of the right ventricle: Toward optimal quantification of right ventricular dysfunction in pulmonary hypertension[END_REF] as well as strainencoded MR imaging [START_REF] Oyama-Manabe | The strain-encoded (SENC) MR imaging for detection of global right ventricular dysfunction in pulmonary hypertension[END_REF] . In these analyses, RV circumferential strain of PAH patients was found to be significantly depressed compared to normal subjects. Our results, which were obtained using a hyperelastic warping method [START_REF] Veress | Measurement of strain in the left ventricle during diastole with cine-MRI and deformable image registration[END_REF] on cine MR images, agree with these studies.

Besides a reduction in RV circumferential strain, we also found that the LV circumferential strain was reduced in the PAH patient compared to the normal subject (-6.8% in PAH vs. -13.2% in Normal). This feature was also found in another echocardiography study [START_REF] Puwanant | Ventricular geometry, strain, and rotational mechanics in pulmonary hypertension[END_REF] in which both the LV lateral wall and septum circumferential strains were found to be significantly reduced by about 5% in pulmonary hypertensive patients.

Myofiber Stress

The use of Laplace's law to estimate RV wall stress is complicated by its irregular geometry [START_REF] Bellofiore | Methods for measuring right ventricular function and hemodynamic coupling with the pulmonary vasculature[END_REF] . Here, we have overcome this complication in our computational analysis, and have found a higher myofiber stress in the RVFW of the PAH patient that is about twice as high as that in the normal subject. The larger RVFW myofiber stress can be attributed to the combination of a higher afterload (3× of the normal) and a substantially bigger RV (2× of the normal) in the PAH patient. In addition, we found that the LV myofiber stress was also about twice as high as that in the normal subject despite having a normal LV afterload. The higher LV myofiber stress could be attributed to the larger LV cavity (1.5 × of the normal). Given that elevated myofiber stress is feature that is associated with ventricular remodeling, this result suggests that the remodeling associated with PAH may not only be confined to the RV, but could also be found in the LV.

Fitted model parameters Ventricular Properties

The fitted model parameters associated with the regional biventricular mechanical properties also suggest significant remodeling in the PAH patient. Specifically, we found that the diastolic passive stiffness C and the contractility T ref of the PAH patient were larger than that of the normal subject in both LV and RVFW. These findings were consistent with previous experimental studies. By comparing the passive tension found in the RV tissue of PAH patients undergoing heart/lung transplantation and donors with non-failing heart, Rain et al. [START_REF] Rain | Right Ventricular Diastolic Impairment in Patients With Pulmonary Arterial Hypertension[END_REF] found that RV diastolic stiffness is increased in PAH patients and is closely associated with the disease severity measured using stroke volume and 6-minute walk test. They attributed this increase to RV fibrosis and the intrinsic stiffening of the RV myocytes. A similar finding was also reported in rats subjected to pulmonary artery banding [START_REF] Hill | Structural and Mechanical Adaptations of Right Ventricle Free Wall Myocardium to Pressure Overload[END_REF] . Interestingly, the study by Rain et al. [START_REF] Rain | Right Ventricular Diastolic Impairment in Patients With Pulmonary Arterial Hypertension[END_REF] also found an increase in the RV end-systolic elastance and the force-generating capacity in isolated RV cardiomyocytes from PAH patients. An increase in RV end-systolic elastance was also found in rats induced with PH via monocrotaline injection compared to the controls [START_REF] De Man | Bisoprolol delays progression towards right heart failure in experimental pulmonary hypertension[END_REF] . As suggested in these studies, the increase in RV contractility associated with PAH is most likely a compensatory mechanism in the RV attempting to cope with the increased afterload.

That we found that the LV diastolic passive stiffness and contractility was also increased suggests the presence of remodeling in the LV of the PAH patient (as with our findings for circumferential strain and myofiber stress). Although we did not find any studies that measure LV mechanical properties directly associated PAH, it was observed that patients with right heart failure often have relaxation abnormalities that may affect the mechanical properties of the LV 36 . Moreover, recent studies have also found evidence of electrophysiological remodeling in the LV of PAH patients 37 .

Myofiber orientation

Histological studies have found that compared to the LV that has muscle fiber running obliquely towards the apex, the deep muscle fibers of the RV are longitudinally aligned from base to apex 38 . Our result also suggests that the myofiber orientation is more longitudinal in the RV than in the LV in both cases, an outcome that is in part due to the smaller circumferential strain found in the RVFW in both cases. On the other hand, studies on rats 7 have found that both myofiber and collagen orientations in the RVFW became more longitudinal as a result of RV remodeling caused by PAH. We, however, did not find this feature in our study.

Septum Curvature

The fact that abnormal septal geometry and motion is a well-recognized feature in PAH patients [START_REF] Marcus | Impaired left ventricular filling due to right ventricular pressure overload in primary pulmonary hypertension: noninvasive monitoring using MRI[END_REF] has motivated us to quantify this feature, as well as investigate its sensitivity to RV preload and afterload. Our results show that the septum curvature from the short axis view at mid-ventricular level was substantially lower in the PAH patient compared to that in the normal subject, which implies a more "flattened" septum.

With varying RV preload (and constant RV afterload), we found that the septum curvature at early-diastole increased linearly with increasing diastolic transseptal pressure gradient (LVP-RVP) (r = 0.997). Septum curvature in systole is not much affected by changes in the preload. Similar findings were also reported in other studies [START_REF] Tanaka | Diastolic bulging of the interventricular septum toward the left ventricle. An echocardiographic manifestation of negative interventricular pressure gradient between left and right ventricles during diastole[END_REF][START_REF] Beyar | Ventricular interaction and septal deformation: a model compared with experimental data[END_REF] , which found that the septum shape and motion at early-diastole are sensitive to the transseptal pressure gradient between the LV and RV. Tanaka et al. [START_REF] Tanaka | Diastolic bulging of the interventricular septum toward the left ventricle. An echocardiographic manifestation of negative interventricular pressure gradient between left and right ventricles during diastole[END_REF] found that diastolic "bowing" of septum toward the LV (negative curvature) in PAH patients resulted from a negative transseptal pressure gradient whereas Beyar et al 41 reported that "bowing" occurs only at a transseptal pressure gradient lower than -5 mmHg. Although we did not observe any "bowing" for the range of transseptal pressure gradient between -1.5 to 2.5 mmHg, our finding of a linear variation between the pressure gradient and septum curvature is consistent with these studies.

Similarly, we found that the end-systolic septum curvature increased almost linearly with increasing systolic transseptal pressure gradient (r = 0.984) while the curvature at diastole is insensitive to the variation in afterload. Our result is quantitatively in good agreement with an echocardiographic study on children with pulmonary hypertension 42 , which found a correlation of r = 0.86 between septum curvature and RV end-systolic pressure. This direct relationship between septum curvature and end-systolic pressure would suggest that the curvature can potentially be used as a non-invasive marker for quantifying RV pressure in PAH patient.

Critique of the model

Overall, predictions of the biventricular computational models are generally consistent with the results from clinical and experimental studies of PAH. Our results have also suggest that PAH may produce significant LV remodeling though caution must be exercised in extrapolating these results, which were based only from one model. Additionally, there are also some limitations that are associated with our models.

First, we have used the biventricular geometry that corresponds to the lowest pressure at the filling phase as the zero-stress state, which may affect the values of the fitted mechanical properties. This, of course, is an approximation as the zero-stress state does not exist in vivo. Other investigators have used the end-systole [START_REF] Walker | MRI-based finite-element analysis of left ventricular aneurysm[END_REF] or mid-systole [START_REF] Kerckhoffs | Patient-specific modeling of the cardiovascular system: Technology-driven personalized medicine[END_REF] as the zero-stress state, or have estimated the zero-stress state from the end-diastole geometry using prescribed mechanical properties 12 . To the best of our knowledge, there appears to be no consensus on the best approach to address this issue.

Second, we have not sampled the entire parameter space in fitting our model. Specifically, we have kept the passive tissue anisotropy fixed based on values obtained in normal humans [START_REF] Genet | Distribution of normal human left ventricular myofiber stress at end-diastole and end-systole-a target for in silico studies of cardiac procedures[END_REF] and the LV myofiber orientation fixed at 60 ○ /-60 ○ . We have done so in order to focus on the key features found in the RV that are associated with PAH.

Figure 1 :

 1 Figure 1: (a) 3D Biventricular geometry reconstructed from MR images with LV and RVFW material regions. (b) Myofiber orientation prescribed using LDRB method with a transmural variation of 60°(endo) to -60°(epi) for the entire model. (c) Coupling the biventricular model to 3-element Windkessel models. (d) Polar coordinate used to calculate local curvature κ in the LV endocardium.

Figure 2 :

 2 Figure 2: Measurements and model prediction of the PV loops for the (a) PAH case and (b) normal case.

Figure 3 :

 3 Figure 3: Circumferential strain E cc for the (a): PAH and (b): normal case. (c): Comparision of regional peak negative E cc between PAH and normal cases.

Figure 4 :

 4 Figure 4: Peak myofiber stress in the PAH patient (left) and the normal subject (right).

Figure 5 :

 5 Figure 5: Curvature κ in the LV endocardial surface of (a) PAH patient and (b) normal subject. (c) Comparison of normalized curvature κ n between PAH patient and normal subject.

Figure 6 :

 6 Figure 6: Effects of afterload on septum curvature within a cardiac cycle. (a): Time-course of κ with different RV afterload. (b): End-systolic curvature κ max,s vs. end-systolic transseptal pressure gradient. (c): Corresponding PV loops.

Figure 7 :

 7 Figure 7: Effects of preload on septum curvature within a cardiac cycle. (a): Time-course of κ with different transseptal pressure gradient. (b): Early-diastolic curvature κ max,d vs. earlydiastolic transseptal pressure gradient. (c): Corresponding PV loops.

Table I :

 I Regional Model Parameters: Passive Stiffness C(kpa), Myocardial Contractility T ref (kpa), Peripheral Resistance R per (mmHg ⋅ min/mL)

	Case	C	LV T ref R per, s	C	RV T ref β	R per, p
	PAH	1.68 262.5 0.005 1.96 160 80 ○ 0.00325
	Normal 0.035 123	0.005 0.12 60 75 ○ 0.000625