Guillaume Guguin

Olivier Allix

Pierre Gosselet

Stéphane Guinard

On the computation of plate assemblies using realistic 3D joint model: a non-intrusive approach

Keywords: assembly, non-intrusive coupling, bolt

Most large engineering structures are described as assemblies of plates and shells and they are computed as such using adhoc Finite Element packages. In fact their computation in 3D would be much too costly. In this framework, the connections between the parts are often modeled by means of simplified tying models. In order to improve the reliability of such simulations, we propose to apply a non-intrusive technique so as to virtually substitute the simplified connectors by a precise 3D nonlinear model, without modifying the global plate model. Moreover each computation can be conducted on independent optimized software. After a description of the method, examples are used to analyze its performance, and to draw some conclusions on the validity and limitation of both the modeling of junction by rigid connectors and the use of submodeling techniques for the estimation of the carrying capacity of bolted plates.

Introduction

The simulation of large structures undergoing complex local nonlinear phenomena is still a major scientific and industrial challenge. One of the main difficulties originates from the difference of length scale between the global response of the structure and the localized phenomena. To address those problems a first type of computational approach is based on homogenization, as FE 2 [1] but it works well as long as the scales are sufficiently separated. To overcome this limitation concurrent multiscale methods have been developed. They are often based on domain decomposition techniques like FETI [START_REF] Wyart | A substructured FE-shell/XFEM-3D method for crack analysis in thin-walled structures[END_REF], FETI-DP [START_REF] Amini | Multi-scale domain decomposition method for large-scale structural analysis with a zooming technique: Application to plate assembly[END_REF] or the LATIN multiscale method [START_REF] Ladevèze | A micro / macro approach for parallel computing of heterogeneous structures[END_REF][START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF] and its optimization for the multiscale treatment of nonlinear problems [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF][START_REF] Saavedra | On a multiscale strategy and its optimization for combined delamination and buckling simulation[END_REF].

Moreover most of large industrial structures are described as an assembly of plates and shells, whereas local phenomena often require 3D models to be properly analyzed. To deal with such problems, several methods have been applied or developed for the coupling of 2D and 3D models, like the Arlequin method [START_REF] Dhia | Problèmes mécaniques multi-échelles: la méthode arlequin[END_REF][START_REF] Dhia | The arlequin method as a flexible engineering design tool[END_REF], transition elements [START_REF] Gmür | Three-dimensional solid-to-beam transition elements for structural dynamics analysis[END_REF][START_REF] Garusi | A hybrid stress-assumed transition element for solidto-beam and plate-to-beam connections[END_REF], MPCs approaches [START_REF] Mccune | Mixed-dimensional coupling in finite element models[END_REF][START_REF] Donaghy | Dimensional Reduction of Analysis Models[END_REF] or Nitsche's method [START_REF] Nguyen | Nitsche's method for two and three dimensional nurbs patch coupling[END_REF].

Most of these methods are quite demanding in terms of software development and therefore they are seldom used in industrial packages. To overcome these drawbacks, non-intrusive approaches have recently been proposed [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF]. They are nowadays the subject of extensions and developments: thermoelasticity with GFEM/FEM coupling [START_REF] Plews | An improved nonintrusive global-local approach for sharp thermal gradients in a standard fea platform[END_REF], crack propagation in XFEM/FEM coupling [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF], stochastic simulations [START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF] and dynamics [START_REF] Bettinotti | A coupling strategy for adaptive local refinement in space and time with a fixed global model in explicit dynamics[END_REF][START_REF] Bettinotti | A fast weakly intrusive multiscale method in explicit dynamics[END_REF].

In [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF] a non-intrusive coupling between plate and 3D models was proposed in the case of linear behaviors. The present paper concerns the extension of this approach to the simulation of bolted assemblies of plates where bolts are described with full 3D nonlinear models. Such structures are good candidates for the iterative global-local non-intrusive strategy for two main reasons. First, the detailed computation of a tightened bolt with frictional contact on all surfaces is a very hard and time-consuming task to perform using commercial software. Second, the construction of the reference model, which would correspond to the assembly of a plate model and a 3D model for the bolt, would be very complex. The non-intrusive framework provides answers to both these problems. First, it allows the use of dedicated software for the local computations (in our case, COFAST a parallel software based on the LATIN domain decomposition method [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF]). Second, it allows the easy coupling of a general Finite Element software (here Code Aster from EDF) for the plate computation, with COFAST, because the global model is unchanged during the iterative process. These properties were exploited in [START_REF] Daghia | A micro-meso computational strategy for the prediction of the damage and failure of laminates[END_REF] for the simulation of damage in composite laminates at the meso and the micro scales using dedicated pieces of software.

The non-intrusive framework aims at solving the reference problem iteratively, by solving at each iteration both the global problem with prescribed residual traction at the interface and the local nonlinear problems submitted to prescribed displacement. Several techniques have been proposed to improve the convergence rate of the method by means of acceleration techniques [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF][START_REF] Liu | A non-intrusive global/local algorithm with non-matching interface: Derivation and numerical validation[END_REF] or improved interface conditions [START_REF] Gendre | A two-scale approximation of the schur complement and its use for non-intrusive coupling[END_REF]. The method has therefore common points with so called nonlinear domain decomposition methods (or nonlinear relocalization techniques) [START_REF] Klawonn | Nonlinear FETI-DP and BDDC methods[END_REF][START_REF] Negrello | Substructured formulations of nonlinear structure problems -influence of the interface condition[END_REF] which proved their efficiency and gain in robustness in the case of buckling [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: Application to post-buckling analyses[END_REF], post-buckling [START_REF] Hinojosa | Domain decomposition methods with nonlinear localization for the buckling and post-buckling analyses of large structures[END_REF] and damage analysis [START_REF] Bordeu | Balancing domain decomposition with nonlinear relocalization: Parallel implementation for laminates[END_REF][START_REF] Allix | A relocalization technique for the multiscale computation of delamination in composite structures[END_REF]. Other proposals have been made to take into account the fact that the phenomena of interest are localized, aiming at a better representation of the target model [START_REF] Düster | pq-adaptive solid finite elements for threedimensional plates and shells[END_REF][START_REF] Düster | Applying the hp-d version of the FEM to locally enhance dimensionally reduced models[END_REF][START_REF] Plews | An improved nonintrusive global-local approach for sharp thermal gradients in a standard fea platform[END_REF].

The paper is organized as follows. In section 2, a summary of the reference problem corresponding to the coupling of 2D and 3D models according to [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF] is presented. The non-intrusive algorithm is presented in section 3. In section 4, the results of the iterative coupling are analyzed in the case of a bolted joint.

These results are compared to those corresponding to the plate solution and to the submodeling technique. In section 5, the influence of some parameters is assessed regarding the rate of convergence of the iterative process and regarding the accuracy of the coupled model compared to a full 3D solution.

The reference problem

We consider the typical problem of two plates connected by a bolt. Starting from a simplified model of the assembly using plates elements and a simple connector (typically a beam), our aim is to perform the non-intrusive substitution of the connector by a full 3D model of the bolt. A sketch of the two models is presented in Fig. 1. In what follows lower case Greek characters are used to describe the geometry of the plates whereas capital Greek characters are used for the geometry of the 3D domains: any 3D domain Ω X possesses a plate counterpart ω X , so does any interface Γ Y whose plate counterpart is written γ Y . The target hybrid model, which will be precisely defined in this section, corresponds to plate models connected by a full 3D model of the bolt. At the convergence of the iterations, the solution of the hybrid model is obtained, as illustrated in Fig. 2. Let us note that since the problem is symmetric (in the (e x , e z ) plane) only half of the bolt is computed and shown on the figure.

The plate model

The global plate model is the assembly of two plates, ω inf and ω sup respectively the lower plate and the top plate, and a rigid connector between them ω conn Fig. 1 (left). The plates are 20 mm thick and 280 mm long. The lower plate is 160 mm wide whereas and the top plate is 80 mm wide. In order to simplify the presentation, the plates are assumed to be made out of homogeneous isotropic linear elastic material with Young modulus E = 200 GPa and Poisson ratio ν = 0.3; the handling of orthotropic composite plates is explained in [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF].

U = 0 on γ inf 0 U = U e x on γ sup L (1)
Note that this tension is applied after preload (tightening) was applied to the bolt (see the description of the patch).

A classical Reissner-Mindlin plate formulation is used, and the contact is not taken into account at this stage. The kinematic is given by:

U (m, z) = V (m) + w(m)e z -e z ∧ θ = V (m) + w(m)z -z   θ y (m) -θ x (m) 0   (2)
where m ∈ ω inf ∪ ω sup is a point of the mid-surface of the plates and e z the normal direction (in the thickness). The stress in the plates is described with the eight generalized forces (N tension, M bending, Q shear):

F = (N xx , N yy , N xy , M xx , M yy , M xy , Q x , Q y ) (3)
Since the connector will ultimately be virtually replaced by the 3D model described in next subsection, a very crude connector is chosen: a perfectly rigid beam element surrounded by a rigid region of the size of the bolt (20 mm) in order to avoid localization around one node, see Fig. 3. Beside its simplicity, this connector has the advantage to be very rigid compared to its 3D counterpart, which is a sufficient condition to ensure the convergence of the coupling algorithm [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Guguin | Stratégie non-intrusive de couplage plaque/3D pour la simulation des assemblages de plaques composites stratifiées[END_REF].

The 3D patch

The local patch is designed to replace the connector. Thus, a full 3D representation of the bolt is used with unilateral frictional contact interfaces between each part of the assembly. The dimensions of the bolt are given in Fig. 4(a). The problem being symmetric, only half of the bolt is computed. The material To solve this model, we use the dedicated contact solver COFAST, developed in the LMT-Cachan by L. Champaney [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF]. This software uses the LATIN method [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics: New Approaches and Non-incremental Methods of Calculation[END_REF] where interfaces between subdomains are the support of the contact conditions, see Fig. Note that before coupling with the global problem, preload is applied by enforcing relative displacements between the nut and the screw. It is adjusted to match realistic values of tension in the screw (200 MPa).

Connections between the models

The plates+connector model describes the entirety of the structure and occupies the (2D+1D) domain ω. The 3D model of the bolt occupies what we call the zone of interest Ω I . In the original method [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF], the coarse modeling of ω I = ω ∩ Ω I was simply replaced by Ω I through iterations; in [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF], it was shown that because of the edge effects which affect plate solutions, it is interesting to introduce a zone of transition between the two models. This idea is sketched in Fig. 5: the 3D model is the only one taken into account in the inner zone of interest Ω I ⊂⊂ Ω I , the plate model is the only one taken into account in the outer complement zone ω C ⊂⊂ ω C , there exists an overlap Ω I ∩ ω C also called buffer zone where the two models are equivalent in a certain sense.

Two interfaces are defined:

γ I = ∂Ω I ∩ ω C and γ C = ∂ω C ∩ Ω I .
In order to simplify the example, regular and ruled quadrilateral meshes are used along the interfaces. Thus, γ I and γ C are straight lines between two conforming meshes, see Fig. 6.

To be more specific, let us assume that a finite element discretization is used: u G is the vector of plate unknowns defined on ω, u L is the vector of 3D unknowns defined on Ω I , f ext stands for the generalized forces (they include the effect of imposed displacements), K for the stiffness matrices, and f int for the internal forces (work of the stress field in the finite element subspace).

The global plate model is enriched by an extra loading δ which is only nonzero on the degrees of freedom of the interface γ C (the values of δ are determined by the coupling algorithm). The global plate model equilibrium thus can be written as:

K G u G = f G ext + δ (6)
Note that the assumption of a linear global model simplifies the explanation and speeds up the solving (because the factorization is done only once) but it is not a requirement of the coupling algorithm.

The local 3D equilibrium submitted to boundary conditions inherited from the plate computation can be written as:

f L int (u L ) = -f L ext in Ω I u L = Ru G + w G (u G ) on Γ I (7)
Matrix R represents the classical 3D rigid section movements. It enables to express the 3D Reissner-Mindlin kinematic from the plate degrees of freedom (u G ).

Because the rigid section kinematic is a very coarse hypothesis, we supplement it with warping displacements w G (u G ), proportional to the plate stress state.

To be more precise, the applied displacement on Γ I has the following form:

u L (m, z) = u(m) + w(m)e z -e z ∧ θ Ru G + 8 i=1 µ I i (z)F C i w G (u G )
with (F C i ) the eight generalized forces extracted along the interface γ I (3). (µ I i (z)) is the basis of warping adapted to the stacking and materials which was determined before the computation by solving preliminary problems [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF]. These preliminary problems conducted on a 3D representative volume element of the composite plate enable us to compute distributions of stress and displacement (warping) associated with any solution of Saint-Venant problems. In turn, Saint-Venant problems are associated with specific components of the plate generalized stresses.

From the global and local computations, it is possible to post-process the plate and 3D nodal reactions on the γ C /Γ C interface:

λ C = K C |ω C u G -f C ext|ω C |γ C λ L = -f L int|Ω I (u L ) + f L ext|Ω I |Γ C (8)
Finally the balance of the nodal reactions is enforced on the static plate quantities, thanks to the transpose operator R T which enables to compute the generalized forces associated to the 3D nodal reactions λ L :

λ C + R T λ L = 0 (9)
In spite of the care taken to recover the 3D boundary conditions on Γ I from plate quantities, small edge effects may occur, in particular in the case of stratified plates where suppressing edges effects require high order theories which are not always available in finite element analysis software or which could difficultly be implement in a non-intrusive manner. The role of the buffer zone (the overlap between Ω I and ω C ) is then to dampen these edge effects and make the evaluation of the 3D stress reliable on Γ C . In the example below, the size of the buffer zone is set to b = 5mm, which was sufficient to absorb the local artificial effects. It corresponds to four elements between the two interfaces Γ C and Γ I . The dimension b is represented on Fig. 6, together with dimension L which characterizes the size of the 3D domain of interest Ω I . L corresponds to the size of the part of the 3D domain which is not strictly necessary to represent the bolt correctly but which was inserted as a way to keep the 3D/2D transition away from the zone dominated to by 3D effects.

The non-intrusive iterative algorithm

The system (6,7,8,9) can be interpreted as finding the traction δ to be imposed to the global plate model on the inner interface γ C in order to generate a reaction λ C in balance with the reaction of the inner zone of interest submitted to the recovery of the plate displacement on its boundary Γ I .

Starting from δ 0 = 0, this can be achieved through the following iterations:

1. Run a global plate analysis with extra load δ:

u G n = K G -1 f G ext + δ 2.
Post-process the reaction on γ C :

λ C n = K C u G n|ω C -f C ext |γ C 3.
Recover the 3D displacement on Γ I :

u L n|Γ I = Ru G n|Γ I + w G (u G ) n 4.
Solve the local problem with imposed displacement on Γ I :

f L int (u L n ) + f L ext = 0 u L n|Γ I given 5.
Post-process the local reaction of Γ C :

λ L n = -f L int (u L n ) + f L ext I |Γ C
6. Compute the residual on γ C :

r n = λ C n + R T λ L n 7.
If residual is small enough then exit, else update the extra load:

δ n+1|γ C = δ n|γ C -r n
and go back to 1.

One important point is that the global step 1 and the local step 4 can be processed with different software. For the test-case developed here, Code Aster is used for the global plate problem, and COFAST3D is used for the local 3D contact problem.

In the end the local/global method is a fixed point algorithm similar to a modified Newton algorithm. The rate of convergence of this algorithm can be quite slow. However acceleration techniques can be applied:

• quasi-Newton acceleration like SR1 algorithm [START_REF] Gendre | A two-scale approximation of the schur complement and its use for non-intrusive coupling[END_REF] or BFGS,

• dynamic relaxation like Aitken's algorithm (∆ 2 ) [24],

• mixed boundary conditions [START_REF] Gendre | A two-scale approximation of the schur complement and its use for non-intrusive coupling[END_REF],

• Krylov solvers (only in the linear case) [START_REF] Guguin | Stratégie non-intrusive de couplage plaque/3D pour la simulation des assemblages de plaques composites stratifiées[END_REF].

Based on [START_REF] Guguin | Stratégie non-intrusive de couplage plaque/3D pour la simulation des assemblages de plaques composites stratifiées[END_REF], it appears that, among those methods, SR1 quasi-Newton acceleration (5.3) technique leads to better performance while being simple to implement in a non-intrusive manner.

Note that running one global analysis (step 1) followed by a local reanalysis with given Dirichlet conditions (steps 3-4) without iterations corresponds to the industrialists' practice called submodeling (or sometimes structural zoom). Such an approach is "purely descending" in the sense that there is no feedback from the local computation towards the global scale.

Analysis of the iterative corrections at the global and local levels

In order to evaluate the convergence of the algorithm, the relative norm of the residual at iteration i is defined as:

r rel i = λ C i + R T λ L i 2 λ C i 2 (10)
In this section, the corrections associated with the coupling along the iterations are analyzed for the fixed-point algorithm (i.e. without acceleration techniques) with convergence threshold set to 10 -6 which is clearly enough for mechanical quantities to have converged.

A typical convergence curve is shown Fig. 7(a).This curve shows that the convergence is more or less independent of the level of nonlinearity in the bolt (weak for the first global time step, strong for the last one). This type of result is typical of nonlinear relocalization methods [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: Application to post-buckling analyses[END_REF][START_REF] Hinojosa | Domain decomposition methods with nonlinear localization for the buckling and post-buckling analyses of large structures[END_REF][START_REF] Bordeu | Balancing domain decomposition with nonlinear relocalization: Parallel implementation for laminates[END_REF][START_REF] Allix | A relocalization technique for the multiscale computation of delamination in composite structures[END_REF]. Let us note that, as shown in Fig. 7(b), four global time steps are used for the plate computation, which appears to be sufficient, see Section 5. Anyhow, sub-stepping is used for the computation of the bolt; balance between the local and global model is only reached at the global time steps in Fig. 7(b). For more involved problems, a better control of the global time steps and of the local sub-stepping process could be obtained using error indicators on the time discretization like in [START_REF] Allix | On the control of the load increments for a proper description of multiple delamination in a domain decomposition framework[END_REF]. Remark. With the proposed technique, for each global increment the first iteration corresponds to a classical submodeling approach: this enables us to easily measure the quality of a that approach, which is a question often raised by engineers. In this application, the level of the error of the submodeling approach is about 25%.

Global effect of the correction

Effect of the preload of the bolt

The global effect of the tightening of the bolt cannot be predicted using a plate theory because the associated forces are equal to zero. Nevertheless this tightening has a non-negligible global effect as can be seen in figure Fig. 8.

Additional global effects of the bolt on the plate solution

After the preload of the bolt, the plate is loaded in tension, in four global time steps. The analysis of the number of global increments is presented in Section 5.

The global solution is modified along the iterations to match with the 3D model of the bolt, as can be seen on Fig. 9 which shows the values of the transverse displacement in the mid-section of the plate for the initial plate solution and the corrected one. One can notice that the correction increases with the time steps, and it becomes significant for the third step and very large for the last one. This is of course due to the fact that for the first two time steps the bolt acts more or less as a 3D rigid connector. Whereas the sliding of the bolt becomes significant for the third time step. During the last stage of the loading loosening of the bolt happens, as can be checked from the local analysis of the solution presented in the next sub-section. Let us recall that the initial plate solution also corresponds to what would be obtained by a submodeling approach since submodeling only improves the local reanalyzed area.

Local effect of the correction

In section the solution within the bolt is compared for three different approaches: a reference full 3D simulation, a submodeling approach and the mixed For example the figure Fig. 10 shows the comparison of the σ xz stress field between the submodeling solution, the hybrid solution and the reference one. It can be seen that the solution associated with the non-intrusive approach matches very well the 3D reference solution while the submodeling gives inaccurate results inside the screw where the shear stress is greatly underestimated. On the contrary, the submodeling approach underestimates the sliding especially for the third and fourth global time steps, as can be seen on figure Fig. 13 where the error on the maximal sliding is more than 30%.

To conclude, the coupling approach seems reliable contrarily to the submodeling approach for which the level of error is 20 -25% on global quantities and can be much more for local quantities of interest. The lack of conservatism of the submodeling extends similar results obtained in previous studies on localized plasticity or buckling for example.

Control of some parameters of the method and acceleration technique

As already discussed several parameters can be tuned for the non-intrusive 2D-3D coupling strategy :

• the number of global time steps,

• the choice of the position of the interface in the global model (dimension L in Fig. 6),

• the width of the buffer zone (dimension b in Fig. 6).

The size of the buffer zone has been chosen from our experience on the 2D-3D coupling involving composite plates and orthotropic plies [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF], in order to minimize the problem of artificial edge effects between the two models. Besides the parameters of the nonlinear solver used for the bolt computation have been tuned in order to ensure a high precision according to dedicated papers like [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF].

Influence of the number of global time steps

The number of time steps has been simply chosen on the basis of the comparison of solutions, starting from one global time step and increasing this number. On the figure Fig. 14, a comparison between 1, 2, and 4 time steps is shown. The solutions with 2 and 4 time steps are considered here as quite close, whereas the one using only one time step deviates largely in the middle of the load

Influence of the position of the interface between the 2D and the 3D models

The question of the position of the interface raises in fact the question of the validity of the plate theory with respect to the 3D theory and has been largely discussed in [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF]. From what is known on the validity of the plate and shell theories, in the case of isotropic materials, one expects the 2D-3D model to be a good approximation of the 3D reference, for an interface situated from the bolt at a distance superior to the thickness of the plate. It is therefore interesting to analyze the influence of the position of the interface and to compare the cases of L = 0 mm and L = 15 mm (see Fig. 6). xz , as shown in the Fig. 16. As already analyzed, if for the first two steps, most of the differences are localized around the nut, during the sliding phase the correction mostly concerns the nut itself. The Fig. 17 shows the evolution of the local tangential jump with respect to the prescribed displacement. This is an interesting quantity in order to observe the initiation of the sliding. The 2D-3D models give close predictions for both values of L. They are much more closer to the reference than what is predicted by the simple submodeling.

On the figure Fig. 15(a), it appears, as expected, that a larger number of iterations has to be carried out when the interface is located on the bolt boundaries (L = 0). This effect of the interface location on the convergence rate is corrected when using a SR1 acceleration technique as can be seen in the next subsection.

L=0 L=15

Step 1

Step 2

Step 3

Step 4

∆ σ xz ( M P a )

Figure 16: Error in the local domain for two positions of the interface and at each time step.

Quasi-Newton acceleration

A final component of the non-intrusive substitution is the implementation of acceleration techniques. In that context, the SR1 quasi-Newton acceleration was tested in [START_REF] Gendre | A two-scale approximation of the schur complement and its use for non-intrusive coupling[END_REF], and the Aitken Delta-2 dynamic relaxation method in [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF].

For this study, we only implement SR1 acceleration. The faster with SR1. Moreover the rate of convergence is almost independent from the load step and the position of the interface. This property can be explained by the fact that the corrections induced by the SR1 acceleration technique are adapted to take into account the main differences between the 2D and 3D models.

Conclusions

In this paper, a non-intrusive coupling between plate models and 3D models [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF] has been extended to deal with the precise computations of bolted plates: the plate model with simplified connector was coupled with a full 3D nonlinear model of the bolt. The flexibility of the method was exploited to easily define the coupled model and to use of two dedicated pieces of software: Code Aster for the plate computations, and COFAST3D for the nonlinear computation of the bolt (including many surfaces of friction).

The proposed technique enabled us to analyze the possibilities and limits of the use of plate connectors and submodeling approach in that case. It appears that, even when the bolt response is globally linear, the 3D effects induced by the bolt largely modify the plate solution itself. This shows that a rigid connector is a poor model to describe the connection between two plates. As expected, when important sliding occurs, such modeling becomes irrelevant. Moreover, the submodeling technique may lead to significant local errors and non-conservative results.

Another important feature for future applications concerning the treatment of multiple bolts in interaction, is that using an interface located at the limit of the bolt leads to acceptable results when compared to the reference solution. Works on that type of problem is in progress.

Another issue concerns the modeling of the bolt itself. In practice, some of the parameters of the bolted assembly are not precisely determined, as the preload of the nut or the friction coefficient. Such types of problems have been analyzed by dedicated techniques for the bolt computation [START_REF] Gant | Modeling of the bolted joint behavior variability with the lack of knowledge theory[END_REF][START_REF] Gant | Updating of uncertain joint models using the lack-of-knowledge theory[END_REF], including multiresolution [START_REF] Allix | A new multi-solution approach suitable for structural identification problem[END_REF]. The use of the proposed non-intrusive techniques allows to extend these studies to the case of 2D-3D structural analyses in a straightforward manner. In addition, the use of model reduction techniques for the global model itself, as proposed in [START_REF] Kerfriden | Local/global model order reduction strategy for the simulation of quasi-brittle fracture[END_REF] should lead to a very significant reduction of the computational time.

Figure 1 :

 1 Figure 1: Left, global model: two plates assembled by one 1D connector. Right, local model: full representation of the bolt with surrounding plates.

Figure 2 :

 2 Figure 2: Deformed shape of the targeted coupled problem.

Figure 3 :

 3 Figure 3: The definition of the connector.

  4(b).

Figure 4 :

 4 Figure 4: 3D representation of the bolt.

Figure 5 :

 5 Figure 5: The domain decomposition with overlap and the interfaces.

Figure 6 :

 6 Figure 6: Definition of the size of the local model and of the buffer zone.

Figure 7 :

 7 Figure 7: Convergence of the fixed-point algorithm: (a) Relative norm of the residual (b) Evolution of one component of the global and local nodal reactions.

Figure 8 :

 8 Figure 8: Effect of the tightening (preload): value of the N xx component of the plate stress tensor.

Figure 9 :

 9 Figure 9: Deformation of the plates, before and after iterations.

Figure 10 :

 10 Figure 10: Transverse shear stress disbribution: reference (top), converged 2D-3D solution (middle), submodeling (bottom). We now compare global quantities associated to the bolt (like the amount of dissipation due to sliding or the global axial component of the force acting on the joint) as well as local quantities (like the pointwise values of the sliding) for the three models. On the figures Fig. 11(a) and Fig. 11(b), local stress quantities are extracted on a point defined on Fig. 11(d). It appears that the relative residual provides an efficient indicator for the errors on global quantities. Typically the submodeling leads to 20% of error on the dissipation. Note that

Figure 11 :

 11 Figure 11: Comparison between full 3D, submodeling and 2D-3D approaches.

Figure

  Figure Fig.12presents the sliding obtained by the converged 2D-3D model; it matches the 3D reference (few percents of deviation on the maximal sliding). On the contrary, the submodeling approach underestimates the sliding especially for the third and fourth global time steps, as can be seen on figure Fig.13where the error on the maximal sliding is more than 30%.To conclude, the coupling approach seems reliable contrarily to the submodeling approach for which the level of error is 20 -25% on global quantities and can be much more for local quantities of interest. The lack of conservatism of the submodeling extends similar results obtained in previous studies on localized plasticity or buckling for example.

Figure 12 :

 12 Figure 12: Sliding at each global time step of the converged solution (mm).

Figure 13 :

 13 Figure 13: Maps of the difference of sliding between converged and submodeling solutions at each global time step (mm).

Figure 14 :

 14 Figure 14: Comparison of solutions for different number of global time steps.

Figure 15 :

 15 Figure 15: Comparison of the results for L = 0 mm and L = 15 mm From a global point of view, the final deformed shapes of the corrected solutions in Fig. 15(c) are very close in the common plate domain. They only slightly differ in the local area of interest which are different for the two models. From a local point of view, on figure Fig. 15(b), both cases give close results. This is confirmed when analyzing the difference of the shear stress within the

  Fig. 18 presents convergence plots. The convergence is roughly 3 times

Figure 17 :

 17 Figure 17: Sliding between the two plates near the nut for different configurations.

Figure 18 :

 18 Figure 18: Effect of the SR1 acceleration on the global-local rate of convergence.

Acknowledgements

This work was partially funded by the French National Research Agency as part of project ICARE (ANR-12-MONU-0002-04).

Competing interests

The authors declare that they have no competing interests.


Jozef Haleš

Alice Héliou

Ján Maňuch

Yann Ponty
email: yann.ponty@lix.polytechnique

Ladislav Stacho

Combinatorial RNA Design Designability and Structure-Approximating Algorithm in Watson-Crick and Nussinov-Jacobson Energy Models

Keywords: RNA structure, Inverse combinatorial optimization, String design

We consider the Combinatorial RNA Design problem, a minimal instance of RNA design where one must produce an RNA sequence that adopts a given secondary structure as its minimal free-energy structure. We consider two free-energy models where the contributions of base pairs are additive and independent: the purely combinatorial Watson-Crick model, which only allows equallycontributing A -U and C -G base pairs, and the real-valued Nussinov-Jacobson model, which associates arbitrary energies to A -U, C -G and G -U base pairs.

We first provide a complete characterization of designable structures using restricted alphabets and, in the four-letter alphabet, provide a complete characterization for designable structures without unpaired bases. When unpaired bases are allowed, we characterize extensive classes of (non-)designable structures, and prove the closure of the set of designable structures under the stutter operation. Membership of a given structure to any of the classes can be tested in Θ(n) time, including the generation of a solution sequence for positive instances.

Finally, we consider a structure-approximating relaxation of the design, and provide a Θ(n) algorithm which, given a structure S that avoids two trivially non-designable motifs, transforms S into a designable structure constructively by adding at most one base-pair to each of its stems.

Introduction

RiboNucleic Acids (RNAs) are biomolecules which act in almost every aspect of cellular life, and can be abstracted as a sequence of nucleotides, i.e., a string over the alphabet {A, U, C, G}. Due to their versatility, and the specificity of their interactions, they are increasingly being used as therapeutic agents [START_REF] Wu | RNAi therapies: Drugging the undruggable[END_REF], and as building blocks for the emerging field of synthetic biology [START_REF] Rodrigo | Full design automation of multi-state RNA devices to program gene expression using energy-based optimization[END_REF][START_REF] Takahashi | A modular strategy for engineering orthogonal chimeric RNA transcription regulators[END_REF]. A substantial proportion of the functional roles played by RNA rely on interactions with other molecules to activate/repress dynamical properties of some biological system, and ultimately require the adoption of a specific conformation. Accordingly, RNA bioinformatics has dedicated much effort to developing energy models [START_REF] Mathews | Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure[END_REF][START_REF] Turner | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF] and algorithms [START_REF] Nussinov | Fast algorithm for predicting the secondary structure of singlestranded RNA[END_REF][START_REF] Zuker | Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information[END_REF] to predict the secondary structure of RNA, a combinatorial description of the conformation adopted by an RNA which only retains interacting positions, or base pairs. Historically, structure prediction has been addressed as an optimization problem, whose expected output is a secondary structure which minimizes some notion of free-energy [START_REF] Nussinov | Fast algorithm for predicting the secondary structure of singlestranded RNA[END_REF][START_REF] Zuker | Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information[END_REF]. The performances of the RNA folding prediction problem have now reached a point where in silico predictions are generally considered reliable [START_REF] Mathews | Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure[END_REF], allowing for large scale studies and fueling the discovery of an increasing number of functional families [START_REF] Griffiths-Jones | RFAM: an RNA family database[END_REF].

Due to the existence of expressive, yet tractable, energy models and conformational spaces, coupled with promising applications in multiple fields (pharmaceutical research, natural computing, biochemistry. . . ), a wide array of computational methods [START_REF] Hofacker | Fast folding and comparison of RNA secondary structures[END_REF][START_REF] Busch | INFO-RNA-a fast approach to inverse RNA folding[END_REF][START_REF] Aguirre-Hernández | Computational RNA secondary structure design: empirical complexity and improved methods[END_REF][START_REF] Dai | RNADesign: Local search for RNA secondary structure design[END_REF][START_REF] Avihoo | RNAexinv: An extended inverse RNA folding from shape and physical attributes to sequences[END_REF][START_REF] Taneda | MODENA: a multi-objective RNA inverse folding[END_REF][START_REF] Zadeh | Nucleic acid sequence design via efficient ensemble defect optimization[END_REF][START_REF] Levin | A global sampling approach to designing and reengineering RNA secondary structures[END_REF][START_REF] Lyngsø | FRNAkenstein: multiple target inverse RNA folding[END_REF][START_REF] Garcia-Martin | RNAiFOLD: a constraint programming algorithm for RNA inverse folding and molecular design[END_REF][START_REF] Höner Zu Siederdissen | Computational design of RNAs with complex energy landscapes[END_REF][START_REF] Reinharz | A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotide distribution[END_REF][START_REF] Zhou | Flexible RNA design under structure and sequence constraints using formal languages[END_REF][START_REF] Bau | RNA design by program inversion via SAT solving[END_REF][START_REF] Esmaili-Taheri | Evolutionary solution for the RNA design problem[END_REF] has been proposed to tackle the natural inverse version of the structure prediction, the RNA design problem. In this problem, one attempts to perform the in silico synthesis of artificial RNA sequences, performing a predefined biological function in vitro or in vivo. Given the prevalence of structure in the function of an RNA, one of the foremost goal of RNA design (sometimes named inverse folding in the literature) is to ensure that the designed sequence folds into a predefined secondary structure, preferentially to any alternative structure. In other words, the chosen conformation should not be challenged by alternative stable structures having similar or lower free-energy.

Despite a rich, fast-growing, body of literature dedicated to the problem, there is currently no exact polynomial-time algorithm for the problem. Moreover, the complexity of the problem remains open (see Section 5 for a discussion). It can be argued that this situation, quite exceptional in the field of computational biology, partly stems from the intricacies of the Turner free-energy model [START_REF] Turner | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF] which associates experimentally-determined energy contributions to ∼ 2.4 × 10 4 structure/sequence motifs. This motivates a reductionist approach, where one studies an idealized version of the RNA design problem, lending itself to algorithmic intuitions, while hopefully retaining the presumed difficulty of the original problem and provides intuitions for future studies of the problem under more sophisticated energy models.

In this work, we introduce the Combinatorial RNA Design problem, a minimal instance of the RNA design problem which aims at finding a sequence that admits the target structure as its unique base pair maximizing structure. After this short introduction, Section 2 states definitions and problems. In Section 3, we state our main results and prove them in Section 4, including an extended weighted version that allows additional types of base pairs. Finally, we conclude in Section 5 with some remarks, open problems and future extensions of this work.

G C A G G A G U C U A G C G A U G C U A G U C A G C U A G C U C A U A A U G A A U U A G G C U A C G A C U A G C G C U G A G A C C C U U 1
G C A G G A G U C U A G C G A U G C U A G U C A G C U A G C U C A U A A U G A A U U A G G C U A C G A C U A G C G C U G A G A C C C U U
(l i , r i )} p i=1 ⊂ [1, n] 2 such that: -∀i ∈ [1, p], l i < r i ;
-Each position is involved in at most one base pair, i.e. ∀i = j ∈ [1, p], l i = l j , l i = r j , r i = r j ; -Base pairs are pairwise non-crossing,i.e. i, j ∈ [1, p], l i < l j < r i < r j , .

The set of unpaired positions U S of a secondary structure S = (n, P ) is the set of indices k ∈ [1, n] that are not involved in any of the base pairs in S. A structure S is called saturated if and only if its positions are all paired, i.e. iff U S = ∅. Conversely, a structure S is empty when none of its positions is paired, i.e. iff S = (n, ∅).

The set of all secondary structures is denoted by S, and its restriction to structures of length n by S n . Secondary structures are typically expressed using a variety of equivalent representations, illustrated by Figure 1 and formally defined further in this section.

Given a sequence w and a structure S = (|w|, P ), let

u i = ε if i ∈ U S and u i = w i ,
• • u i |, |u 1 • • • u j |) | (i, j) ∈ P })
, where w is any sequence of length n and u i 's are defined as above.

A maximal subset B = {(i, j), (i+1, j -1), . . . , (i+ -1, j -+1)} of P for some integer i, j, is called a band (also referred to as helix or stem in related works) of size = |B|, of S = (n, P ). Note that every base pair belongs to exactly one band. In other words, the base-pairs of a secondary structure can be unambiguously partitioned into a set of bands. Dot-parentheses notation. A well-parenthesized sequence s ∈ {(, ), .} * can be used to represent a secondary structure. There exists a one-to-one correspondence between secondary structures and such well-parenthesized sequences: any base pair (l, r) ∈ P becomes a pair of corresponding opening and closing parentheses in s at position l and r respectively (s l = ( and s r = )), and any unpaired position i corresponds to a dot (s i = .). This representation is illustrated by Figure 1.c. A concatenation of two structures S and S , denoted by S.S or simply SS wherever unambiguous, is the structure corresponding to the well-parenthesized sequence obtained by concatenating the well-parenthesized sequences of S and S . k-stutter. The k-stutter of a sequence s, denoted by s [k] is the result of an independent copy k-times of each of the characters in s. For instance, the 3-stutter of a sequence AUUC is AAAUUUUUUCCC. This operation also applies to an RNA structure S, and S [k] denotes the RNA structure obtained by applying the usual k-stutter to the dot-parentheses representation of S.

Tree representation. Alternatively, the tree representation, denoted by T S , for S = (n, P ) is a rooted ordered tree whose vertex set V S consists of intervals [l, r] for any base pair (l, r) ∈ P , and [k, k] for every k ∈ U S . A virtual root [0, n + 1] is added for convenience. Any node labeled by [k, k] is called unpaired, and any other node (including the virtual root) are considered as paired. The children of an interval I ∈ V S are the maximal proper subintervals I ∈ V S of I ordered by the left points of the intervals. The degree of a vertex I ∈ V S is the total number of its paired neighbors, including its parent (if any). We denote by D(S) the maximal degree of nodes in T S . Figure 1.d shows the tree representation of a typical secondary structure.

Proper, greedy and separated coloring of the tree representation. Consider the tree representation T S of structure S. A coloring of T S associates a color, chosen among black, white, or gray, to each paired node of T S that is different from the root. This coloring is called proper if: i) Each node has at most one black child, at most one white child, and at most two grey children; ii) Any c-colored node has at most one c-colored child; iii) Black nodes shall not have a white child, and white nodes shall not have a black child.

A greedy coloring of T S is the coloring obtained by recursive application of the following rule starting from the root and continuing towards leaves: if the node is black, color the first paired child black and the remaining paired children gray, if the node is white, color the first paired child white and the remaining paired children gray, otherwise (the gray node or the root), color the first paired child black, second white and the remaining paired children gray. It is easy to check that if the degree of each node is at most four then the greedy coloring is a proper coloring.

Given a proper coloring of T S , let the level of each node be the number of black nodes minus the number of white nodes on the path from this node to the root. A proper coloring is called separated if the two sets of levels, associated with gray and unpaired nodes respectively, do not overlap.

Statement of the generic RNA design problem

Consider an energy model M, which associates a free-energy E M (w, S) ∈ R -∪ {+∞} to each secondary structure S ∈ S |w| for a given RNA sequence w. The minimum free-energy (MFE) structure prediction problem is typically defined as follows:

RNA-FOLD M problem Input: RNA sequence w Output: S M (w) := argmin S ∈S |w| E M (w, S ) ,
where argmin returns a single structure S M (w), arbitrarily chosen amonst those having minimum free-energy.

The existence of alternative competing structures, i.e. one or several secondary structure(s) having (almost) minimal free-energy for a given RNA, impacts the efficacy of the folding process. The detection of such situations is therefore of interest, and can be rephrased as the following problem:

UNIQUE-FOLD M problem

Input: Sequence w + Energy distance ∆ > 0 Output: True if, for every S ∈ S |w| \ {S M (w)}, one has:

E M (w, S ) ≥ E M (w, S M (w)) + ∆ .

False otherwise.

We can now define the combinatorial RNA Design problem as:

RNA-DESIGN M,Σ problem Input: Secondary structure S + Energy distance ∆ > 0 Output: RNA sequence w ∈ Σ -called an (M, Σ, ∆)-design for S -such that:

RNA-FOLD M (w) = S and UNIQUE-FOLD M (w, ∆),
or ∅ if no such sequence exists.

The structures for which there exists an (M, Σ, ∆)-design are called (M, Σ, ∆)designable. Let Designable(M, Σ, ∆) be the set of all such structures. If it is clear from the context, we will usually drop M, Σ and/or ∆.

✭ ✭ ✳ ✮ ✭ ✳ ✳ ✮ ✮ • • ❆ ❈ ❆ • • ❯ ❈ ❆ ❈ ❆ • • ❯ ❯ ❈ ❯ ❛✳ ❚❛r❣❡t s❡❝✳ str✳ S ❜✳ ■♥✈❛•✐❞ s❡q✉❡♥❝❡ ❢♦r S ❝✳ ❉❡s✐❣♥ ❢♦r S
Fig. 2 The combinatorial RNA design problem: Starting from a secondary structure S (a.), our goal is to design an RNA sequence which uniquely folds, with maximum number of base pairs, into S. The sequence proposed in b. is invalid due to the existence of an alternative structure (lower half-plane, red) having the same number of base pairs as S. The right-most sequence (c.) is a design for S.

Base pair sum energy models

In this work, we will consider two types of base pair sum energy models, where the free-energy of a structure is simply obtained by sum, over all base pairs, individual independent contributions associated with each pair.

Definition 1 (Base pair sum energy model M) Let w be an RNA sequence and S a secondary structure in S |w| . Then

E M (w, S) = (l,r)∈S E M (w l , w r ) ,
where E M (x, y) is the energy induced by a base pair x -y.

To define a model of interest, it is sufficient to specify the energies of base pairs. First, we consider a minimal energy model, named Watson-Crick model due to its similarity with the DNA base-pairing rules. The model is purely combinatorial, as it associates a homogenous -1 energy contribution to each valid base-pair, and only allows G -C and A -U to pair.

Definition 2 (Watson-Crick energy model W)

E W (x, y) = -1 if {x, y} = {G, C} or {x, y} = {A, U} +∞ otherwise.
A more general model, named the Nussinov-Jacobson model, allows G -U basepairs to occur, and associates arbitrary weights to the base pairs depending on their content. It is named after the authors of the first polynomial-time algorithm for predicting the MFE under a similar energy model [START_REF] Nussinov | Fast algorithm for predicting the secondary structure of singlestranded RNA[END_REF].

Definition 3 (Nussinov-Jacobson energy model N ) E N (x, y) =          α if {x, y} = {G, C} β if {x, y} = {A, U} γ if {x, y} = {G, U} +∞ otherwise.
where α, β, γ ∈ R -and α, β < γ.

Note that the last condition of the above definition is typically satisfied by empirical estimates of base pair energies. Namely, G -U base pairs, also named Wobble base pairs, are much weaker than its alternatives. They are mediated by a single hydrogen bond, as opposed to 2 (resp. 3) bonds for A -U (resp. G -C).

We say that the structure is compatible with a sequence w, according respect to an energy model M, if E M (w, S) < +∞.

It is worth noting that minimizing E W (w, S) is equivalent to maximizing |S|, thus RNA-FOLD W is a classic base pair maximization problem. Moreover, both RNA-FOLD W and RNA-FOLD N can be solved in polynomial time using dynamic programming, historically in O(n 3 ) complexity [START_REF] Nussinov | Fast algorithm for predicting the secondary structure of singlestranded RNA[END_REF], or in O(n 3 / log(n)) current best time complexity [START_REF] Frid | A simple, practical and complete o(n 3 / log n)-time algorithm for RNA folding using the Four-Russians speedup[END_REF]. A backtracking procedure reconstructs the structure having minimal energy, and can be easily adapted to provide a Θ(n 3 ) algorithm for the UNIQUE-FOLD M problem.

Statement of the results

In this section, we characterize sets of secondary structures which can or cannot be designed using (a subset of) {A, C, G, U}, a (restricted) set of base pairs and a desired energy distance ∆ within an energy model M. The proofs of our statements are largely interconnected, and have been regrouped in Section 4.

First, let us remark that the empty secondary structures are the only ones that are designable for arbitrary large energy distances ∆.

Theorem 1 For any M ∈ {W, N }, and any energy distance ∆ such that

∆ > -E M (X, Y), ∀X, Y ∈ {A, C, G, U} 2 , (1)
only the empty secondary structures are designable.

Proof For any non-empty secondary structure S having energy E on some sequence w, removing any base pair X -Y from S yields an alternative structure S whose energy is

E = E -E M (X, Y) < E +∆.
In other words, S is a competing structure at distance less than ∆ of S, i.e. w is not a valid ∆-design for S. Moreover, any empty structure of length n is designable. Indeed, none of the models allows for pairs of the form X -X, X ∈ {A, C, G, U}, so any sequence X n admits the empty structure, having 0 energy, as its only secondary structure having finite free-energy, i.e. X n is a design for the empty structure for any finite ∆ > 0.

Watson-Crick model W (∆ = 1)

We provide (partial) characterizations for the sets Designable(Σ) of designable structures over partial alphabets Σ in the W model. From Theorem 1, combined with the purely combinatorial nature of the energy model, we observe that nontrivial structures are designable only when ∆ ∈ (0, 1], and that the set of designable structures is unaffected by the precise value of ∆ on the segment. Therefore we set ∆ = 1 without loss of generality, and consider the structures whose designed sequences lose at least one base-pair when forming alternative structures. We obtain the following meta-theorem.

❛ ❜

Fig. 4 Examples of two saturated structures, one satisfying conditions of R4 (a.) and one that does not (b.). The structure on the left has pair degree 4 and is designable, and the structure on the right has pair degree 5 and is not designable.

Theorem 2 In the Watson-Crick energy model W, and assuming an energy distance of ∆ = 1, results R1 through R8 hold.

Let Σ c,u be an alphabet with c pairs of complementary bases and u bases without a complementary base. Without loss of generality in the W model, we will assume that Σ 1,0 = {G, C} and Σ 1,1 = {G, C, A}.

Designability over restricted alphabets.

First, we provide a complete characterization of designable (secondary) structures using an alphabet Σ of restricted cardinality:

R1 For every u ∈ N + , Designable(Σ 0,u ) = {(n, ∅) | ∀n ∈ N}; R2 Designable(Σ 1,0 ) = {S ∈ S | S is saturated and D(S) ≤ 2} ∪ {(n, ∅) | ∀n ∈ N}; R3 Designable(Σ 1,1 ) = {S ∈ S | D(S) ≤ 2}.
Figure 3 shows examples of secondary structures satisfying conditions of these three results.

Designability over the complete alphabet

Σ 2,0 = {A, U, C, G}.
We first characterize the set of designable saturated structures, i.e. structures whose positions are all paired to some other position.

R4 {S ∈ Designable(Σ 2,0 ) | S is saturated} = {S ∈ S | D(S) ≤
Ŝ: ✭ ✳ ✳ ✮ ✭ ✳ ✳ ✮ ✭ ✳ ✳ ✮ S: ✭ ✭ ✳ ✳ ✮ ✮ ✭ ✭ ✳ ✳ ✮ ✮ ✭ ✭ ✳ ✳ ✮ ✮ • • ❯ ❯ ❈ ❈ ❈ ❈ ❆ ❆ • • ❆ ❆ • • ❯ ❯ Fig. 6
An example of non-designable (left) and designable structure (right).

Figure 4 shows two saturated secondary structures one with pair degree 4 and one with pair degree 5. By R4, the former is Σ 2,0 -designable, while latter is not.

When unpaired positions are allowed in the target structure, our characterization is only partial: R5 Let m 5 represent "a node having degree more than four", and m 3 • be "a node having one or more unpaired children, and degree greater than two" (cf. Figure 5), then

Designable(Σ 2,0 ) ∩ {S ∈ S | S contains m 5 or m 3 • } = ∅ ;
R6 Let Sep be the set of structures for which there exists a separated (proper) coloring of the tree representation, then Sep ⊂ Designable(Σ 2,0 ); R7 The set of Σ 2,0 -designable structures is closed under the k-stutter operations:

∀S ∈ S, ∀k ∈ N + : S ∈ Designable(Σ 2,0 ) =⇒ S [k] ∈ Designable(Σ 2,0 ) .
We note however that reverse implication is not true: S [k] ∈ Designable(Σ 2,0 ) does not imply that S ∈ Designable(Σ 2,0 ). For instance, it is easily verified that Ŝ [START_REF] Avihoo | RNAexinv: An extended inverse RNA folding from shape and physical attributes to sequences[END_REF] is Σ 2,0 -designable, while Ŝ is not, as shown in Figure 6. Membership to the classes described in R1-R5 can be tested by trivial linear-time algorithms. These algorithms can also be easily adapted into linear-time algorithms for the production of a concrete design, thereby offering partial solutions to the RNA-DESIGN M,Σ problem.

Structure-approximating algorithm.

Unfortunately, avoiding m 5 and m 3 • , while necessary, is generally not sufficient to ensure designability. For instance, consider Ŝ in Figure 6 clearly does not contain

Ŝ: ✭ ✳ ✳ ✮ ✭ ✳ ✳ ✮ ✭ ✳ ✳ ✮ S: ✭ ✭ ✳ ✳ ✮ ✮ ✭ ✭ ✳ ✳ ✮ ✮ ✭ ✭ ✳ ✳ ✮ ✮ • • ❯ ❯ ❈ ❈ ❈ ❈ ❆ ❆ • • ❆ ❆ • • ❯ ❯ Fig. 7
Application of the structure-approximating algorithm to the non-designable structure Ŝ in Fig. 6: A base pair (circled black node) is inserted in the greedily colored tree, offsetting the levels of white and unpaired nodes (crosses) to even and odd levels respectively, so that the resulting tree is proper/separated, representing a designable structure.

m 5 or m 3 • , yet it cannot be designed. In such cases, unwanted interactions can be somehow penalized by duplicating some base pairs. For instance, duplicating a single base pair in Ŝ yields a Σ 2,0 -designable structure S, as shown by Figure 7.

R8 Any structure S avoiding m 5 and m 3 • can be transformed in Θ(n) time into a Σ 2,0 -designable structure S . This is done by duplicating a subset of the base pairs of S, at most one per band, such that the greedy coloring of the resulting structure is proper and separated, as illustrated by Figure 7.

Nussinov-Jacobson energy model N (∆ ≤ min(|α|, |β|)])

We consider the validity of the above results in the Nussinov-Jacobson model. Note that the consideration of G -U base pairs, by loosening the notion of complementarity, forces us to abuse our notation for Σ i,j . Namely, Σ 2,0 is taken to represent the full alphabet, even though it now strictly allows three types of base pairs.

Theorem 3 For any ∆ ∈ (0, min(|α|, |β|)], statements R1 through R8 hold in the Nussinov-Jacobson model.

Proofs

Watson-Crick model (∆ = 1)

R1 is trivial since, in the absence of complementary letters, empty structures are the only one whose energy is not infinite.

Theorem 4 (⇒ Result R2, R3 and R4) A saturated secondary structure S is Σ c,0 -designable if and only if D(S) ≤ 2c.

Proof First, we will show that the degree condition is necessary. Assume to the contrary that D(S) > 2c and S has a design w. Let [a, b] be a vertex with degree

d ≥ 2c + 1 in T S . Let {[l i , r i ]} d i=1 be the (paired) children of [a, b] and the node [a, b] if [a, b] is not the root. Let L i = l i and R i = r i if [l i , r i ] is a child of [a, b], and
L i = r i and R i = l i if it is [a, b].
Then among bases w L 1 , . . . , w L d must be a pair of repeated letters. Let w L i = w L j be such a pair with L i < L j . It is easy to check that S \ {(l i , r i ), (l j , r j )} ∪ {(L i , R j ), (R i , L j )} is a structure compatible with w with the same number of base pairs as S, a contradiction with the assumption that w is a design for S.

To show that the degree condition is also sufficient, we need further definitions and claims. First, we say that a sequence w ∈ Σ * is saturable if there is a saturated structure compatible with w. Note that the concatenation of two saturable sequences is also saturable. Then the following claim characterizes the cases when a saturable sequence can be split into saturable sequences. Claim 4.1 Let w = uv be a saturable sequence of length k. If u is saturable, then so is v.

Proof Consider a saturated structure S compatible with sequence w and a saturated structure S u compatible with u. We will construct a saturated structure S v compatible with v.

Consider a graph G with vertex set {1, . . . , k} and edge set defined by pairs in S ∪ S u . Obviously, this graph is a collection of alternating paths (alternating between pairs from S and from S u , starting and ending with positions in v) and alternating cyclic paths, and it has a planar embedding such that all vertices lie on a line in their order: pairs in S are drawn as non-crossing arcs above the line and pairs in S u as non-crossing arcs below the line. Note that every position in v is an end-point of exactly one path in the collection.

Define set of base pairs S v by pairing the end-points of the paths in G, cf. Figure 8. We will show that S v is a structure. Consider a graph G constructed by adding pairs in S v to G. This graph is a collection of cyclic paths. Consider an embedding of G into plane that extends the planar embedding of G by adding arcs corresponding to the pairs in S v below the line containing all the vertices. If two base pairs b, b ∈ S v cross then the cyclic path containing b and the cyclic path containing b intersect in exactly one point. By Jordan's curve theorem, this is a contradiction. It follows that S v is a saturated structure, and hence v is also saturable.

We define w to be an atomic saturable sequence if no proper prefix of w is saturable. Every saturated structure compatible with an atomic saturable sequence w contains the base pair (1, |w|), since otherwise it contains the base pair (1, j), with j < |w|, and consequently w [1,j] is a saturable proper prefix of w. On the other hand, by Claim 4.1, if every saturated structure compatible with w contains the pair (1, |w|), then w is an atomic saturable sequence. A design w that is also an atomic saturable sequence will be called an atomic saturable design. A concatenation of two or more atomic saturable designs is obviously not an atomic saturable sequence and it is not necessarily a design. However, we have the following claims.

Claim 4.2

The concatenation of t atomic saturable designs w 1 , . . . , w t for structures S 1 , . . . , S t , such that w

i 1 = w j 1 , ∀1 ≤ i < j ≤ t, is a design for the concate- nated (saturated) structure S = S 1 • • • S t .
Proof Assume that W := w 1 • • • w t is not a design, then there exist a saturated structure S = S for W . We show that positing such an alternative structure leads to a contradiction. Recall that each S i is saturated and contains a pair (1, |w i |). If S pairs the first and last letters in each w i , i ∈ [1, t], then S = S since each w i -While there exists an unprocessed internal node [a, b] whose parent has been processed (if there is no such node, stop and return w). Let {[l i , r i ]} d i=1 be the children of [a, b]. Assign to each w l i , w r i complementary bases such that ∀1 ≤ i ≤ d : w l i = w b and ∀1 ≤ i < j ≤ d : w l i = w l j .

Note that since the alphabet contains c pairs of complementary bases, the assignment at each step of the algorithm is possible. We will show that the returned sequence w is a design for S. We will show by tree induction on the size subtrees that w i • • • w j is an atomic saturable design for every internal node [i, j]. It is easy to check that this is satisfied at the leaves. Consider an internal node u. By the induction hypothesis, sequences for each child subtree of u are atomic saturable designs. Furthermore, by the choice of bases at children nodes of u, all assumptions of Claim 4.3 are satisfied, hence, the sequence for node u is also an atomic saturable design. The claim holds. Finally, we can apply Claim 4.2 at the root, which yields that w is a design.

Corollary 5 (Result R2) A structure S is Σ 1,0 -designable if and only if it does not contain any base pairs, or it is saturated and D(S) ≤ 2.

Proof If S contains a base pair and an unpaired position, then it can be easily checked that S is not Σ 1,0 -designable. Hence, any Σ 1,0 -designable structure is either empty, and trivially designable using a single letter, or saturated. In the latter case, by Theorem 4, we know that designable structures are exactly those that are saturated, and such that D(S) ≤ 2. The claim follows.

Corollary 6 (Result R3) A structure S is Σ 1,1 -designable if and only if D(S) ≤ 2.
Proof First, suppose S is Σ 1,1 -designable and let w be a design for Conversely, suppose that D(S) ≤ 2. Construct a design for S as follows. Since Paired(S) is saturated, by Theorem 4, there is a design w for Paired(S) over Σ 1,0 ⊂ Σ 1,1 . Construct w from w by inserting the base without a complementary base at every unpaired position of S. Let S be an MFE structure for w. Obviously, all unpaired positions in S are also unpaired in S . We must have Paired(S ) = Paired(S), otherwise we have an alternative structure for w, a contradiction. Hence, S = S, i.e., w is a design for S.

Result

Let L i = l i and R i = r i if [l i , r i ] is a child ❆ ❇ ✲ ✲ ✲ ✲ ✲ ❆ ❇ ❆ ❇ ✲ ✲ ✲ ❆ ❇ ❈ ❉ ✲ ✲ ✲ ❆ ❇ ❈ ❉ ❇ ✲ ✲ ❆ ❇ ❈ ❉ ❆ ✲ ✲ ❆ ❇ ❈ ❉ ❈ ✲ ✲ ❆ ❇ ❈ ❉ ❉ ✲ ✲ ❆ ❇ ❇ ❆ ✲ ✲ ✲ ❆ ❇ ❇ ❆ ❇ ✲ ✲ ❆ ❇ ❇ ❆ ❈ ✲ ✲ ❆ ❇ ❇ ❆ ❆ ✲ ✲ ❆ ❇ ❇ ❆ ❈ ❆ ❇ ❆ ❇ ❇ ❆ ❈ ❇ ❆ ❆ ❇ ❇ ❆ ❈ ❈ ❉ ❆ ❇ ❇ ❆ ❈ ❉ ❈
L i = r i and R i = l i if it is [a, b]. If among bases w L 1 , .
. . , w L 3 there is a pair of repeated letters, then we can construct an alternative MFE structure for w (see the first paragraph in the proof of Theorem 4). Assume that these three bases are different. Then for some i = 1, 2, 3, w u equals either w l i or w r i , say it equals w l i . Then S \ {(l i , r i )} ∪ {(u, r i )} is an MFE structure for S, a contradiction with the assumption that w is a design for S.

Theorem 8 (Result R6) If the tree representation of a structure S admits a separated coloring then S is Σ 2,0 -designable.

Proof Given a sequence w, we define the level L(i) of position i as

L(i) = |w [1,i] | G - |w [1,i] | C .
Claim 8.1 Consider any structure S compatible with sequence w that contains some A -U base pair between positions at different levels, then there exists a position G or C that is left unpaired in S.

Proof Assume that the A -U base pair occurs at position (a, b), and note that the bases of the substring w [a+1,b-1] can only base pair among themselves without introducing crossings. We will show that G's and C's are not balanced in this substring. Since

w b ∈ {A, U}, L(b) = L(b -1)
. Hence, by the definition of L, we have that

|w [a+1,b-1] | G -|w [a+1,b-1] | C = L(b -1) -L(a) = L(b) -L(a) = 0 .
Therefore, at least one G or C in the substring remains unpaired in this structure.

❈ ❯

• ❆ Fig. 10 The compatibility graphs of the Watson-Crick (blue edges) and the Nussinov-Jacobson (blue and red edges) energy models are bipartite.

Consider a separated coloring of the tree representation of S. We will use this coloring to construct a design w for S, by specifying a nucleotide at each position of w. First, for each unpaired position i, set w i = A. Second, apply a modified version of the algorithm described in the proof of the Theorem 4 to set the bases of paired positions in which black nodes are assigned to base pair G -C, white nodes to C -G and grey nodes to A -U or U -A. The algorithm ignores unpaired nodes in the tree representation of S. Since the coloring is proper such assignment is always possible at every step of the algorithm. We claim that for any node [i, j] (paired or unpaired), the level of position i is the same as the level of the node [i, j]. To verify this, observe that the substring of w corresponding to any subtree has the same number of G's and C's. Hence, for any node [i, j], the level of position i depends only on nodes on the path from this node to the root. It is easy to check that the level of i is equal to the level of the node. Note that if [i, j] is a grey node then the level of position j is the same as the level of i, i.e., the same as the level of [i, j].

We will show that the constructed w is a design for S. Since all C's and U's of w are paired in S, S is an MFE structure for w. We need to show that it is the only MFE structure for w. Consider an MFE structure S for w different from S. Since w has the same number of G's and C's, S must pair all G's, C's and U's of w. We will show that all unpaired positions in S are also unpaired in S . Assume to the contrary that position i is unpaired in S, but it is paired to j in S . We must have w i = A and w j = U. Since the coloring is separated, the unpaired node [i, i] has a different level than the grey node containing j, and hence, the level of i is different from the level of j. It follows by Claim 8.1 that some G or C is unpaired in S , a contradiction. Consider the paired restrictions of S, S and w. Both Paired(S) and Paired(S ) are saturated and compatible with Paired(w, S) and they are different since S and S are different and agree on the unpaired positions. Furthermore, Paired(w, S) can be produced by the algorithm described in the proof of Theorem 4 for the input structure Paired(S), and hence, by Theorem 4, Paired(w, S) is a design for Paired(S), which contradicts the existence of Paired(S ). Hence, w is a design for S.

Next, we show the closure of the set of designable structures under the kstutter operation. To that purpose, we introduce the compatibility graph of an energy model M, whose vertices are the four nucleotides {A, C, G, U}, and whose edges correspond to valid base-pairs in M, i.e. having finitely-valued contributions.

Definition 4 (Bipartite energy model) An energy model M is bipartite if and only if its compatibility graph is bipartite.

The Watson-Crick energy model W and the Nussinov-Jacobson energy model N are bipartite, as can be seen in Figure 10.

Theorem 9 (Result R7) For any bipartite energy model M and any energy distance ∆, if w is a ∆-design for a structure S, then for any integer k ≥ 1, w [k] is also a ∆-design for S [k] .

Proof Consider a designable structure S and let w = w 1 • • • w n be a design for S. We will show that w [k] is a design for S [k] . Let us use the i th block in S [k] (resp. w [k] ) as a shorthand for the subset [1 + i • k, 1 + (i + 1) • k) of its positions. Note that the positions involved in the i th block of S [k] correspond to the i th position in S (resp. w).

Consider a valid structure S = S [k] for w [k] . Define an interaction multigraph I(S ) = (V I(S ) , E I(S ) ) of S as follows: the vertex set V I(S ) is the set of positions {1, . . . , n} in w, and there are as many edges between i and j in I(S ) as there are base-pairs in S between the i th and j th blocks. Clearly, I(S ) is a multigraph of maximal degree k. Moreover any edge between the i th and j th block in I(S ) corresponds to a valid base-pair (i, j) for w. Therefore the sequence of nucleotides read along any path in I(S ) must constitute a valid path in the compatibility graph. Since the energy model is bipartite, then any cyclic path cannot have odd length, and so I(S ) is also bipartite.

Since I(S ) is a bipartite multigraph of maximal degree k, then by König's theorem [START_REF] Kőnig | Gráfok és alkalmazásuk a determinánsok és a halmazok elméletére[END_REF] it is k edge-colorable (see [23, page 52] for an English version of the proof). In other words, we can color the base-pairs of S , using less than k colors, such that each block in S is involved in at most one base-pair of each color. Therefore we can partition the base-pairs of S into k structures S 1 , S 2 , • • • , S k that are compatible with w. Note that the base-pairs of S i are pairwise non-crossing since S itself is non-crossing.

The sequence w is a design for S, thus one has

E W (w, S i ) ≥ E W (w, S), for every 1 ≤ i < k.
Moreover, one has S = S [k] so there exists a structure S j such that S j = S, and therefore E W (w, S j ) ≥ E W (w, S) + ∆. It follows that

E W (w [k] , S ) = k-1 i=0 E W (w, S i ) ≥ k • E W (w, S) + ∆ = E W (w [k]
, S [k] ) + ∆.

We conclude that S [k] is the sole MFE structure for w [k] , and has energy at least ∆ less than its foremost competitor, so w [k] is a ∆-design for S [k] .

Result R7 immediately follows from Theorem 9, by reminding the bipartite nature of the Watson-Crick energy model.

Theorem 10 (Result R8) Each structure S without m 5 and m 3 • can be transformed into a Σ 2,0 -designable structure S by inflating a subset of its base pairs (at most one per band). Furthermore, this transformation can be done in Θ(n) time.

Proof We start with the greedy coloring of T S . Since S does not contain m 5 and m 3 • , it is a proper coloring and there is no node having both a grey child and an unpaired child. We will insert base pairs within S so that the grey nodes and any unpaired node end up at levels of different parities. If the root has a grey child, assign even parity to the grey nodes, otherwise (if the root has an unpaired child, or no grey and no unpaired children), assign even parity to the unpaired nodes. Now we proceed from the children of the root towards leaves adjusting parity level for grey and unpaired nodes to keep one type even and the other one odd. We repeatedly apply the following simple operation on T S : If the node N does not match its intended parity level. Denote N P the parent of N (N P is not the root as all children of the root already have the correct parity level) and N P P the parent of N P . Insert a new paired node N N between N P P and N P , assign it with the color of N P , and apply the greedy algorithm on N N . Observe that N P always takes either black or white color changing the parity level of all its descendants (including N ). Note that the children of N P may get recolored, we can even get one more grey child but after this operation the parity levels of all children of N are correct and we do not change parity levels outside the subtree rooted at N . After fixing all nodes, we get a separated proper coloring (which is actually the greedy coloring) of T S . Hence, by Theorem 8, S is designable. Figure 7 illustrates this process. First, we will establish that R1, R2, R3 and R7 hold in the N model. R1 concerns an alphabet that does not allow base pairs to occur, so their weighting is unconsequential. R2 is equally trivial since the uniform weighting of every occurrence of a base pair type does not affect the relative order of structures.

R3 requires a clarification, as the introduction of two partners A and G for U somehow assign an unambiguous semantics to Σ 1,1 notation, leading to a disjunctive discussion. Let us assume that Σ 1,1 represents {A, C, G} (resp. {A, C, U}), then the argument used for R2 holds, since A (resp. C) cannot form base pairs. R7 is a direct corollary of Theorem 9 which is applicable for any bipartite energy model. Since the Nussinov-Jacobson model is bipartite, as shown in Figure 10, then R7 also holds in the N model. Definition 5 Let X ⊆ {C, G, A, U}. A design w for a structure S is X-unpaired if and only if the bases of w found at unpaired positions in S belong to X. If X = {b} is a singleton, the notation is shortened to b-unpaired.

Let n GC (w, S) (resp. n AU (w, S), n GU (w, S)) be the number of G-C (resp. A-U, G -U) base pairs of S on w. Note that

E W (w, S) = -n GC (w, S) -n AU (w, S) and E N (w, S) = α.n GC (w, S) + β.n AU (w, S) + γ.n GU (w, S). Proposition 11 Let ∆ W = 1 and 0 < ∆ N ≤ min(|α|, |β|), if a structure is A- unpaired and (W, Σ 2,0 , ∆ W )-designable then it is also (N , Σ 2,0 , ∆ N )-designable.
Proof Let w be an A-unpaired (W, Σ 2,0 , ∆ W )-design for S. Since in W there are no G -U base pairs, we have

n GC (w, S) = |w| G = |w| C , n AU (w, S) = |w| U ,
and for any other structure S ∈ S |w| , n GC (w, S ) ≤ n GC (w, S) n AU (w, S ) ≤ n AU (w, S) , with at least one of the inequalities being strict. We will show that w is also a (N , Σ 2,0 , ∆ N )-design for S. Consider any alternative structure S ∈ S |w| . If n GU (w, S ) = 0, then

E N (w, S ) = α n GC (w, S ) + β n AU (w, S ) ≥ αn GC (w, S) + β n AU (w, S) + min{|α|, |β|} ≥ E N (w, S) + ∆ N , as required. Otherwise (n GU (w, S ) > 0), first observe that n GC (w, S ) + n GU (w, S ) ≤ |w| G and n AU (w, S ) + n GU (w, S ) ≤ |w| U . (2)
Now, we have

E N (w, S ) = α n GC (w, S ) + β n AU (w, S ) + γ n GU (w, S ) ≥ (2) α |w| G -α n GU (w, S ) + β |w| U -β n GU (w, S ) + γ n GU (w, S ) = α |w| G + β |w| U + (γ -α -β)n GU (w, S ) ≥ E N (w, S) + γ -α -β.
Moreover, since α, β, γ < 0 and α, β < γ, then

γ -α -β ≥ max(|α|, |β|) ≥ min(|α|, |β|) ≥ ∆ N .
We conclude that E N (w, S ) ≥ E N (w, S) + ∆ N as required.

Corollary 12 Results R4, R5, R6 and R8 hold in any Nussinov-Jacobson energy model N with 0 < ∆ N ≤ min(|α|, |β|).

Proof The validity of Results R6 and R8 in N follows directly from a close inspection of the constructive proofs in the Watson-Crick energy model, both establishing the existence of A-unpaired designs. Proposition 11 therefore applies, and extends the validity of those designs to any suitable ∆ N . R5 follows from the fact that, in the proof of Theorem 7, our counterexamples 'locally' trade one base pair in S for another in the alternative structure, and that the two base pairs are of the same type. Therefore, both structures have the same energy in the Nussinov-Jacobson energy model.

From this, we conclude on the validity of R4. Indeed, any failure to the degree condition D(S) ≤ 4 implies the existence of m 5 in S, i.e. such structures are undesignable and the degree condition is therefore necessary.

Conclusion, discussion and perspectives

In this work, we introduced the Combinatorial RNA Design problem, a minimal instance of the RNA design problem which aims at finding a sequence that admits the target structure as its unique base pair maximizing structure using A -U and G -C base pairs. First, we provided complete characterizations for the structures that can be designed using restricted alphabets. Then we considered the RNA design under a four-letters alphabet, and provide a complete characterization of designable saturated structures, i.e., free of unpaired positions. Turning to those target structures that contain unpaired positions, we provided partial characterizations for classes of designable/undesignable structures, and showed that the set of designable structures is closed under the stutter operation. Finally, we introduced structure-approximating version of the problem and, assuming that the input structure avoids two motifs, provided a structure approximating algorithm of ratio 2 for general structures. We showed that our results also hold in the more realistic Nussinov-Jacobson energy model, which allows G -U base pairs to occur, and associates arbitrary negative free-energy to each base pair type (G -U being the weakest).

An important question that is left open by this work is the computational complexity of the RNA design problem. Schnall-Levin et al. [START_REF] Schnall-Levin | Inverting the Viterbi algorithm: an abstract framework for structure design[END_REF] established the NP-hardness of a more general problem, called the inverse Viterbi algorithm, which takes as input a stochastic grammar (representing the energy model) and a targeted parse tree (representing the structure), and outputs a sequence (design) whose most probable parsing should match the target. However this result does not settle the complexity of the RNA design, essentially because the proposed reduction relies critically on an encoding of 3-SAT instances within the input grammar. While the hypothetical perfect grammar/energy model for RNA folding probably differs from the currently accepted Turner model, it should ultimately reflect the laws of physics and should certainly not depend on the instance. As the reduction [START_REF] Schnall-Levin | Inverting the Viterbi algorithm: an abstract framework for structure design[END_REF] requires a different grammar (i.e., energy model) for each instance, it does not seem easily adaptable into a proof that holds for a fixed energy model. Consequently, despite two decades of work on the subject, the computational tractability of RNA design is still open, either in its general instance and in our purely combinatorial version.

In our opinion, this exceptional resistance of the RNA design problem to any attempt so far at characterizing its computational complexity can be attributed to two main reasons:

-The inverse nature of the problem: While polynomial, the direct computation MFE folding for an RNA requires dynamic programming, and runs in Θ(n 3 ) time (up to polylogarithmic factors, see [START_REF] Zakov | Reducing the worst case running times of a family of RNA and cfg problems, using valiant's approach[END_REF] for a complete state-of-the-art).

Unfortunately, the optimal-substructure property of the direct problem does not transpose to the inverse problem. Therefore, solving the RNA design problem somehow requires inverting a nontrivial -yet polynomially computable -function. It is tempting here to draw a parallel with some areas of cryptography, where multiple protocols are based on a -sometimes difficult to establish -disymmetry between the complexities of the direct and inverse computation.

-The intricacies of the objective function: Current state-of-the-art implementations of MFE folding prediction algorithms rely on a sophisticated energy model, the Turner model [START_REF] Turner | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF]. This model associates energy contributions to as much as 24 000 different types of structure/sequence motifs, and vastly increases the complexity the characterization of the space, and energies, of competing structures. On the other hand, oversimplified statements for the problem, as would result from a relaxation of the uniqueness condition, can be trivially solved in linear time. Such problems are not only largely unrealistic from a biological perspective, but they also probably do not retain the potential difficulty of the general problem.

Besides complexity issues, natural extensions of this work may include the consideration of more sophisticated energy models such as those based on stacking pairs or, ultimately, to the full Turner energy model [START_REF] Turner | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF]). One could also consider incorporating additional constraints, expressed as the presence/avoidance of motifs [START_REF] Zhou | Flexible RNA design under structure and sequence constraints using formal languages[END_REF], GC-content [START_REF] Reinharz | A weighted sampling algorithm for the design of RNA sequences with targeted secondary structure and nucleotide distribution[END_REF]. . . or the design under other objectives, such as the Boltzmann probability [START_REF] Zadeh | Nucleic acid sequence design via efficient ensemble defect optimization[END_REF]. In the Nussinov-Jacobson model, our result could be completed by the consideration of more liberal values for ∆, although it should be noted that considering larger such values would only gradually deplete the sets of designable structures until it becomes empty when the conditions of Theorem 1 are met. More precise bounds for the ratio of the structure-approximating could also be established. Finally, the structure-approximation problem could be revisited in an optimization setting, in which one would attempt to minimize the number of modifications made to the structure, so that a given structure becomes designable (or, more modestly, belongs to an identified class of designable structures). We plan to address some of these questions in future works.

Fig. 1

 1 Fig. 1 Four equivalent representations for an RNA secondary structure of length 68, consisting of 20 base pairs forming 7 bands: outer-planar graph (a.), arc-annotated representation (b.), parenthesized expression (c.), and tree representation (d)

  otherwise, where ε is the empty sequence. Define the S-paired restriction of w, denoted by Paired(w, S), as the subsequence of w consisting of the paired positions only, i.e. Paired(w, S) = u 1 • • • u |w| . Similarly, define the paired restriction of S = (n, P ), denoted by Paired(S), as the substructure of S consisting of paired positions only, i.e. Paired(S) = (|Paired(w, S)|, {(|u

Fig. 3

 3 Fig. 3 Examples of structures satisfying conditions of R1 (a.), R2 (b.) and R3 (c.).

  4 and S is saturated}.

Fig. 5

 5 Fig. 5 Examples of two structures containing motifs m 5 (a.) and m 3 • (b.), respectively. By R5 these two structures are not designable.

  S. Then Paired(w, S) is a design for the paired restriction Paired(S) of S. Since Paired(S) is saturated, Paired(w, S) is over alphabet Σ 1,0 ⊂ Σ 1,1 , and by Theorem 4, D(Paired(S)) ≤ 2. Hence, D(S) = D(Paired(S)) ≤ 2.

  R4 follows immediately from Theorem 4 by taking c = 2. Lemma 7 (Result R5) Any structure that contains m 5 or m 3 • is not Σ 2,0designable. Proof Assume that S is Σ 2,0 -designable and let w be a design for S. Then Paired(w, S) is a design for Paired(S). Since Paired(S) is saturated, by Theorem 4, D(S) = D(Paired(S)) ≤ 4, hence, S cannot contain motif m 5 . Now, assume that S contain motif m 3 • appearing at node [a, b] of T S . Let {[l i , r i ]} 3 i=1 be some paired children of [a, b] and the node [a, b] if [a, b] is not the root, and [u, u] an unpaired child of [a, b].

Fig. 9

 9 Fig. 9 Exhaustive search and systematic counter-examples for the design of m 3 • . A -B and C -D respectively represent the first and second pairs of letters found in the design in prefix order, allowing the factorization of trivial symmetries.

4. 2

 2 Nussinov-Jacobson energy model N We will show that, for a value of ∆ ∈ (0, min(|α|, |β|)] our results for the Watson-Crick model transpose to the more general Nussinov-Jacobson model.

Acknowledgements The authors thank Cédric Chauve for fruitful discussions and constructive criticisms. YP is supported by the Pacific Institute for the Mathematical Sciences (PIMS), and the French Agence Nationale de la Recherche (ANR-14-CE34-0011).

is a design, a contradiction. Let w i be the leftmost sequence such that w i 1 is not paired with w i |w i | in S . Since S must be also saturated, w i 1 must be paired. Let w j k , j ≥ i, be the partner of w i 1 in S , and let u :

, then j > i and, by complementarity, w i 1 = w j 1 which contradicts the preconditions. Hence, we can assume that k < |w j |. Since u and each of the w i , . . . , w j-1 are saturable, by iterated application of Claim 4.1, we conclude that v = w j [1,k] is saturable as well. This contradicts the precondition that w j is an atomic saturable design, since v is a proper prefix of w j . We conclude that no alternative saturated folding exists for W , i.e., W is a design for S.

Claim 4.3 Consider t atomic saturable designs w

Proof We will first show that W is an atomic saturable sequence. Assume to the contrary that there is a proper prefix of W that is saturable. Consider the shortest such prefix aw 1 • • • w i w i+1

[1,j] with 1 ≤ j < |w i+1 | and 1 ≤ i < t. Obviously, a has to be paired with w i+1 j , otherwise we can find a shorter saturable prefix. This implies that b = w i+1 j and that w 1 • • • w i w i+1

[1,j-1] is saturable as well. By repeated application of Claim 4.1, we have that w i+1

[1,j-1] is saturable. Since it is a prefix of atomic saturable sequence w i+1 , it must be the empty sequence, i.e., j = 1. Therefore, b = w i+1 1 , a contradiction with the assumptions of the claim. Thus, W is an atomic saturable sequence. Now we will show that W is a design. Consider any MFE (saturated) structure S for W . Since W is atomic saturable, a is paired with b in S. By Claim 4.2, w 1 • • • w t is a design. It follows that W is a design as well.

To prove the sufficiency of the degree condition, consider the following algorithm, which takes as input a saturated structure S with D(S) ≤ 2c, and returns a design w for S:

-Let {[l i , r i ]} d i=1 be the children of the root. Assign to each w l i , w r i complementary bases such that ∀1 ≤ i < j ≤ d : w l i = w l j .