Guillaume Pagès

Sergei Grudinin
email: sergei.grudinin@inria.fr

Analytical symmetry detection in protein assemblies. II. Dihedral and Cubic symmetries

Keywords: Point-Group Symmetry, Cubic Groups, Group Theory, Discrete Optimization, Protein Structure, Continuous Optimization

Protein assemblies are often symmetric, as this organization has many advantages compared to individual proteins. Complex protein structures thus very often possess high-order symmetries. Detection and analysis of these symmetries has been a challenging problem and no efficient algorithms have been developed so far. This paper presents the extension of our cyclic symmetry detection method for higher-order symmetries with multiple symmetry axes. These include dihedral and cubic, i.e., tetrahedral, octahedral, and icosahedral, groups. Our method assesses the quality of a particular symmetry group and also determines all of its symmetry axes with a machine precision. The method comprises discrete and continuous optimization steps and is applicable to assemblies with multiple chains in the asymmetric subunits or to those with pseudo-symmetry.

We implemented the method in C++ and exhaustively tested it on all 51,358 symmetrical assemblies from the Protein Data Bank (PDB). It allowed us to study structural organization of symmetrical assemblies solved by X-ray crystallography, and also to assess the symmetry annotation in the PDB. For example, in 1.6% of the cases we detected a higher symmetry group compared to the PDB annotation, and we also detected several cases with incorrect annotation. The method is available at http://team.inria.fr/nanod/software/ananas. The graphical user interface of the method built for the SAMSON platform is available at http://samson-connect.net.

Introduction

Symmetrical protein complexes are very common in nature, as has been highlighted in our previous work on symmetry detection in cyclic protein assemblies [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF], and many of these are deposited to the Protein Data Bank (PDB) [START_REF] Rose | The rcsb protein data bank: integrative view of protein, gene and 3d structural information[END_REF]. As function of proteins is very often determined by their structure, it appears that complex function requires complex structures [START_REF] Levy | 3d complex: a structural classification of protein complexes[END_REF][START_REF] Levy | Assembly reflects evolution of protein complexes[END_REF]. High-order symmetries are thus essential to build large and complex protein assemblies. Dihedral and cubic groups are overrepresented among large protein assemblies with some specific structural functions, for example those of viral capsids. Also, high-order symmetry drastically reduces the complexity of de novo design of self-assembling nanomaterials [START_REF] King | Computational design of self-assembling protein nanomaterials with atomic level accuracy[END_REF][START_REF] King | Accurate design of co-assembling multicomponent protein nanomaterials[END_REF][START_REF] Bale | Accurate design of megadalton-scale two-component icosahedral protein complexes[END_REF][START_REF] Hsia | Design of a hyperstable 60-subunit protein icosahedron[END_REF].

To assess the quality of symmetry for such assemblies, a cyclic symmetry measure is necessary, as the cyclic axes constitute the basic bricks from which one can reconstruct high-order symmetry groups. However, considering each symmetry axis separately would result in a globally incorrect assessment, as there are strict geometrical constraints between different axes of symmetry in high-order symmetry groups. This motivated us to develop a symmetry detection method specifically suited for dihedral and cubic groups. Indeed, the need for this symmetry detection method exists, as some approximate methods, i.e. those from BioJava [START_REF] Prlić | Biojava: an opensource framework for bioinformatics in 2012[END_REF], are massively used to display the symmetry axes on the PDB website [START_REF] Rose | The rcsb protein data bank: integrative view of protein, gene and 3d structural information[END_REF].

Inspired by the quaternion arithmetic applied to the best superposition of a set of points [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF][START_REF] Diamond | A note on the rotational superposition problem[END_REF][START_REF] Kearsley | On the orthogonal transformation used for structural comparisons[END_REF] together with our recent developments [START_REF] Popov | Rapid determination of rmsds corresponding to macromolecular rigid body motions[END_REF][START_REF] Neveu | RapidRMSD : Rapid determination of RMSDs corresponding to motions of flexible molecules[END_REF][START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF], here we propose a new symmetry measure and an analytical method to find the best symmetry axes of a symmetrical assembly possessing multiple symmetry axes. The method guaranties that the detected axes are consistent with the symmetry constraints. Similar to the case of cyclic symmetry detection [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF], our method produces results with a machine precision, its cost function is solely based on 3D Euclidean geometry, and most of the operations are performed analytically. This makes it extremely fast and particularly suitable for exhaustive analysis of PDB data. Below we provide details about the high-order symmetry measure and the computation of the symmetry axes for an assembly possessing any point symmetry group. The method first perceives the topology between different chains, and is able to deal with complex subunits that are composed of multiple chains. Then it iteratively solves a constrained quadratic optimization problem using a set of analytical solutions.

Methods

Notations

Similarly to the analysis of the cyclic groups [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF], in this paper we will be mainly dealing with 3 × 3 matrices and 3-vectors. Therefore, bold upper case letters (i.e. A) will denote matrices, bold lower case letters (i.e. b) will denote vectors, and normal weight lower case letters (i.e. c) will denote scalars. For trigonometric operations and illustrations we will also use an arrow notation for 3-vectors, such as v.

All amino acids, except glycine, are chiral. Hence, symmetry groups that can be present in protein assemblies cannot contain any reflection, inversion, or improper rotation. The only remaining finite point groups are the cyclic (C n for the cyclic group of order n), dihedral (D n for the dihedral group of order n), tetrahedral, octahedral and icosahedral (respectively T , O and I), the three cubic groups. Symbol Γ ∈ {C n } n>1 {D n } n>1 {T, O, I} will denote one of these point groups. Its cardinality, i.e. the number of its elements, will be denoted as |Γ|.

Root mean square deviation

As in our previous work on cyclic groups [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF], we will express the symmetry measure of a molecular assembly using the root mean square deviation (RMSD). Given two sets of N points each, A = {a i } N and B = {b i } N , the RMSD between them is given as

RMSD(A, B) 2 = 1 N 1≤i≤N |a i -b i | 2 .
(1)

Group Theory

Firstly, we will give a brief introduction to the group theory used in the paper. A group is a set of elements equipped with an operation that combines any two elements to form a third element. Formally, it can be written in the form of (Γ, ⋆), where Γ is a set of elements supplied with a group operation ⋆. A homomorphism is a function that takes a group element as input and returns an element of another group as output, preserving the group structure. Given two groups (Γ 1 , ⋆) and (Γ 2 , * ), a homomorphism f from Γ 1 to Γ 2 satisfies the following property,

∀g, g ′ ∈ Γ 1 : f (g ⋆ g ′ ) = f (g) * f (g ′ ).
(

) 2
A homomorphism is called bijective if

∀g 2 ∈ Γ 2 ∃!g 1 ∈ Γ 1 such that f (g 1 ) = g 2 , (3)
meaning that it provides a one-to-one association between the elements of the two groups. In this paper, we will consider three different types of groups :

• Rotation groups, where the elements are rotations and the group operation is the composition of rotations. The group of all rotations in 3 dimensions will be noted SE(3), and rotations will be noted r.

• Permutation groups, where the elements are permutations and the group operation is the composition of permutations. A permutation σ of the set {1, .., n} will be noted (σ(1), ..., σ(n)). The group of all permutations of n elements will be noted S n .

• Point groups Γ that will be composed of abstract elements. We know how these elements are combined (meaning that we know the results of composition of any two elements), and we will take the freedom to use the elements either as permutations or as rotations, thanks to the bijective homomorphisms that will be computed.

The considered point groups are only those that can be described with rotation operators (no reflections or inversions). We say that a set of points (a protein assembly will be represented with a set of points) has a Γ symmetry if the rotations in Γ keep this set globally invariant. More precisely, any rotation operator in Γ will displace each point that is located outside of the rotation axis, but another point will take its place. It is natural to see Γ as a rotation group, but also as a permutation group, since a rotation applied to a set of points will permute them. For example, Figure 1 illustrates a C 5 point group, which has a bijective homomorphism with a rotation group and a bijective homomorphism with a permutation group.

Problem definition

For a molecular point group Γ ∈ {C n } n>1 {D n } n>1 {T, O, I}, and an assembly A = {a i,j } Ns,Na consisting of N s subunits, each composed of N a atoms, we want to minimize the following loss function,

Loss 2 = 1 |Γ|N s N a min σg,rg g∈Γ Ns i=1 Na j=1 a σg(i),j -r g (a i,j ) 2 , (4)
such that g → σ g and g → r g are the bijective homomorphisms from Γ to subsets of S n and SE(3). The loss function is the sum of RMSDs between the original assembly and the rotated assemblies for every rotation in the group Γ. We should mention that this loss function is very natural, since it is only based on Euclidean 3D distances, no adjustable parameters are required and all the rotations r g have equal importance.

Workflow

Minimization of the loss function 4 requires optimization over the group of rotations, which is a continuous optimization, and over the group of permutations, which is a discrete optimization. In practice, we do not know how to do both simultaneously, so we first apply a heuristic approach to determine the correspondences σ g between the subunits, and then we optimize the rotations. Overall, we solve the optimization problem in three steps, 1. Subunit definition 2. Estimation of the permutations 3. Optimization of the rotations [START_REF] King | Computational design of self-assembling protein nanomaterials with atomic level accuracy[END_REF] The third step is an analytical continuous minimization. It gives the expected result with a machine precision, as we have already demonstrated for cyclic symmetries [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF]. The first two steps are heuristics that assume the assembly to be symmetric enough, which allows to estimate the best correspondences between the subunits. This problem of estimating the best correspondence is discrete. Any error during the estimation of the correspondences typically leads to the solution of optimization problem 4 with the Loss comparable to the distance between the center of masses (COMs) of the subunits. It is therefore straightforward to verify whether the result of our discrete optimization is correct. Without loss of generality, we assume that the molecular point group called Γ is given. If it is not the case, we can exhaustively search over three cubic groups and also try all the dihedral groups below a certain maximum order. Below we discuss the individual steps of our optimization algorithms in more detail.

Finding the subunits

We call subunit a minimum part of the symmetrical assembly, from which the entire assembly can be reconstructed by replicating the subunits according to the symmetry operator. Note that the total number of subunits in a complete symmetrical assembly has to be equal to |Γ|. In most of the practical cases, subunits are actually the individual chains of the molecular assembly. In this case, it is straightforward to define them. In some cases, however, subunits can be composed of several chains. To find these multichain subunits, we create several sets of chains, each forming a Γ symmetrical assembly. Then, the subunits are defined by appropriately choosing one chain from from each of these sets as shown in Fig. 2. All the sets are Γ-symmetric assemblies, and thus contain precisely |Γ| chains. Since we know the initial number of chains in the assembly, we also know that the number of sets has to be equal to the number of chains divided by |Γ|. These sets are constructed using a penalty function defined for a pair of chains with four contributions. These are obtained by computing a pairwise sequence alignment followed by a structural superposition for all the pairs of chains. For two chains i and j in an assembly composed of n chains, we define their penalty function as

d(i, j) = d seq (i, j) d max seq + 40 + d struct (i, j) d max struct + 3 Å + d angle (i, j) d max angle + 0.05rad + d center (i, j) d max center + 3 Å , (6) where
• d seq (i, j) is the minus BLOSUM62 score [START_REF] Henikoff | Amino acid substitution matrices from protein blocks[END_REF] of the sequence alignment between chains i and j, computed using the MUSCLE package [START_REF] Edgar | Muscle: multiple sequence alignment with high accuracy and high throughput[END_REF].

• d struct (i, j) is the RMSD between two chains i and j, after the best superposition of the corresponding alpha-carbons.

• Let us define an angle α(i, j) between two chains i and j as the angle returned by the superposition procedure, i.e. the rotation angle between them after the COMs were superposed. In a perfectly symmetrical assembly, only a few values of these angles are possible. More precisely, if a group Γ contains multiple symmetry axes n j ...n k of orders j...k, correspondingly, then the pairs of chains will be mutually rotated by angles 2mπ/l, where 1 ≤ m < l with l = j...k. We define d angle (i, j) as the absolute value of the difference between α(i, j) and the closest listed angle.

• d center (i, j) is the distance between the initial COM of the whole assembly and its position after applying the rigid-body transformation that superposes chain i with chain j.

In all the contributions in the equation above, d max stands for the largest value of the corresponding contribution. This way, we ensure that all the terms in the above equation have weights of approximately the same magnitude, and none of these terms are bigger than 1. To improve the discrimination of the penalty function for the assemblies with nearly perfect symmetry, we also add constants to d max to make the values of denominators sufficiently large.

After having computed all the pairwise penalties d(i, j), we apply an algorithm to cluster the chains into several sets, as it is listed in Alg. 1. To merge two sets, we define a new set containing all the chains from the two other sets. The clustering algorithm is based on the computation of pairwise distances between the sets, where each distance is defined as an average of the penalties between all pairs of chains, such that one chain belongs to the first set and another chain belongs to the second set.

Once the sets are computed, we proceed with the construction of the subunits. To do so, we first compute for each of the sets the symmetry axes corresponding to the group Γ. The detailed axes computation procedure is explained below. Then, we group chains from different sets to construct the subunits. We first match axes computed for all the sets, then we choose one arbitrary reference chain from each set and assemble them to create the first subunit. For all the other subunits, we choose a group operator and assemble together all the chains from the different sets that correspond to the reference chain to which this operator is applied. The correspondence estimation method is detailed below.

foreach

Correspondences between subunits after a rotation

The number of bijective homomorphisms g → σ g from Γ to subsets of S n grows exponentially with the size of Γ, it is thus not feasible to do an exhaustive search for |Γ| bigger than 10. Therefore, to estimate the correspondence between the subunits we set up heuristics.

Let us start by imagining a perfectly symmetrical system consisting of n subunits, with the center of the symmetry located at the origin. Let {c i } {1≤i≤n} be the COMs of the subunits. Remark that c i are located on a sphere. The convex hull of {c i } {1≤i≤n} , which is a polyhedron with n vertices, possesses the following properties: each axis of symmetry of order s > 2 crosses two faces of the convex hull, which are the regular polygons with s vertices, and each axis of symmetry of order 2 crosses an edge of the convex hull. We first create a reference graph from the edges of the convex hull that belong to the regular polygons with one of symmetry axes of order s > 2 passing through them. The topology of this graph does not depend on the choice of the selected perfectly symmetrical system (see Figure 3). For Examples of finding subunits in a C 5 complex containing 10 chains, which leads to 2 sets and 5 subunits. Each chain is shown as a shape, with a local coordinate frame computed with a structural superposition. A correct clustering in this example should determine a set with the chains 1,2,3,4,5 (light grey) and another one with the chains 6,7,8,9,10 (dark grey). The correspondence between the subunits after a rotation is represented by the colored arrows. Chains belonging to the same subunits are drown with the same outline colors. After we define the two sets, the first subunit is created by taking arbitrary chains from each cluster (1 and 7 in (A), 1 and 6 in (B), 1 and 9 in (C)).

The other subunits are created using the correspondences between the chains in the sets. Note that the assembling of subunits is not unique. In (A) and (B) the subunits seem to be assembled more naturally, as they contain chains that are spatially proximate to each other. However (C) is a perfectly valid assembly too. All the three subunit assemblies give exactly the same result in terms of symmetry axes and the RMSD loss function.

the dihedral assemblies we also include into the graph the edges crossing the 2-fold axes. The main idea behind the correspondence estimation method is to fit our non-perfect system to a perfect canonical example, and then to use the known correspondence from this canonical example to deduce the correspondence of our system. To create such a fitting, the COMs of the subunits are projected on a sphere centered in the COM of the assembly.

Then the polyhedral graph P [START_REF] Grünbaum | Convex polytopes: Prepared by volker kaibel, victor klee, and günter ziegler[END_REF] of the convex hull is computed, and we seek for subgraphs of P that are isomorphic to the reference graph. This problem is generally known as the subgraph isomorphism problem, which has been well studied in literature [START_REF] Ullmann | An algorithm for subgraph isomorphism[END_REF]. Since our example is not perfect, it happens that we obtain either zero or several matches. If we get zero matches, we connect the two most spatially proximate yet unconnected vertices in our graph with an edge and restart the subgraph isomorphism procedure until we obtain some result. If we get several results, we use two geometric criteria to select the best one. More precisely, the first criterion is the variance of the lengths of the graph edges that belong to the same regular polyhedra. The second criterion is the difference between the reference subunits' angles α(i, j) and the observed angles.

Once the polyhedral graph computed from the input structure has been mapped to the reference graph, we use the precomputed correspondences of the reference graph to map them on our graph and obtain the correspondences between subunits for each element of Γ. This way, we ensure that the function g → σ g is, by construction, a bijective homomorphism.

Graph representation of the group generation

A dual representation of a symmetry group given as a set of permutations, will be a set of rotations. All of these can be obtained as a combination of two generator rotations, one r 3 being a rotation about a 3-fold axis v 3 by an angle 2π 3 for cubic groups, (respectively r n being a rotation about a nfold axis v n by an angle 2π n for dihedral groups) and the other r 2 being a rotation about a 2-fold axis v 2 by an angle π. We should emphasize that these two rotations are present in the 3 possible cubic groups. Then, we may represent all the elements of a point group symmetry as vertices in a Cailey graph [START_REF] Cayley | Desiderata and suggestions: No. 2. the theory of groups: graphical representation[END_REF], whose edges correspond to the two types of the generator rotations. A group element here can be seen as a certain rotation of the symmetrical assembly. For example, a tetrahedral symmetry would have 12 elements (or rotations), an octahedral symmetry would have 24 elements, and an icosahedral symmetry would have 60 elements. The same representation

A B C

Figure 3: Three convex hulls of COMs of perfectly symmetrical tetrahedral assemblies.

The black edges correspond to the maximum common topology between all the possible tetrahedral configurations. The green and blue edges are present only in some convex hull topologies. (A) and (C) show two typical topologies, while (B) shows a degenerated example, where some of the convex hull faces are rectangular instead of being triangular. By applying some random noise to the assembly shown in (B), we can obtain a convex hull with randomly chosen green and blue edges.

holds for dihedral symmetry groups. These, however, will have one generator being r n , a rotation about the n-fold axis by an angle 2π n instead of r 3 as the first generator. The number of elements is dependent on their order, 2n elements for D n . It is easy to demonstrate that all the elements in the groups can be obtained from any initial element by successively applying a combination of the two generator rotations, such that the Cailey graphs are connected, as it is shown in Figure 4. We can also see that a combination of the two generator rotations only produces the group elements, such that the graphs are finite.

Geometry of multiple axes of symmetry

A rotation operator can be uniquely represented by an axis and an angle of rotation [START_REF] Popov | Rapid determination of rmsds corresponding to macromolecular rigid body motions[END_REF]. In our case, to determine the rotation operators, it is convenient to only work with their axes, since the angles are already constrained by the symmetry group. These axes are always placed in an identical configuration with respect to each other, and the position of two axes is sufficient to determine the positions of all other axes [START_REF] Ritchie | Spherical polar fourier assembly of protein complexes with arbitrary point group symmetry[END_REF]. Thus we construct our basis with the axes of the two generator rotations defined above. The angle a between the two basis axes is uniquely defined by the type of symmetry group, as it is shown in Figure 5. It is more convenient, however, to use cosine of the angle α = cos(a), such that

α Tetrahedral = 1 √ 3 α Octahedral = 2 3 α Icosahedral = φ √ 3 , ( 7
)
where φ is the golden number ( √ 5 + 1)/2. The angles a will be then

a Tetrahedral ≈ 54.7 • a Octahedral ≈ 35.3 • a Icosahedral ≈ 20.9 • . (8)
It is interesting to note that the dihedral symmetry group D n can also be rigorously described with two generator axes. One is a n-fold axis v n defining the C n symmetry and the second is a 2 fold axis v 2 perpendicular to it, such that α Dihedral = 0 and a Dihedral = 90 • . For any 3D rotation r g , its rotation axis can be expressed in a basis spanned by 3 axes v 2 , v 3 , and v 3 × v 2 (respectively v 2 , v n , and v n × v 2 for a dihedral group D n ), such that the associated rotation quaternion Qg is written as

Qg ≡ [s g , q g ] = [s g , a g v 3 + b g v 2 + c g v 3 × v 2 ] (9)
for cubic groups and

Qg ≡ [s g , q g ] = [s g , a g v n + b g v 2 + c g v n × v 2 ] (10)
for dihedral groups. The rotation quaternion Q2 Qg obtained by applying r 2 after r g has the following coefficients,

s = -(αa g + b g ) a = -(c g ) b = s g -(αc g ) c = -(a g ). (11)
Similarly, the rotation quaternion Q3 Qg obtained by applying r 3 after r g has the following coefficients,

s = - s g 2 + √ 3 2 (a g + αb g ) a = - a g 2 + √ 3 2 (s g + αc g ) b = - b g 2 + √ 3 2 c g c = - c g 2 + √ 3 2 b g . (12)
Finally, for a dihedral group D n , the rotation quaternion Qn Qg has the following coefficients,

s = cos 2π n s g + sin 2π n a g a = cos 2π n a g + sin 2π n s g b = cos 2π n b g + sin 2π n c g c = cos 2π n c g + sin 2π n b g . (13)

Optimization of the rotations

Here we will use the same notations as in our previous work on cyclic symmetries [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF] . From now on, for simplicity, we will only write equations for the cubic group. Indeed, the equations for the dihedral group are obtained by substituting the index n for the index 3. Our goal is to minimize the loss function defined in equation 4. For each element g of the chosen symmetry group, the contribution to the loss function is the RMSD between r g (A) and A g = {a σg(i),j }. According to the RMSD master equation (3) from [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF] with B = A g , we can say that A and B have the same COM, so the translational part of RMSD becomes null and we obtain

RMSD 2 (r g (A), A g ) = 4 N q T I g q + 4sq T x g⊥ + x gs , (14)
where

I g =  
(y i,j y σg(i),j + z i,j z σg(i),j ) -(x σg(i),j y i,j + x i,j y σg(i),j )/2 -(x σg(i),j z i,j + x i,j z σg(i),j )/2 -(x i,j y σg(i),j + x σg(i),j y i,j )/2 (x i,j x σg(i),j + z i,j z σg(i),j ) -(y σg(i),j z i,j + y i,j z σg(i),j )/2 -(x i,j z σg(i),j + x σg(i),j z i,j )/2 -(y i,j z σg(i),j + y σg(i),j z i,j )/2 (x i,j x σg(i),j + y i,j y σg(i),j )

  , (15)
and

x g⊥ = i,j a σg(i),j × a i,j /N x gs = i,j (a i,j -a σg(i),j ) 2 /N. (16)
Our aim will be to minimize the sum of squared RMSDs over all elements g of the group Γ. Let us first assume that we know the value of one of the two axes v 3 or v 2 , for example, v 3 . In practice, we first compute v 3 axis as a cyclic axis using the method from [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF], then we alternate the computations of v 2 and v 3 considering the other axis as known. This method converges to machine precision in about 10 iterations. Thanks to the RMSD master equation, we can write the loss function as a function of the axis v 2 as follows,

g∈Γ RMSD 2 g (v 2 ) = v T 2 g∈Γ b 2 g 4 N I g + 2 g∈Γ b g c g 4 N I g [v 3 ] × + g∈Γ c 2 g [v 3 ] × T 4 N I g [v 3 ] × v 2 + 2 g∈Γ a g b g v T 3 4 N I g + 2 g∈Γ a g c g v T 3 4 N I g + 4 g∈Γ s g b g x T g⊥ v 2 + g∈Γ a 2 g v T 3 4 N I g v 3 + 4 g∈Γ s g b g x T g⊥ v 3 + g∈Γ x gs . (17)
We can rewrite this equation as the following minimization problem with respect to v 2 , arg min

v 2 v T 2 Av 2 + b T v 2 + c s.t. v T 2 v 2 = 1 v T 3 v 2 = α. (18)
The two constraints come from the unit norm of the rotation axes and the geometry of the generator axes. The above equations has the following coefficients,

A = g∈Γ b 2 g 4 N I g + 2 g∈Γ b g c g 4 N I g [v 3 ] × + g∈Γ c 2 g [v 3 ] × T 4 N I g [v 3 ] × b T = 2 g∈Γ a g b g v T 3 4 N I g + 2 g∈Γ a g c g v T 3 4 N I g + 4 g∈Γ s g b g x T g⊥ c = g∈Γ a 2 g v T 3 4 N I g v 3 + 4 g∈Γ s g b g x T g⊥ v 3 + g∈Γ x gs . (19)
Similar equations can be written for the optimization of the loss function with respect to v 3 .

2D trust-region optimization problem

The optimization problem ( 18) can be efficiently solved by reducing it to the standard form of the trust-region subproblem. However, in our particular case, we can use one of the constraints in eq. ( 18) to project the optimization problem to a two-dimensional subspace. This allows us to solve it analytically, as we explain below.

First of all, it is convenient to chose an orthonormal basis (v x , v y , v 3 ) and rewrite the vector v 2 in this basis as

v 2 = αv 3 + xv x + yv y . ( 20
)
Then, the optimization problem ( 18) reduces to arg min

x,y

x 2 v T x Av x + 2xy v T x Av y + y 2 v T y Av y +x 2αv T x Av 3 + b T v x + y 2αv T y Av 3 + b T v y +α 2 v T 3 Av 3 + b T v 3 + c s.t. x 2 + y 2 = 1 -α 2 . ( 21
)
To solve it, we find stationary points of the corresponding Lagrangian L(x, y, λ), L(x, y, λ) = kx 2 + 2lxy + my 2 + 2px + 2qy + λ(x 2 + y 2 -1 + α 2 ), [START_REF] Dryzun | Quantitative symmetry and chiralitya fast computational algorithm for large structures: Proteins, macromolecules, nanotubes, and unit cells[END_REF] with the following coefficients

k = v T x Av x l = v T x Av y m = v T y Av y p = αv T x Av 3 + 1 2 b T v x q = αv T y Av 3 + 1 2 b T v y . (23)
Assigning the partial derivatives of the Lagrangian to zeros, we arrive to the following system of equations,

     kx + ly + p + λx = 0 lx + my + q + λy = 0 x 2 + y 2 = 1 -α 2 . ( 24
)
After eliminating λ we obtain

lx 2 + (m -k)xy -ly 2 + qx -py = 0 x 2 + y 2 = 1 -α 2 . ( 25
)
Finally, we exclude the last equation by changing the variables and introducing the new optimization variable t,

x = 2t √ 1 -α 2 1 + t 2 ; y = (1 -t 2 ) √ 1 -α 2 1 + t 2 . ( 26
)
Then, making the change of variables and multiplying the first equation by non-zero (1 + t 2 ) 2 we obtain,

-l(1 -α 2 ) + pt √ 1 -α 2 t 4 + 2 (1 -α 2 )(k -m) + t √ 1 -α 2 q t 3 + 6(1 -α 2 )lt 2 + 2 1 -α 2 )(-k + m) + √ 1 -α 2 q t -(1 -α 2 )l -p = 0.
(27) This is our final fourth-order algebraic equation, whose roots can be found analytically [START_REF] Neumark | Solution of cubic and quartic equations[END_REF]. After finding all of its roots, we discard the complex ones, then compute the corresponding values of x and y, substitute them in the original quadratic function [START_REF] Neumark | Solution of cubic and quartic equations[END_REF] and choose the pair of x and y that gives the smallest value. We also additionally test the case of y = -√ 1 -α 2 and x = 0 that has been excluded during the change of variables in eq. (26).

Results and Discussion

Examples

Figure 6 presents an example of symmetry axes detection for each of the cubic groups, i.e. tetrahedral, octahedral and icosahedral, and for a dihedral group of order 6. These assemblies do not possess any particular computational difficulty. Indeed, their asymmetric subunits are composed of a single chain, which makes the first step of optimization problem 5 trivial. Some assemblies contain more chains than the number of asymmetric subunits expected from their point group symmetry. Each subunit thus must be composed of several chains. For example, Figure 7 shows the 5t0v structure, which is an octahedral assembly with 48 chains and a stoichiometry of A24B24. This example demonstrates that our method determines symmetry axes in assemblies where the asymmetric subunits are composed of multiple chains. We should also note that in this case it is important to rigorously take into account all the chains, since the angular difference in the axis determination can be as large as 1 • if only chains A or B are considered. While scanning the PDB, we found several assemblies that are classified with a low-order symmetry group, but can alternatively possess a higher symmetry group. For example, Figure 8 shows the 1ocw structure, which is a perfect C 4 assembly with a stoichiometry of A4B4 and the RMSD loss of 0 Å. Our algorithm also detects a D 4 pseudo-symmetry with the RMSD loss of 2.68 Å, which is rather low. The visual inspection of this protein confirms this possibility (see Fig. 8). Similarly, we also discovered some assemblies with cubic symmetries that were labelled as cyclic in the PDB database. Figure 9 shows two of such examples. One is the 4itv protein labelled as C 2 (RMSD loss of 4.44 Å), but also possessing a tetrahedral symmetry with the RMSD loss of 10.94 Å. The other is the 5hpn protein labelled as C 5 (RMSD loss of 0.68 Å), but also possessing an icosahedral symmetry with the RMSD loss of 0.56 Å.

Comparison with other methods

We compared our approach with two other published methods following the comparison strategy from our previous work on symmetry detection in cyclic protein assemblies [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF]. More precisely, we compared it to the results published by David Avnir and colleagues [START_REF] Dryzun | Quantitative symmetry and chiralitya fast computational algorithm for large structures: Proteins, macromolecules, nanotubes, and unit cells[END_REF][START_REF] Pinsky | Continuous symmetry analyses: Cnv and dn measures of molecules, complexes, and proteins[END_REF]. We will refer to it as to A B C CSM (Continuous Symmetry Measure). We also compared our method to the the one from Emmanuel Levy [START_REF] Levy | 3d complex: a structural classification of protein complexes[END_REF], and will refer to it as to Levy. Please refer to the first part of our paper [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF] for more details.

For the comparison, we have selected all dihedral assemblies presented in the original CSM publications [START_REF] Dryzun | Quantitative symmetry and chiralitya fast computational algorithm for large structures: Proteins, macromolecules, nanotubes, and unit cells[END_REF][START_REF] Pinsky | Continuous symmetry analyses: Cnv and dn measures of molecules, complexes, and proteins[END_REF]. These are listed in Table 1. We have also complemented these assemblies with three examples of cubic groups, 5x47 with tetrahedral symmetry, 4p18 with octahedral symmetry, and 4zor with icosahedral symmetry. Table 1 lists the execution time and the symmetry measure (RMSD value) for the three tested methods. As in the cyclic case [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF], it clearly shows that our method scales with the size of the input assembly much better than the two other methods. This is especially notecable for large assemblies. Regarding the accuracy of the obtained results, it is typically much better than in the Levy method for high-order symmetries. As we have mentioned in the first part of this work, comparison to CSM is trickier because this method considers more atoms than we do. Therefore, the additional atoms add more freedom to the CSM method when it chooses the correspondences between these, which can explain small differences in the computed RMSD values. For example, in the 1f52 case CSM reports a smaller RMSD measure than we do (0.15 Å vs. 0.19 Å).

90°B

A Figure 8: The 1ocw protein colored in blue for the A chains and red for the B chains. A -as seen from the 4-fold axis. B -as seen from a 2-fold axis computed with our method.

Exhaustive analysis of symmetrical structures in the PDB

To demonstrate the efficiency of our approach, we exhaustively analyzed all the structures labelled as symmetrical in the PDB. To do so, we downloaded their biological assemblies (about 40,800 cyclic, 9,800 dihedral and 1,300 cubic examples as for January 2018) and assessed the symmetry for each of these. Figure 10 plots the distribution of the RMSD symmetry measures for assemblies with different types of symmetry. We should note that there are many structures with a very low RMSD value (< 0.001 Å), which is the precision of the pdb format. These are typically obtained by replicating subunits with crystallographic symmetry or BIOMT transforms, so they have a perfect symmetry. Regarding all other structures, we can see that all the three distributions of cyclic, dihedral, and cubic groups follow the same law in the log-log scale. The maxima of the distributions belong to the range of 0.2-0.5 Å, and there are no noticeable differences between the shapes of all of these.

Another interesting question we are able to answer using our tool is whether the degree of asymmetry is related to the size of the assembly under consideration. In other words, we can study if the RMSD symmetry measure is related with the radius of gyration of the symmetrical assemblies. A geometrical intuition would suggest that as the angular uncertainty should stay constant with the size of the assembly, and of protein assembly grows larger, the imperfections of its symmetry become more pronounceable. Visually, we would expect a linear correlation between the RMSD symmetry A B a AnAnaS and Levy times were measured on a Windows laptop equipped with an Intel i7 @ 3.1 GHz.

b CSM times and CSM symmetry measures were taken from [START_REF] Dryzun | Quantitative symmetry and chiralitya fast computational algorithm for large structures: Proteins, macromolecules, nanotubes, and unit cells[END_REF] and [START_REF] Pinsky | Continuous symmetry analyses: Cnv and dn measures of molecules, complexes, and proteins[END_REF] with a different, 7 year older, CPU. However, we believe that the order of magnitude of these timings is still correct. c For these structures, the biological assembly was used. measure and the radius of gyration of the assemblies. However, it is not the case in reality. Indeed, Figure 11 does not demonstrate any clear relation between the size and the imperfection of the PDB assemblies, and the correlation between these two variables is only about 0.1. Interestingly enough, large assemblies are very well organized with sufficiently small values of the RMSD measure. This is one of the reasons behind our choice of RMSD as the symmetry measure instead of its normalization by the size of the structure (as it is often done in other methods [START_REF] Dryzun | Quantitative symmetry and chiralitya fast computational algorithm for large structures: Proteins, macromolecules, nanotubes, and unit cells[END_REF][START_REF] Pinsky | Continuous symmetry analyses: Cnv and dn measures of molecules, complexes, and proteins[END_REF]). We should specifically add that in the case of very small assemblies, we consider them symmetric only if the corresponding RMSD measure is smaller than half of the radius of gyration of the assembly.

How good are symmetry annotations in the PDB?

Our tool also allows to assess the overall quality of annotations of symmetrical assemblies in the PDB. More precisely, we compared the highest symmetry group suggested by our method with the group provided in the PDB. If these two groups are different, there are two types of possible errors. First, one of the two groups can be a subgroup of the other one (e.g. C 4 is a subgroup of D 4 ). This type of errors simply results from a difference of sensibility between the annotation methods. We call the groups compatible. Second, the two groups may also be incompatible (e.g. C 4 and D 5 ). This case means that one of the two results is wrong and a careful visual inspection is generally required.

Table 2 lists the results for 51,358 PDB structures. In 50,378 cases (98.1% of all the cases), the symmetry group annotated by the PDB is the one found by our method. These cases are located at the green diagonal of the table. Red cells show the incompatible groups, while white cells show the compatible groups. Our method is generally more sensitive compared to the PDB annotation. Indeed, there are 845 structures (1.6%) for which it finds a higher order compatible group, while only in 125 cases (0.2%) the PDB annotated compatible group has a higher order. Finally, there are only 13 cases (0.03%) that present incompatible groups. We have visually inspected all of these structures. The two of these annotated as T and detected as C 5 are 4aod and 4aoe, for which the biological assemblies are indeed C 5 . The 11 other cases have uncertainties between C 2 and C 3 annotation. In all of these cases, both symmetries are detected by our method, and the difference of RMSD between the two symmetries is smaller than 1 Å. Moreover, some of these examples have less than 5 amino acids in each chain, and are at the limit of the usability of the annotation techniques. We can also mention two particular cases. One is 3alz, for which both perfect C 3 and C 2 axes are detected, and is actually a part of a D 3 assembly. The other is 3aqq, which is annotated as C 2 in the PDB, but looks much more like a partial C 3 assembly.

The first column of Table 2 lists 75 structures for which AnAnaS was not able to detect symmetry. There are 4 reasons that explain this:

• For the 6 icosahedral structures, we ran out of memory at the discrete optimization step. Thus, no results were outputted and we considered these cases as assymetric.

• Some structures have missing or additional chains that are not supported by our program. For example, 2zl2 has a D 7 symmetry but contains 24 chains, 10 of them being very small peptides. AnAnaS expects a multiple of 14 chains as input to test a D 7 symmetry and, therefore, does not test it. However, if we remove these small peptides, we detect a D 7 symmetry with an RMSD of 0.35 Å.

• Some structures are at the edge of the threshold that we set up for the assemblies to be symmetric. More precisely, as we explain it below, RMSD must be smaller than 7 Å and also smaller than half of the radius of gyration.

• Finally, some structures do not possess the symmetry annotated in the PDB. For example, 2ol9 is the structure of two identical peptides translated with respect to each other, and these are annotated as C 2 , while a C 2 symmetry would necessarily require a rotation between the two peptides.

Computational details

We implemented the AnAnaS method using the C++ programming language. It is available as a standalone executable and also as a module with graphical user interface for the SAMSON software platform. We can also provide the source code upon request.

We have exhaustively assessed our program with all the structures labeled as symmetric in the PDB. This demonstrates the reliability and robustness of our method overall, and its heuristic for the discrete optimization steps in particular. Running the tests on all of these structures took us about 10 hours on a Windows laptop equipped with an Intel Core i7 @ 3.1 GHz CPU. For all the examples we tested, the running time was largely dominated by the multiple sequence alignment, which is required to compare the relevant alpha carbons in different subunits. Only in one case (2qzv) with a D 48 symmetry, the computational bottleneck turned out to be the graph matching step.

We should also say that if no symmetry group is specified by a user, then the program exhaustively tests all the symmetry groups that are consistent with the number of chains in the input assembly. Also, we label an assembly as symmetric only if the corresponding RMSD measure is smaller than 7 Å and smaller than half of its radius of gyration. The second condition is added to filter out very small asymmetric assemblies.

Conclusions

This work extends our previous cyclic symmetry detection method [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF] for high-order point groups. It required to develop a robust heuristic algorithm that perceives the correspondence between asymmetric subunits, and also to extend the constrained quadratic optimization problem from [START_REF] Pagès | Analytical symmetry detection in cyclic protein assemblies[END_REF] to multiple symmetry axes with mutual constraints. Using the quaternion arithmetic, we expressed the constrained optimization problem as a 2D trust-region subproblem and found its solution analytically. We have compared out method with two other published techniques that can detect symmetry in high-order symmetrical assemblies and demonstrated that it is generally much more robust and efficient.

We have demonstrated the efficiency of our method on all the structures marked as symmetric in the PDB, including those with multiple chains per asymmetric subunit or with pseudo-symmetry. It allowed us to verify symmetry annotations in the PDB and detect several inconsistencies in the annotations. For example, in 1.6 % of the cases, we detected a higher symmetry group compared to those provided in the PDB. We have also compared structural organization of protein assemblies with different point group symmetries and concluded that these follow the same distribution laws. Finally, we have detected that the angular impurity in symmetry does not scale with the size of the assemblies. More precisely, very often these are the largest and high-order symmetry systems that are organized the most regularly.

The method is available at https://team.inria.fr/nano-d/software/ananas/. The SAMSON GUI-assisted module is available at http://samson-connect.net/.

Aknowledgements

Figure 1 :

 1 Figure 1: Point group C 5 is illustrated. (a) Set of points with a cyclic symmetry of order 5; (b) All the rotations that keep this set of points fixed; (c) Associations between the elements of C 5 , permutations of the points and rotations.

Figure 2 :

 2 Figure2: Examples of finding subunits in a C 5 complex containing 10 chains, which leads to 2 sets and 5 subunits. Each chain is shown as a shape, with a local coordinate frame computed with a structural superposition. A correct clustering in this example should determine a set with the chains 1,2,3,4,5 (light grey) and another one with the chains 6,7,8,9,10 (dark grey). The correspondence between the subunits after a rotation is represented by the colored arrows. Chains belonging to the same subunits are drown with the same outline colors. After we define the two sets, the first subunit is created by taking arbitrary chains from each cluster (1 and 7 in (A), 1 and 6 in (B), 1 and 9 in (C)). The other subunits are created using the correspondences between the chains in the sets. Note that the assembling of subunits is not unique. In (A) and (B) the subunits seem to be assembled more naturally, as they contain chains that are spatially proximate to each other. However (C) is a perfectly valid assembly too. All the three subunit assemblies give exactly the same result in terms of symmetry axes and the RMSD loss function.

Figure 4 :Figure 5 :

 45 Figure 4: Cailey graph representation of the different cubic groups, (A) the tetrahedral, (B) the octahedral, (C) the icosahedral, and (D) the dihedral D 8. Each vertex represents a group element, which is also a rotation of the symmetrical assembly. The directed edges are the generator rotations applied to the group elements. The blue edges are the r 3 generator rotations, the red edges are the r 2 generator rotations, and the green edges are the r n generator rotations. Each vertex has indegree and outdegree of two, such that the graph is balanced.

Figure 6 :

 6 Figure6: Four examples of symmetrical assemblies with their axes. All of these are seen from a 3-fold axis except for the last one, seen from a 6 fold axis. The order n of each axis is represented with a regular n-gone, except of order 2 represented with a rhombus. A -A tetrahedral assembly (1d0i) with the RMSD loss of 0.36 Å. B -An octahedral assembly (1bfr) with the RMSD loss of 0.22 Å. C -A perfect icosahedral assembly (1stm) with the RMSD loss of 0.0 Å. D -A dihedral D 6 assembly (1f52) with the RMSD loss of 0.20 Å. This illustration and all the illustrations below were produced in SAMSON (www.samson-connect.net).

Figure 7 :

 7 Figure 7: The 5t0v octahedral assembly. The homologous chains are colored with the same color. A -The chains of the first type form an octahedral assembly with the RMSD loss of 2.94 Å. B -The chains of the second type also form an octahedral assembly with the RMSD loss of 2.67 Å. The axes are slightly different from the first assembly, with about 1 • of difference. C -y The axes are computed for the full assembly, with the RMSD loss of 2.83 Å.

Figure 9 :

 9 Figure 9: A -The 4itv protein classified in PDB as C 2 (RMSD loss of 4.44 Å), also has a tetrahedral symmetry with the RMSD loss = 10.94 Å. B -The 5hpn protein classified in PDB as C 5 (RMSD loss of 0.68 Å), also has an icosahedral symmetry with the RMSD loss = 0.56 Å.

Figure 10 :

 10 Figure 10: Distribution of the RMSD symmetry measure for different types of symmetry shown in a log-log scale.

Figure 11 :

 11 Figure 11: Distribution of the RMSD symmetry measure with respect to the radius of gyration for non-perfectly symmetrical assemblies from PDB.

  chain do Create a new set; Add this chain to the set; end while Number of sets is bigger than expected do Find the two closest sets whose sum of sizes is smaller or equal to the expected set size; if Such a pair of sets is found then Merge the two sets;

	else
	Fail;
	end
	end
	Algorithm 1: Clustering algorithm.

Table 1 :

 1 Comparative results between AnAnaS, CSM and Levy methods for dihedral and cubic molecular assemblies.

Table 2 :

 2 PDB AnAnaS C 1 C 2 C 3 C 4 C 5 C 6 C 7 C 8 D 2 D 3 D 4 D 5 D 6 D 7 D 8 TSummary of the symmetry groups annotated in the PDB (rows) against the ones discovered by AnAnaS (columns). Red cells mark incompatible groups, while white cells mark compatible groups and green cells mark identical groups. For example, first cell shows that there are 54 structures annotated as C 2 in the PDB for which we found a C 1 symmetry (i.e. no symmetry).

	O	I Total

The authors thank Nikolay Mayorov for his 2D trust-region algorithm developed during Google Summer of Code 2015, and Sergei Khashin from Ivanovo State University for his forth order polynomial solver available at http://math.ivanovo.ac.ru/dalgebra/Khashin/poly/index.html. The authors also thank Elvira Kinzina and Andrei Kazennov from MIPT Moscow for their support at the initial stage of the project.

Funding

This work has been supported by L'Agence Nationale de la Recherche (grant number ANR-15-CE11-0029-03).


Philippe Martin Arteris
email: philippe.martin@arteris.com

Design of a Virtual Component Neutral Network-on-Chip Transaction Layer

Research studies have demonstrated the feasibility and advantages of Network-on-Chip (NoC) over traditional bus-based architectures but have not focused on compatibility communication standards. This paper describes a number of issues faced when designing a VC-neutral NoC, i.e. compatible with standards such as AHB 2.0, AXI, VCI, OCP, and various other proprietary protocols, and how a layered approach to communication helps solve these issues.

NoC Layered Architecture

A true Network-On-Chip uses a layered communication approach, similar to the OSI model which is at the foundation of LANs and WAN, but simpler. Arteris NoC defines transaction, transport, and physical layers.

The NoC transaction layer defines communication primitives available to IP blocks that are plugged into the NoC. A Network Interface Unit (NIU) is responsible for converting the foreign IP protocol to the NoC transaction layer.

The transport layer defines information format and transport rules between NIUs. Arteris NoC packet-based transport layer defines how packets are routed, quality of service used for prioritizing packets, etc. The transport layer is completely transaction unaware, and conversely, transaction level is transport unaware (for example, wormhole or store-and-forward packet handling makes no difference at the transaction level).

The physical layer defines how packets are physically transmitted -much like the Ethernet defines the MII, 10Mb/s, 1Gb/s physical interfaces. Again, the physical layer is independent from transaction and transport layers.

Layer independence enables a separate optimization of NoC features: the transaction layer focuses on compatibility with existing IP blocks and busses, transport layer focuses on quality of service and scalability, physical layers on implementation related features such as achieving raw bandwidth, matching clocks, off-chip communication, etc.

VC compatibility handling

To create System-On-Chip (SoC), VCs with mixed sockets ideally plug into an interconnect according to Interconnects use a multi-master and multi-slave protocol and include bridges to the various VC socket protocols.

In practice, the intertwining of transaction, transport and physical levels within standard interconnects, and the interconnect flexibility necessary to handle many application designs makes this ideal approach very difficult. Designs mixing different VC standards look more like Fig 2 , where the interconnect has its own reference socket standard. Bridges to the reference standard are used plug the IP blocks (an alternative is to redesign the IP to support directly the interconnect reference standard, but this is time-consuming and prone to errors, and not possible for off-the-shelf IP). Bridges introduce area and latency penalties, but worse, they also do not support the full set of VC transactions because they are limited by the interconnect protocol and physical design. Compared to this approach, the layered architecture provides a strong advantage to NoCs with respect to compatibility handling:

-The switching fabric managing transport being transaction-unaware, transactions can be customized to the actual set of VCs that plug into the NoC, without altering the transport and physical layers.

-Adding a socket-specific feature to a NIU is easy because the NoC protocol layering drives the process with a simple set of questions:

1. Does the feature require some specific transaction state to be stored in the NIU? If yes, add the state to the standard NIU state lookup tables (which track for example that a Load request is waiting for a response).

2. Does it require information to be exchanged between NIUs? If yes, add it to the packet format.

Since neither adding bits to the packets nor state in the NIUs impacts transaction or physical layers, supporting VC-specific features in the NoC only impacts the corresponding NIU. However, the various VC sockets have many incompatibilities in their basic features that need to be carefully addressed upfront in the NoC transaction layer.

Incompatibility examples and solutions

Many VC compatibility issues must be carefully analyzed and resolved to reach a truly VC neutral NoC. Among these are: ordering model, exclusive accesses and atomicity, bursts, endianness. We will discuss ordering model and exclusive accesses below.

AHB, PVCI, BVCI VCs are fully-ordered between requests and responses. OCP compatible VCs are fully ordered within a thread, but may be multi-threaded with no ordering constraint between threads. AXI and AVCI provide transaction IDs with requests and response, also allowing out-of-order responses. Some of these protocols (such as OCP) support WRITEs without responses, and others (such as AXI) have independent READ and WRITE channels, further obscuring ordering constraints.

The various flavors of ordering models and pipelining schemes between these protocols create a challenge to define a transaction layer compatible with all configurations, while keeping a low NIU gate count.

Arteris NoC protocol uses a packet destination field (called SlvAddr), a packet source (called MstAddr), and a Tag. The ordering model adapts to the fully-ordered AHB, the multi-threaded OCP and the ID-based AXI ordering models using a careful assignment policy of these fields from the OCP or AXI ones such as ThreadID and TID. Further, this policy is flexible and allows NIUs to support one or many simultaneously outstanding transactions and/or targets, scaling their gate count to their expected performance within the system.

As outlined above, the NoC switch fabric itself is unaware of actual NIU field assignment policies -it sees uniform packets with SlvAddr, MstAddr, Tag information, allowing the seamless plugging of OCP, AXI, AHB, etc. IP blocks to the same NoC.

OCP and AXI have introduced new transactions called « lazy synchronization » and « exclusive access » respectively.

They implement non-blocking synchronization between several masters, contrary to the older READEX and LOCK transactions. To implement all these in the NoC, we found out that READEX/LOCK transactions impact transport level (switches take specific decisions when they see LOCK-related packets) but that handling AXI and OCP exclusive access only requires adding a single user-defined bit in the packets, and state information in the NIU. This optional packet bit becomes simply part of a family of similar "NoC services" that can be activated in a particular NoC configuration.

Summary

Arteris NoC transaction layer is compatible with AHB, VCI flavors, AXI, OCP transaction levels and thus allows seamless mixing of IP blocks using these sockets. The layered structure of NoCs allows IP compatibility to be independent from transport and physical layer specification, allowing easy plug-in of new socket protocols to the NoC through appropriate NIUs, while at the same time the transport and physical layers can be adapted to specific QoS requirements or physical environments.

  Fig 1:

Fig 1 :

 1 Fig 1: Ideal system with mixed VC protocols.

Fig 2 :

 2 Fig 2: Usual system with mixed VC protocols

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE


Arshad Jhumka
email: arshadjhumka@hotmail.com

Stephan Klaus
email: klaus@iss.tu-darmstadt.de

Sorin A Huss
email: huss@iss.tu-darmstadt.de

A Dependability-Driven System-Level Design Approach for Embedded Systems

The objective of this paper is to introduce dependability as an optimization criterion in the system-level design process of embedded systems. Given the pervasiveness of embedded systems, especially in the area of highly dependable and safety-critical systems, it is imperative to directly consider dependability in the system level design process. This naturally leads to a multi-objective optimization problem, as cost and time have to be considered too. This paper proposes a genetic algorithm to solve this multi-objective optimization problem and to determine a set of Pareto optimal design alternatives in a single optimization run. Based on these alternatives, the designer can choose his best solution, finding the desired tradeoff between cost, schedulability, and dependability.

Introduction

Computer systems are becoming increasingly pervasive such that our reliance on their continual provision of correct services, in spite of external perturbations such as node crashes, has correspondingly increased. This means that we, as users, want these systems to be dependable [START_REF] Laprie | Dependability: basic concepts and terminology[END_REF]. These computer systems are being deployed in a wide range of areas, from safety-critical systems such as nuclear power plants, fly-by-wire systems to consumer-oriented products such as cars, mobile phones, washing machines etc.

Focusing on consumer-oriented products, cost is a critical issue. Where for safety-critical systems, dedicated processors can be used to execute safety-critical applications, in consumer-oriented products such as cars, safety-critical applications, such as the braking system, may well execute on the same processor as a non-safety-critical application, e.g., temperature control, so as to keep the cost factor tractable. Such a design approach is already taking place in so-called Integrated Modular Avionics (IMA), where safety-critical applications are being executed on the same platform as non-safety-critical applications. Designing dependable systems using such a methodology entails ensuring that non-safety-critical applications do not interfere with the safety-critical ones. Rushby [START_REF] Rushby | Partitioning in avionics architectures: Requirements, mechanisms, and assurance[END_REF] pointed out that, for this approach to work, system partitioning along two dimensions is needed: (i) spatial, and (ii) temporal. In this work, we focus on temporal partitioning, and we assume spatial partitioning.

The concern of temporal partitioning is to ensure that activities in one partition do not interfere with the timing of events in other partitions. In this paper, since tasks of different 'criticalities' may execute on a given processor, we assume that each partition has to contain one task. As suggested by Rushby in [START_REF] Rushby | Partitioning in avionics architectures: Requirements, mechanisms, and assurance[END_REF], static scheduling can be used at the level of partitions, while dynamic scheduling can be used within partitions. Since a partition consists of a single task, in this paper, static scheduling is used for temporal partitioning.

Designing for dependability in the latter stages of the design phase may lead to a complete redesign of the system in the worst case, or can be very costly (as we will show in Sec. 5). This means that dependability aspects must be considered during the earlier part of the design process. Hence, there are certain dependability-driven factors that need to be taken into consideration earlier on. One of these factors is replication. Replication of the more critical tasks is needed to ensure that, in case of a failure of a node upon which a critical task resides, there is at least one other replica of that task still executing to take over.

Replication of the more critical tasks makes the system more dependable, however this may come at the expense of cost (e.g., more resources may be needed to accomodate the extra tasks), and schedulability. This leads to a multi-objective optimization problem (between cost, time and dependability). To address this, we propose a genetic algorithm to solve the multi-objective optimization problem, whereas a set of Pareto optimal design alternatives is generated in a single optimization run. The system designer can then choose the solution that best suits the system's re-1 quirements.

The remainder of the paper is organized as follows. Sec. 2 discusses the related work. In Sec. 3, our system model is presented and the approach to increase the dependability of embedded systems is introduced. Sec. 4 discusses the proposed genetic design space exploration. Sec. 5 provides some experimental results and some insights. We conclude in Sec. 6.

Related Work

For reasons of space, we will focus in the sequel on work related to the provision of dependabilty in embedded systems.

It has been shown that making a given system dependable requires that the given system has sufficient redundancy [START_REF] Gaertner | Redundancy in space in faulttolerant systems[END_REF]. There are two dimensions to be considered, namely: (i) redundancy in space, and (ii) redundancy in time. For example, to support a single node failure, there should be at least two nodes to support continuous service provision (space), or in order to mask a value failure, a reexecution might be needed (time). In this work, we consider the dependability aspect during the design phase and incorporate the necessary redundancy (spatial and temporal) so as to make the system dependable.

[16] developed a fault containment-based approach, whereby applications are first decomposed into their constituents tasks and procedures, which are subsequently transformed to create fault-containment blocks. These blocks are later re-integrated to form dependable applications, which are then allocated to nodes, and scheduled. This approach works when the source code is available, but it is not readily usable at the system design level. A similar approach has been adopted by [START_REF] Dave | Cofta: Hardware-software cosynthesis of heterogeneous distributed embedded systems for low overhead fault tolerance[END_REF]. Up to now, very little work has addressed the problem of integrating dependability as an optimization criterion. [START_REF] Gonzalez | Adaptive fault tolerance and graceful degradation under dynamic hard real-time scheduling[END_REF] developed an approach in which processes/tasks are replicated until the overall dependability utility of the system is maximized. However, this comes at the expense of cost, which is an important constraint in embedded system design. Also, they fail to model factors such as importance.

System-level synthesis approaches based on task graphs can be found in [START_REF] Eles | System level hardware/software partitioning based on simulated annealing and tabu search[END_REF] and [START_REF] López-Vallejo | On the hardware-software partitioning problem: System modeling and partitioning techniques[END_REF]. A comparison of simulating annealing, genetic algorithm and tabu search for the functional HW/SW partitioning using a single optimization criterion is presented in [START_REF] Axelsson | Architecture synthesis and partitioning of realtime systems: A comparison of three heuristic search strategies[END_REF]. [START_REF] Dick | MOGAC: A Multiobjective Genetic Algorithm for Hardware-Software Co-synthesis of Hierarchical Heterogeneous Distributed Embedded Systems[END_REF] and [START_REF] Haubelt | Hierarchical Synthesis of Embedded Systems Using Evolutionary Algorithms[END_REF] address genetic algorithms for multi-objective system-level synthesis. But all these approaches do not consider dependability at all. This is done by [START_REF] Haubelt | ReCoNet: Modeling and Implementation of Fault Tolerant Distributed Reconfigurable Hardware[END_REF], who introduces fault tolerance as a parameter for the design of reconfigurable hardware. Overall, there is a dearth of work that addresses dependability as an optimization criterion in embedded systems design.

A good survey and a detailed description of specification models for embedded systems can be found in [START_REF] Cortés | A survey on hardware/software codesign representation models[END_REF]. Next to state-oriented models and to abstract models like UML, task graphs are a widely spread means for the specification of data-flow dominated embedded systems. Beside the dataflow, most embedded systems contain at least a few controlflow decisions. Therefore, a lot of research was done extending task graphs by control-flow information [START_REF] Eles | Scheduling of conditional process graphs for the synthesis of embedded systems[END_REF][START_REF] Klaus | A Novel Specification Model for IP-based Design[END_REF].

Models

System Model: We assume the following model in this paper. A system S consists of a set of n applications, A 1 . . . A n , and can be represented by a directed graph. Each node in the graph represents an application. An edge from node A j to node A k represents data flowing from application

A j to A k . Each application A i consists of a set of k i tasks, t 1 i . . . t ki i .
It can be represented by a task graph, where each node in the graph represents a task, and the edges between the nodes represent data dependencies between these tasks. By unwrapping each application node with their respective task graph, the system can be represented as a task graph, where nodes denote tasks and edges represent data between tasks. Thus, S consists of m tasks, T 1 . . . T m , respectively.

We denote by R the set of available resources. A cost factor v i is associated with every resource r i ∈ R. Every time an additional resource r j ∈ R is used, the overall system cost increases by v j . Each task t i can be represented as a tuple (R i , E i , C i ), where

• R i ⊆ R is the set of resources task T i can use. • E i ⊆ R i × Z
is the execution time of task T i when using a resource in R i (modeling heterogeneity).

• C i represents the importance of the task. Importance is necessary when dependability is taken into account. The higher the importance of a task, the more critical it is for providing correct services.

There are two possible ways to model replication. First, we can associate a replication factor with every application of the system. The higher the replication factor, the more critical the application is for the system. Suppose a system designer decides to assign a replication factor of 2 to a given application. We can then duplicate every task of that application. Second, since a system can be represented as a task graph, we can directly associate a replication factor with every task. In this work, we adopt the second approach. Since every task in an application is possibly not critical, it is not necessary to replicate all the tasks of that application. Hence, the ability to selectively replicate tasks is crucial.

When replicating a given task, then the initial task graph (hence called task graph) needs to be extended into a repli-cation task graph. The replication task graph is constructed from the task graph as follows (Fig. 1):

• For every edge (n 1 , n 2 ) pointing to a node n 2 in the task graph, there is a an edge (n 1 , n 2 ) pointing to every replicated node n 2 in the replication task graph.

• For every edge (n 1 , n 2 ) starting at node n 1 in the task graph, there is an edge (n 1 , n 2 ) starting at every replicated node n 1 in the replication task graph. When replication is introduced in the system-level design process, then another constraint regarding the binding has to be taken into account. If a task is duplicated in order to increase the dependability of the system, then it makes no sense to bind the replicated task to the same resource as the original task. If a failure occurs, then the replicas should not be affected by this failure. So, each replication task of T i has to be bound to a resource r ∈ R i in such a way that no other replica of T i is bound to the same resource.

Based on failure probability, importance, and replication, the dependability is calculated as follows:

Dependability = T i=1 C i * (1 -p i fi ) (1)
whereas: T: number of tasks; f i : replication factor of task t i ; C i : importance of t i ; p i : failure probability of t i . We assume that the replication factor, and importance of a task to be provided by the system designer. We further assume that the failure probability of a task is equal to the failure probability of the node it is running on, and the probability of failure of the node can be obtained from field data.

Fault Model: Since we are concerned with enhancing dependability via replication, the underlying fault model is resource crashes.

System-Level Synthesis

The purpose of system-level design or design space exploration is to determine a set of suited implementation alternatives. This process can be detailed as the selection of the necessary resources (allocation) and the mapping of the functional units onto the selected architecture in space (binding) and time (scheduling). Naturally, this is a multi-objective optimization problem, since const, time, and especially dependabilty are taken into account. Therefore, no universal definition of an optimum can be given. One solution for an appropriate definition of optimality in multi-objective optimization was given by Vilfredo Pareto in 1896. This definition expresses basically that a solution is Pareto-optimal if there exists no feasible other solution, which would decrease some objective without causing a simultaneous increase in at least one other objective.

Figure 2 illustrates a solution space S for two optimization objectives f 1 and f 2 , respectively. The solution space is restricted by the design constrains (dashed line) and by the Pareto-front (solid line), which consists of all Paretooptimal solutions. Among the Pareto-optimal solutions S1 to S3, S1 is the best one according to f 1 and S3 according to f 2, respectively. S2 dominates S1 with respect to f 2 and S3 regarding to f 1.

Genetic Design Space Exploration

Genetic algorithms are choosen to solve this complex multi-objective optimization problem, since a complete search is in general not feasible due to the size of the solution space. This optimization heuristic [START_REF] Goldberg | Genetic Algorithms in Search Optimization and Machine Learning[END_REF] is inspired by Darwin's theory of evolution and is very well-suited for problems with large, non-convex search spaces and multiobjective optimization goals.

Basically, a set of possible solutions (the population) is modified by crossover, mutation, and selection operations. Thereby, the selection is controlled by a fitness function, which measures the quality of each individual (one possible solution) in the population. The general flow is illustrated in Algorithm 1. First, the population is initialized by assigning randomly some values to each individual. This population is evaluated by a Pareto ranking, and then a subset of the population is mutated according to a given probability. Subsequently, new individuals are generated by crossover operations, which replace randomly individuals in the current population.

Algorithm 1 Genetic Design

Representation

Each individual represents one possible solution for the optimization problem of embedded system design. Therefor, a coding has to be found, which is suited to represent all design alternatives of the envisaged solution space. First of all, an individual must hold information about a complete implementation of an embedded system (i.e., allocation, binding, and scheduling). Allocation is being represented by a set of resources on top of which the system can be implemented. The binding assigns a resource instance to each task considering both the current allocation and the possible binding alternatives of the task. The schedule is represented by unique priorities, which are utilized in the evaluation phase to determine the start and end times for all tasks. In addition, information about replication factors of each task is assigned to each individual.

Figure 3 depicts the representation of the replication factors for the example of Figure 1. It can be easily seen, that the list contains a replication factor for each task.

Genetic Operators

The genetic operators for initialization, mutation, and crossover are divided into parts for allocation, binding, scheduling, and replication factors too. The operators for allocation, binding, and scheduling are working on the extended task graph including the replicas according to the replication factors. The population is initialized by randomly generating some individuals, whereas only feasible solutions are generated. For the initialization of the allocation a subset of all avaliable resources is chosen. The binding assigns randomly a resource instance to each task considering the current allocation. The schedule is initialized by simply giving each task an unique priority and the replication factors recive a value between 1 and a maximal replication count.

The mutation and the crossover operators for the replication factors are illustrated in Figure 4. Mutation chooses randomly a position and changes the replication factor for this task. For the crossover the lists of replication factors are splitted into to sets: In the first set, the factors according to the first parent are taken. The second set is generated accordingly from the data of the second parent. For the second child the roles of the parents are swapped.

Figure 4. Genetic Operators

The allocation is mutated by removing or adding randomly some resources. As a result, the allocated set may be too small, so that tasks exist, for whom no resource is avaliable. In contrast to "punishing" such individuals, the allocation is repaired by adding randomly missing resources. The binding is mutated by assigning a different resource from the allocation to some tasks, if possible. For mutating the schedule, the priorities of two tasks are randomly exchanged.

The crossover of the allocation is done by choosing randomly resources from one parent and by completing this set by resources from the other parent. For the crossover of the binding the sets are divided into two parts: In the first set, it is tried to bind each task according to the binding of the first parent. If this is not possible, then the second parent is taken. In the second part the same is done with swapped roles of the parents. Crossover of schedules is done by dividing the tasks into three parts. Hence filling the middle with the priorities of one parent, and the first and last part with the priorities of the other parent, if they were not used yet. The rest of the data structure is filled with the remaining priorities.

Evaluation

In the evaluation phase the costs, the dependability and both the start and end times of all tasks are calculated. The cost of a single solution is determined by summing up the costs of all used resources. The dependability is calculated according to eq. ( 1) considering importance, failure probability, and replication factors of the tasks. A scheduling based on the determined binding and on the calculated priorities results in exact execution times for each task. So, the overall execution time and possible deadline misses of a single task can be derived. The scheduling is done by globally ordering the tasks on each resource according to the data dependencies and the priorities. This circumvents the disadvantage of list scheduling algorithms. These algorithms unveil the problem of timing anomalies. Such a anomaly occurs in the case that the overall performance of a system increases, if a task is not started immediately after this tasks becomes ready and the resource is available.

The evaluation is based on a Pareto ranking, which sorts all solutions according to the amount of individuals which are not dominated by this solution. The resulting scores of each solutions are depicted as attributes in Fig. 2. E.G., S5 is not dominated by S3, but by all other solutions. At the end of the design space exploration the complete set of Pareto-optimal solutions is passed to the designer, the socalled human decider, who can now select his best design solution.

Results

This section demonstrates the figures of merit of the proposed approach, by means of a small but expressive example. Here, we take |R| = 1. The advantages of directly considering dependability in the system level design process are illustrated by means of the example of Figure 5. The results of a classical design approach are directly compared in the following to our comprehensive system-level synthesis approach.

Figure 6 depicts an implementation variant generated by a classical design space exploration without considering dependability at all. This implementation exploits the resources R1 and R2 resulting in overall costs of 270. Faster, but more expensive solutions exist when exploiting R2 and R3. Late replication v/s our approach: Dependability issues are being addressed by many designers of new systems. Conceptually, dependability is treated as an "add-on", e.g., it is accomplished by doubling the complete implementation. The resulting system has a dependability value of 6 according to eq. ( 1), but the cost is doubled too. This situation occurs when dependability is considered late in the system design phase. On the other hand, Figure 7 shows an implementation that is derived from applying our proposed approach, assuming all tasks to have the same importance value (hence, same replication factor). A dependability value of 6 is achieved too, but by exploiting only three resources (rather than 4). This results in system cost of 470 instead of 540 for the naive, i.e., doubled version. All Pareto optimal solutions generated by the proposed genetic algorithm are presented in Table 1. As it can easily be seen, the solution space spans from of expensive and highly dependable systems to cheap and fast implementations without any dependability. The implementation presented in Figure 6 is not part of this list, since a dependability value of 4.25 can be reached without increasing the costs or the overall execution time of the system (alternative 9). This result shows that it pays off to consider dependability issues early in the design phase, as opposed to a-posteriori add-ons. Different replication factor: A design space exploration was performed by applying a maximal replication factor of 4. The importance parameter C i of tasks T 4 and T 5 is set to 2, as these tasks are considered as the most critical once. and therefore should have a higher replication factor. The replication factor for each task are given in Table 2. The resulting implementation exploits four resources, i.e., two R1, one R2 and R3. The overall cost sums up to 570, the dependability value results in 6.375, and the overall execution time is 1300. Compared to previous cases (for the same resource cost)we only pay an extra 13% time for a gain of more than 16%. This result shows again that early consideration of dependability is important, as well as an individual assignment of replication factors to tasks, rather than to applications (finer grained). This holds especially for safety-critical systems, where resources are abundant. Obviously replicating non-critical tasks does not necessarily lead to a significant increase in dependability.

Alt. Execution Time

Table 2. Replication Factors

Conclusion

This paper presents a novel system-level design method for embedded systems which considers dependability as an optimization criterion. Replication is introduced at task level to increase the dependability of the system. Further, instead of determining a single solution only, a set of Pareto optimal design alternatives is calculated by the proposed genetic algorithm in a single optimization run to help the system designer in finding the required tradeoff between runtime, costs, and dependability. As demonstrated by the de-tailed example, taking dependability directly into account leads to superior implementations. Currently, we allow tasks of different applications to be potentially distributed across different processors. In other cases, this may not be feasible, since it does not provide fault containment. So, a subject of future work will be system design under more complex failure models.

Figure 1 .

 1 Figure 1. Initial and Replication Task Graph

Figure 2 .

 2 Figure 2. Illustration of Pareto-Front

Figure 3 .

 3 Figure 3. Genetic Representation of Replication Factors

Figure 6 .

 6 Figure 5. Example Task Graph

  Figure 7. Dependable Implementation

Table 1 . Implementation Alternatives

 1

			Cost Dependability
	1	1550	270	4.5
	2	1025	670	5.5
	3	1000	390	4.50
	4	1150	590	5.50
	5	1150	390	5
	6	1350	440	5.75
	7	1300	540	6.0
	8	1550	470	6.0
	9	1300	270	4.25
	10	1550	320	5.25

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE