Patrick Foulon
email: foulon2@cirm-math.fr

Vladimir S Matveev
email: vladimir.matveev@uni-jena.de

Dmitri Burago

ZERMELO DEFORMATION OF FINSLER METRICS BY KILLING VECTOR FIELDS

Keywords: Mathematics Subject Classification. Primary: 53B40. Secondary: 53C35 Finsler, Zermelo deformation, Killing vector fields, Jacobi fields, curvature, symmetric spaces

We show how geodesics, Jacobi vector fields, and flag curvature of a Finsler metric behave under Zermelo deformation with respect to a Killing vector field. We also show that Zermelo deformation with respect to a Killing vector field of a locally symmetric Finsler metric is also locally symmetric.

Introduction

Let F be a Finsler metric on M n and v be a vector field such that F (x, -v(x)) < 1 for any x ∈ M n . We will denote by F the Zermelo deformation of F by v. That is, for each point x ∈ M , the unit F -ball Bx := {ξ ∈ T x M n | F (x, ξ) < 1} is the translation in T x M n along the vector v(x) of the unit F -ball B x := {ξ ∈ T x M n | F (x, ξ) < 1} (see Figure 1).

v ξ J ξ + v J Figure 1.
The unit ball of F (dashed line) is the v-translation of that of F (bold line). If a vector J is tangent to the unit ball of F at ξ, it is tangent to the unit ball of F at ξ + v

Equivalently, this can be reformulated as F (x, ξ) = F (x, ξ -F (x, ξ)v(x)).

Indeed, Equation ( 1) is positively homogeneous, and for any ξ such that F (x, ξ) = 1 we have F (x, ξ -v(x)) = 1. The first result of this note is a description of how geodesics, Jacobi vector fields, and flag curvatures of F and of F are related if the vector field v is a Killing vector field for F , that is, if the flow of v preserves F .

Theorem 1. Let F be a Finsler metric on M n admitting a Killing vector field v such that F (x, -v(x)) < 1 for all x ∈ M n . We denote by Ψ t the flow of v and by F the v-Zermelo deformation of F . Then, for any F -arc length parametrized geodesic γ of F , the curve t → Ψ t (γ(t)) which we denote by γ(t) is an F -arc length parametrized geodesic of F .

Moreover, for any Jacobi vector field J(t) along γ such that it is orthogonal to

γ(t) in the metric g (γ(t), γ(t)) := 1 2 d 2 ξ F 2 (γ(t), γ(t))
, the pushforward J(t) = Ψ t * (J(t)) is a Jacobi vector field for γ(t) and is orthogonal to

γ(t) in g(γ(t), γ(t)) := 1 2 d 2 ξ F 2 (γ(t), γ(t))
. Moreover, flag curvatures K and K of F and F are related by the following formula: for any x ∈ M and any "flag" (ξ, η) with "flagpole ξ ∈ T x M with F (x, ξ) = 1" and transverse edge η ∈ T x M , we have K(x, ξ, η) = K(x, ξ + v, η) provided that ξ + v and η are linearly independent.

We do not pretend that the whole result is new but rather suggest that certain parts of it are known. The first statement of Theorem 1 appears in [START_REF] Katok | Ergodic properties of degenerate integrable Hamiltonian systems[END_REF]. We recall the arguments of A. Katok in Remark 1. The third statement was announced in [START_REF] Foulon | Ziller-Katok deformations of Finsler metrics[END_REF] and follows from the recent paper [START_REF] Huang | On the flag curvature of a class of Finsler metrics produced by the navigation problem[END_REF]. Special cases when the metric F is Riemannian were studied in detail in, e.g., [START_REF] Bao | Zermelo Navigation on Riemannian manifolds[END_REF][START_REF] Shen | Finsler manifolds of constant positive curvature[END_REF]. Though we did not find the second statement of Theorem 1, the one about the Jacobi vector fields, in the literature, we think it is known in folklore.

Unfortunately, in all these references, the proof is by direct calculations, which are sometimes quite tricky and sometimes require a lot of preliminary work. One of the goals of this note is to show the geometry lying below Theorem 1 and to demonstrate that certain parts of Theorem 1 at least are almost trivial.

Our second result shows that Zermelo deformation with respect to a Killing vector field preserves the property of a Finsler metric to be a locally symmetric space. We will call a Finsler metric locally symmetric, if for any geodesic γ the covariant derivative of the Riemann curvature (= Jacobi operator) vanished:

D γ R γ = 0. (2)
Here D γ stays for the covariant derivative along the geodesic: D γ = ∇ γ γ . Both most popular Finslerian connections, Berwald and Chern-Rund connections, can be used as the Finslerian connection in the last formula, see, e.g., [10, §7.3], whose notation we partially follow.

Remark 1. In Riemannian geometry there exist two equivalent definitions of locally symmetric spaces: according to the "metric" definition, a space is locally symmetric if for any point there exists a local isometry such that this point is a fixed point and the differential of the isometry at this point is minus the identity. By the other "curvature" definition, a space is locally symmetric if the covariant derivative of the curvature tensor is zero. The equivalence of these two definitions is a classical result of E. Cartan. We see that our definition above is the generalization to the Finsler metrics of the "curvature" definition; it was first suggested in [START_REF] Foulon | Locally symmetric Finsler spaces in negative curvature[END_REF].

In the Finsler setup, both "metric" and "curvature" definitions are used in the literature, but they are not equivalent anymore: the symmetric spaces with respect to the "metric" definition are symmetric spaces with respect to the "curvature" definition, but not vice versa.

In fact, the "metric" definition is much more restrictive; in particular locally symmetric metrics in the "metric" definition are automatically Berwaldian [8, Theorem 9.2] and are clearly reversible. On the other hand, all metrics of constant flag curvature, in particular all Hilbert metrics in a strictly convex domain, are locally symmetric in the "curvature" definition, and are symmetric with respect to the "metric" definition if and only if the domain is an ellipsoid.

In view of this, the name "locally symmetric" is slightly misleading, since locally symmetric manifolds may have no (local) isometries. We will still use this terminology because it was used in the literature before.

Theorem 2. Suppose that F is a locally symmetric Finsler metric and v a Killing vector field satisfying F (x, -v(x)) < 1 for all x. Then, the v-Zermelo deformation of F is also locally symmetric.

All Finsler metrics in our paper are assumed to be smooth and strictly convex but may be irreversible.

Proofs.

2.1. Proof of Theorem 1. Let γ(t) be an arc length parametrized F -geodesic, we need to prove that the curve t → Ψ t (γ(t)) is an arc length parametrized Fgeodesic. In order to do it, observe that for any F -arc length parametrized curve x(t) the t-derivative of Ψ t (x(t)) is given by Ψ t * ( ẋ(t)) + v(Ψ t (x(t))). Since the flow Ψ t preserves F and v, it preserves F and therefore

F (Ψ t (x(t)), Ψ t * ( ẋ(t)) + v(Ψ t (x(t)))) = F (x(t), ẋ(t) + v(x(t))) = F (x(t), ẋ(t)) .
The last equality in the formula above is true because, for any ξ such that F (x, ξ) = 1, we have F (x, ξ + v(x)) = F (x, ξ) by the definition of the Zermelo deformation. Thus, if the curve x(t) is F -arc length parametrized, then the curve Ψ t (x(t)) is F -arc length parametrized.

This also implies the integrals F (x(t), ẋ(t))dt and F Ψ t (x(t)), (Ψ t (x(t))) dt coincide for all F -arc length parametrized curves x(t). Since geodesics are locally the shortest arc length parametrized curves connecting two points, for each arc length parametrized F -geodesic γ the curve t → Ψ t (γ(t)) is an F -arc length parametrized geodesic as we claimed.

Remark 2. Alternative geometric proof of the statement that for each arc length parametrized F -geodesic γ the curve t → Ψ t (γ(t)) is an F -arc length parametrized geodesic is essentially due to [START_REF] Katok | Ergodic properties of degenerate integrable Hamiltonian systems[END_REF]: consider the Legendre-transformation T : T * M n → T M corresponding to the function 1 2 F 2 and denote by F * the pullback of F to T * M n , F * := F • T . Next, view the vector field v as a function on T * M by the obvious rule η ∈ T *

x M n → η(v(x)). It is known that the Hamiltonian flow corresponding to the function v is the natural lift of the flow of the vector field v to T * M . Since v is assumed to be a Killing vector field, the Hamiltonian flows of F * and of v commute. Next, consider the function

F * := F * + v. If v satisfies F (x, -v(x)) < 1,
then the restriction of F * to T x M n is convex, consider the Legendre-transformation T : T M n → T * M corresponding to the function 1 2 ( F * ) 2 and the pullback of F * to T M n , it is a Finsler metric which we denote by F . It is a standard fact in convex geometry that the Finsler metric F is the v-Zermelo-deformation of F . Since the Hamiltonian flows of F * and of v, which we denote by ψ t and d * Ψ t , commute, the Hamiltonian flow of F * is simply given by

ψt = d * Ψ t • ψ t . (3)
Then, for any point (x, ξ) ∈ T M with F (x, ξ) = 1, the projections of the orbits of ψt and of ψ t starting at this point are arc length parametrized geodesics γ of F and γ of F . By (3) we have γ(t) = Ψ t • γ(t) as we claimed.

Let us now prove the second statement of Theorem 1. Consider a Jacobi vector field J(t) which is orthogonal to γ. We need to show that the pushforward J(t) = Ψ t * (J(t)) is a Jacobi vector field for the F -geodesic γ(t) := Ψ t (γ(t)). By the definition of Jacobi vector field there exists a family γ s (t) of geodesics with γ 0 = γ such that J(t) = d ds |s=0 γ s (t), since J(t) is orthogonal to γ we may assume that all geodesics γ s (t) are arc length parametrized. As we explained above, Ψ t (γ s (t)) is a family of F -geodesics; taking the derivative by s at s = 0 proves what we want.

Let us now show that J is orthogonal to γ. First observe that the condition that J(t) is orthogonal to γ(t) is equivalent to the condition that J(F ) := r J r ∂F ∂ξr vanishes at γ(t) for each t. Indeed, consider the one-form U ∈ T γ(t) M n → g γ(t), γ(t) ( γ(t) , U ). Because of the (positive) homogeneity of the function F we have that at a point γ(t)

∈ T γ(t) M n g (γ(t), γ(t)) ( γ(t) , U ) = U r ∂F ∂ξ r . (4)
Next, take Equation (1) and calculate the differential of the restriction of F to the tangent space: its components are given by ∂ F

∂ξ i = 1 1 + v(F ) ∂F ∂ξ i . (5)
In this formula, the derivatives of the function F are taken at ξ ∈ T x S 2 , and the derivatives of the function F are taken at ξ -F (x, ξ)v. By v(F ) we denoted the function r ∂F ∂ξr v r . In view of (5), J( F ) := r J r ∂ F ∂ξr vanishes at γ(t)+v(γ(t)), so J(t) is orthogonal to γ(t) + v(γ(t)) (the orthogonality is understood in the sense of g(γ(t), γ(t)+v(γ(t))) ). Then, J(t) = Ψ t * (J(t)) is g(γ(t), γ(t)) orthogonal to γ(t). Remark 3. Geometrically, the just proved statement that J is orthogonal to γ, after the identification of T γ(t) M n and T Ψt(γ(t)) M n by the differential of the diffeomorphism Ψ t , corresponds to the following simple observation: if J is tangent to the unit F -sphere at the point ξ = γ, then it is also tangent to the unit F -sphere at the point ξ + v, see Figure 1.

Let us now prove the third statement of Theorem 1, we need to show that K(x, ξ, η) = K(x, ξ + v, η). We consider an F -geodesic γ(t) with γ(0) = 0 and γ(0) = ξ and the corresponding F -geodesic γ(t) = Ψ t (γ(t)); since Ψ 0 = Id, we have γ(0) = γ(0) + v = ξ + v. Observe that by combining [2, Eqn. 6.16 on page 117] and [2, Eqn. 6.3 on page 108] we obtain 1 2

d 2 dt 2 g(J(t), J(t)) = g(D γ(t) D γ(t) J(t), J(t)) + g(D γ(t) J(t), D γ(t) J(t)) = -K(γ(t), γ(t), J(t))(g( γ(t), γ(t))g(J(t), J(t)) (6)
-g( γ(t), J(t)) 2 ) + g(D γ J, D γ J).

We will assume that J(0) := η is g-orthogonal to γ(0). As explained above this implies that γ(0) is g-orthogonal to J(0). Then, by [START_REF] Javaloyes | Zermelo navigation in pseudo-Finsler metrics[END_REF], the minimum of

d 2
dt 2 g(J(t), J(t)) |t=0 taken over all Jacobi vector fields J along γ which are equal to η at t = 0, is equal to -K(x, ξ, η)g(ξ, ξ)g(η, η). An analogous statement is clearly true also for γ, g and J: namely, the minimum of d 2 dt 2 g( J(t), J(t)) |t=0 taken over all Jacobi vector fields J along γ which are equal to η at t = 0, is equal to -K(x, ξ + v, η)g(ξ + v, ξ + v)g(η, η). Here we used the relation J(0) = Ψ t * (J(t)) |t=0 = J(0). Finally, in order to show that K(x, ξ, η) = K(x, ξ + v, η), it is sufficient to show that the function t → g(J(t), J(t)) is proportional, with a constant coefficient, to the function t → g( J(t), J(t)).

In order to prove this, let us first compare

g (γ(t), γ(t)) = 1 2 d 2 ξ F 2 (γ(t), γ(t)) and g(γ(t), γ(t)) = 1 2 d 2 ξ F 2 (Ψt•γ(t),Ψt * ( γ(t))+v(Ψt(γ(t)))) .
It is convenient to work in coordinates (x 1 , . . . , x n ) such that the entries of v are constants, in these coordinates for each t the differential of the diffeomorphism Ψ t is given by the identity matrix, so in these coordinates J(t) = J(t) and γ(t) = γ(t) + v(γ(t)). Differentiating (5), we get the second derivatives of F . They are given by

∂ 2 F ∂ξ i ∂ξ j = 1 1 + v(F ) ∂ 2 F ∂ξ i ∂ξ j - 1 (1 + v(F )) 2 r ∂ 2 F ∂ξ i ∂ξ r v r ∂F ∂ξ j + ∂ 2 F ∂ξ j ∂ξ r v r ∂F ∂ξ i . (7)
Again, all derivatives of the function F are taken at ξ, and of the function F are taken at ξ -F (x, ξ)v(x). Note that one term in the brackets in (7) appears because we differentiate 1 1+v(F ) , and the other appears because the derivatives of ∂F ∂ξj are taken at ξ -F (x, ξ)v(x). When we differentiate it, we also need to take into account the additional term -F (x, ξ)v(x). Now, in view of the formula

1 2 d 2 ( F 2 ) = F d 2 F + d F ⊗ d F we obtain from (7) the formula for gij : gij = F 1 + v(F ) g ij - F (1 + v(F )) 2 r ∂ 2 F ∂ξ i ∂ξ r v r ∂F ∂ξ j + ∂ 2 F ∂ξ j ∂ξ r v r ∂F ∂ξ i + 1 (1 + v(F )) 2 ∂F ∂ξ i ∂F ∂ξ j . (8)
Let us now compare the length of J in g(x, ξ) with that of in g(x, ξ + v). We multiply (8) by J i J j and sum with respect to i and j. Since by assumptions J(F ) = r J r ∂F ∂ξr vanishes at ξ, all terms in the sum but the first vanish. We thus obtain that the length of J in g is proportional to that of in g with the coefficient which is the square root of F 1+v(F ) . But along the geodesic both F and v(F ) are constant. Indeed, v(F ) is the "Noether" integral corresponding to the Killing vector field. Theorem 1 is proved.

Remark 4. After this paper was accepted we learned that our Theorem 1 also follows from [6, Theorems 1.2 and 1.3], which at the time of writing is in preprint.

2.2.

Proof of Theorem 2. First, observe that a Finsler metric is locally symmetric if and only if for any geodesic γ and any Jacobi vector field J along γ the vector field D γ J is also a Jacobi vector field. Indeed, the equation for Jacobi vector fields is

D γ D γ J + R γ (J) = 0. (9)
D γ -differentiating this equation, we obtain

D γ (D γ D γ J + R γ (J)) = D γ D γ (D γ J) + R γ (D γ J) + (D γ R γ ) (J) = 0.
If D γ J is a Jacobi vector field, D γ D γ (D γ J) + R γ (D γ J) vanishes so the equation above implies (D γ R γ ) (J) = 0, and since it is fulfilled for all Jacobi vector fields we have D γ R γ = 0 as we claimed. Thus, we assume that for any geodesic and for any Jacobi vector field for F its D γ derivative is also a Jacobi vector field, and our goal is to show the same for F . Clearly, it is sufficient to show this only for Jacobi vector fields which are g-orthogonal to γ. Note that for such Jacobi vector fields D γ J is also orthogonal to γ, since both g (γ, γ) and γ are D γ -parallel.

Take a (arc length parametrized) F -geodesic γ and a point P = γ(0) on it. Consider the geodesic polar coordinated around this point, let us recall what they are and their properties which we use in the proof.

Consider the (local) diffeomorphism of T P M n \ {0} to M n which sends ξ ∈ T P M n \ {0} to exp(ξ) := γ ξ (1), where γ ξ is the geodesic starting from P with the velocity vector ξ. As the local coordinate systems on T P M n \ {0} we take the following one: we choose a local coordinate system x 1 , . . . , x n-1 on the unit F -sphere {ξ ∈ T P M n | F (P, ξ) = 1} and set the tuple F (P, ξ), x 1 1 F ξ , . . . , x n-1 1 F ξ to be the coordinates of ξ. Combining it with the diffeomorphism exp, we obtain a local coordinate system on M n . By construction, in this coordinate system each arc length parametrized geodesics starting at P , in particular the geodesic γ, is a curve of the form (t, const 1 , . . . , const n-1 ). Next, consider the following local Riemannian metric ĝ in a punctured neighborhood of P : for a point σ(t) of this neighborhood such that σ is a geodesic passing through P we set ĝ := g (σ(t), σ(t)) . It is known that in the polar coordinates the metric ĝ is block-diagonal with one 1 × 1 block which is simply the identity and one (n -1) × (n -1)-block which we denote by G:

ĝ = 1 G .
It is known, see, e.g., [START_REF] Shen | Differential Geometry of Spray and Finsler Spaces[END_REF]Lemma 7.1.4], that the geodesics passing through P are also geodesics of ĝ, and that for each such geodesic the operator ∇ γ , where ∇ is the Levi-Civita connection of ĝ, coincides with D γ . Next, consider analogous objects for the metric F . As the local coordinate system on the unit F -sphere we take the following: as the coordinate tuple of ξ with F (P, ξ) = 1 we take the coordinate tuple (x 1 (ξ), . . . , x n-1 (ξ)), where ξ := ξ -v(P ). (Recall that the v-parallel-transport sends the unit F -sphere to the unit F -sphere.) By Theorem 1, in these coordinate systems each Jacobi vector field J = (J 0 (t), . . . , J n-1 (t)) along γ which is orthogonal to γ is also a Jacobi vector field along γ, which is the F -geodesic such that γ(0) = P and γ(0) = γ(0) + v(P ), and is orthogonal to γ. By [START_REF] Matveev | The Binet-Legendre metric in Finsler geometry[END_REF], the corresponding block G is given by G = 1 1+v(F ) G. Since the function v(F ) is constant along geodesics, the coefficients Γ k ij of the Levi-Civita connection ∇ of ĝ such that i = 0 or j = 0 coincide with that of for the analog for F . A direct way to see the last claim is to use the formula Γ k ij = 1 2 g ks ∂gis ∂xj + ∂gjs ∂xi -∂gij ∂xs , where all indices run from 0 to n -1 and the summation convention is assumed. Then, in our chosen coordinate system, the formula for the covariant derivative in ∇ along γ for vector fields which are orthogonal to γ simply coincides with that of the formula for the corresponding objects for F . Then, for any F -Jacobi vector field J orthogonal to γ we have that D γ J is again a Jacobi vector field. Theorem 2 is proved.

The authors thank Sergei Ivanov for useful comments. V. M. was partially supported by the University of Jena and by the DFG grant MA 2565/4.


Sergio Peignier

Pauline Schmitt

Federica Calevro

Data-driven Gene Regulatory Network Inference based on Classification Algorithms

Keywords: Bioinformatics, Gene Regulatory Network Inference, Classification

Different paradigms of gene regulatory network inference have been proposed so far in the literature. The data-driven family is an important inference paradigm, that aims at scoring potential regulatory links between transcription factors and target genes, analyzing gene expression datasets. Three major approaches have been proposed to score such links relying on correlation measures, mutual information metrics, and regression algorithms. In this paper we present a new family of data-driven inference approaches, inspired on the regression based family, and based on classification algorithms. This paper advocates for the use of this paradigm as a new promising approach to infer gene regulatory networks. Indeed, the implementation and test of five new inference methods based on well-known classification algorithms shows that such an approach exhibits good quality results when compared to well-established paradigms.

Introduction

Gene Regulatory Networks (GRNs) describe the complex interactions between specialized genes such as transcription factors (TFs) and their target genes (TGs). Such interactions mediate to a large extent the regulation of the expression of genes, and the adaptation of biological systems to different conditions. Understanding such networks, and analysing their organization and their dynamics, are therefore important steps towards the comprehension of complex mechanisms that shape living organisms. The development of high-throughput technologies has motivated the creation of computational methods to reverse-engineer the underlying GRNs. Inferring GRNs from high-throughput data is a challenging problem, and different methods tackling this have been proposed so far in the system biology literature [START_REF] Sanguinetti | Gene regulatory network inference: an introductory survey[END_REF]. Among the major categories of methods, the so-called data-driven approaches are among the most popular techniques due to their simplicity, their computational efficiency and their accuracy [START_REF] Sanguinetti | Gene regulatory network inference: an introductory survey[END_REF]. These methods aim at scoring each possible regulatory link, by estimating the dependency between genes from experimental high-throughput data. Depending on the method used to infer the link between genes, data-driven methods have been classified in three major families [START_REF] Sanguinetti | Gene regulatory network inference: an introductory survey[END_REF], gathering respectively methods based on i) correlation metrics, ii) mutual information, and iii) feature importance scoring based on regression algorithms. In practice, scoring the dependency between TFs and TGs using feature importance, could similarly be achieved using classification algorithms, however this paradigm has been understudied in the literature.

In order to assess the effectiveness of this new paradigm, we present in this paper a computational framework including 5 data-driven inference methods, based on classification algorithms. The proposed methods have been assessed and compared to state-of-the-art approaches using the benchmark datasets and the evaluation procedure described in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF]. Our comparative analysis revealed that the proposed methods allow to obtain better results than state-of-the-art techniques, and thus, methods based on this paradigm are interesting tools for the analysts. We also included 8 well-known preprocessing techniques that were also tested to study their impact on the inference quality. For the sake of reproducibility, our framework as well as the result tables are available online 1 .

The rest of the paper is organized as follows. Section 2 presents the related work. Section 3 describes more formally the GRN inference problem. Section 4 introduces major preprocessing techniques that were tested. Section 5 introduces our data-driven GRN inference paradigm, based on classification algorithms. Section 6 and 7 describe respectively the experimental assessment protocol and the results. We conclude with a summary and some perspectives in Section 8.

2 State-of-the-art

Gene regulatory network inference methods

Gene regulatory network inference techniques from gene expression matrix can be classified in three major families, namely Model-Based, Data-Driven and Multi-Network approaches. This section reviews the founding principles of these families, as well as previous works on GRN inference based on classification algorithms.

Model-Based

These approaches rely on a predetermined GRN model, based on specific hypothesis and parameters. In this context, inferring a GRN, consists in fitting the parameters of the model with respect to experimental data. Once fitted, the GRN model can be used to simulate and study the biological system in-silico. Two major families of model-based approaches have been reported in [START_REF] Sanguinetti | Gene regulatory network inference: an introductory survey[END_REF], namely Probabilistic Models and Dynamical Models. The former family mainly includes Bayesian networks and Gaussian Graphical Models. While the latter one has been used to study time series of gene expression data, incorporating modeling techniques that can integrate temporal considerations, such as Dynamic Bayesian Networks, Ordinal Differential Equations, Boolean Networks, Probabilistic Boolean Networks and Neural Networks.

We refer the reader to [START_REF] Chai | A review on the computational approaches for gene regulatory network construction[END_REF] for a thorough survey of model-based GRN inference and modeling methods.

Data-Driven methods These methods aim at scoring each possible regulatory link between TFs and TGs, by estimating their level of dependency, using experimental data. According to [START_REF] Sanguinetti | Gene regulatory network inference: an introductory survey[END_REF] data-driven methods are very popular techniques due to their simplicity, their speed and their accuracy. These algorithms have been classified in three major families in [START_REF] Sanguinetti | Gene regulatory network inference: an introductory survey[END_REF], depending on the method used to estimate the dependency between genes. 1) The first family of methods is based on the assumption that a TG and a TF regulating its expression should exhibit correlated gene expressions, and thus rely on correlation scores to infer regulatory links (e.g., [START_REF] Zhang | A general framework for weighted gene co-expression network analysis[END_REF]). 2) The second family of methods relies on information theory scores such as Mutual Information, to capture more complex relationships between TFs and TGs, which cannot be apprehended by linear correlations (e.g., [START_REF] Faith | Large-scale mapping and validation of escherichia coli transcriptional regulation from a compendium of expression profiles[END_REF]). 3) The last family reported in [START_REF] Sanguinetti | Gene regulatory network inference: an introductory survey[END_REF] is based on feature importance scores assigned by regression algorithms that are trained to predict the expressions of a TG from those of TFs (e.g., [START_REF] Irrthum | Inferring regulatory networks from expression data using tree-based methods[END_REF][START_REF] Haury | Tigress: trustful inference of gene regulation using stability selection[END_REF]).

Multi-Network methods Unlike the previous families, the methods belonging to this family tackle the GRN inference problem by combining heterogeneous data sources such as gene expression data, TF binding site motifs, or Chromatin Immuno-Precipitation data. For instance a recent method called SCENIC [START_REF] Aibar | Scenic: single-cell regulatory network inference and clustering[END_REF], refines the output of the GE-NIE3 data-driven method [START_REF] Irrthum | Inferring regulatory networks from expression data using tree-based methods[END_REF], using cis-regulatory TF binding site motif analysis. Another recent method, called PANDA [START_REF] Glass | Passing messages between biological networks to refine predicted interactions[END_REF][START_REF] Glass | High performance computing of gene regulatory networks using a message-passing model[END_REF] reconstructs GRNs by integrating heterogeneous sources of data, using a message passing approach.

Classification applied to GRN inference So far, some studies have used classification algorithms to infer GRNs in a supervised way. For instance, Support Vector Machine (SVM) classifiers (e.g., [START_REF] Mordelet | Sirene: supervised inference of regulatory networks[END_REF][START_REF] Gillani | Comparesvm: supervised, support vector machine (svm) inference of gene regularity networks[END_REF]), have been used to reverse-engineer GRNs. In these previous works, the GRN inference problem is seen as a binary classification task: given a TG and a TF, the task consists in classifying whether their interaction is true or not. Such methods depend to a large extent on their training datasets, and even if an increasing number of datasets describing regulatory interactions is available, unsupervised approaches remain necessary [START_REF] Cerulo | Learning gene regulatory networks from only positive and unlabeled data[END_REF]. The paradigm presented in this paper differs from these previous approaches, in the sense that it is not a supervised approach requiring training datasets. Indeed, in our framework classification algorithms are simply used to score the dependency between TFs and TGs from gene expression data, using an approach analogous to the regression-based one.

Problem Statement

Gene expression dataset Let X ∈ R I×J be a gene expression matrix (e.g., derived from RNAseq or Microarray experiments), such that Xi,j represents the level of expression of gene i in condition j. Moreover let Xi,. denote the vector of the gene i expression level across all conditions, and let X.,j be the vector of the expression level of all genes for condition j. Let J and I denote respectively the number of conditions (i.e., columns) and the number of genes (i.e., rows), considered in X. In most of the cases, tens of thousands of genes (I) are described in at most a few hundred conditions (J), and thus often J I.

Gene regulatory network GRNs are traditionally represented as an oriented graph, which nodes represent genes and edges denote regulatory links between TFs and their TGs. More formally, let T G = {tg 1 , . . . , tg I } denote the set of I genes of a given organism, and let T F ⊂ T G be the subset of genes corresponding to TFs. G = T G, E denotes a GRN. The set of nodes of G is simply the set of genes T G. And the set of oriented edges E represent the set of regulatory links, such that (tf , tg) ∈ E indicates that the transcription factor tf ∈ T F regulates the gene expression of target gene tg ∈ T G. Notice that a TF can also be the TG of another TF.

Data-Driven GRN inference Methods belonging to this family aim at using gene expression data to score all possible regulatory links, and then selecting the most promising edges only. Let E f ull = {(tf , tg) ∈ T F × T G | tf = tg} be the set of all possible links between TFs and TGs (excluding self-loops). Let us consider a function w : R I×J , T F , T G → R such that w(X, tf , tg) is the dependency score associated to the regulatory link (tf , tg) ∈ E f ull , inferred from the gene expression dataset X.

Then, a subset of regulatory links is usually chosen to define a putative GRN, by selecting for instance the k links with the highest scores. Therefore, the output of a data-driven inference algorithm strongly depends on the scoring function it relies on. Data-Driven GRN inference is a particularly complex task, since it suffers from the problem known as the "curse of dimensionality" [START_REF] Banf | Computational inference of gene regulatory networks: approaches, limitations and opportunities[END_REF]. Indeed, the number of genes I and, consequently, the number of possible interactions between them, i.e., |E |, is often orders of magnitude larger than the number of experimental conditions J, which leads to a large number of possible GRNs that can explain the experimental gene expression data [START_REF] Banf | Computational inference of gene regulatory networks: approaches, limitations and opportunities[END_REF].

Preprocessing techniques

Standardization

As in many machine learning problems, an important preliminary step consists in standardizing the dataset. For instance in [START_REF] Cheadle | Analysis of microarray data using z score transformation[END_REF], the authors have shown the importance that standardization techniques, such as the well-known z-score transformation, may have on gene expression data analysis. In this work we tested three standardization methods based on z-scores, described more formally hereafter.

• Z-score rows replaces each entry of the gene expression matrix as follows Xi,j ←

X i,j -µ i σ i
. Where µi = j Xi,j/J represents the average gene expression of gene i and σi = j (Xi,j -µi) 2 /(J -1) denotes its standard deviation. This standardization ensures that all genes have comparable levels of expression.

• Z-score columns replaces each entry of the gene expression matrix as follows:

Xi,j ← X i,j -µ j σ j

. Where µj = i Xi,j/I denotes the average gene expression of genes in condition j and σj = i (Xi,j -µj) 2 /(I -1) represents its standard deviation. This standardization technique ensures that all conditions have comparable ranges of expression.

• Polishing standardization [START_REF] Olshen | Successive normalization of rectangular arrays[END_REF] applies iteratively the z-score standardization along columns and rows until convergence. In practice, a few iterations are usually sufficient to converge. Here we only ran 10 iterations, which revealed to be sufficient to reach convergence.

These standardization techniques have been implemented using Pandas Python 3.7 library [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF].

Discretization

The discretization of gene expressions is another common preprocessing step for several GRN inference algorithms [START_REF] Li | Comparative study of discretization methods of microarray data for inferring transcriptional regulatory networks[END_REF]. This procedure aims at partitioning continuous variables to discretized intervals. In this work we included five well-known discretization techniques for gene expression, studied in [START_REF] Jung | Evaluation of data discretization methods to derive platform independent isoform expression signatures for multi-class tumor subtyping[END_REF], and described formally hereafter.

• EFD (Equal Frequency Discretization) partitions, for each gene, the full range of observed gene expression values in K subsets of equal sizes using K-quantiles as boundaries. Specifically, let Xi,. be the expression of gene i, EFD computes the K-quantiles of the Xi,..

Let Q k denotes the k-th K-quantile of Xi,. such that P r[x < Q k ] ≤ k/K, for x ∈ Xi,.. Then ∀j ∈ {1, . . . , J} Xi,j ← k if Q k-1 ≤ Xi,j < Q k for k ∈ {1, . . . , K}.
• EWD (Equal Width Discretization) splits, for each gene, the full range of observed gene expression values in K bins with equal size. More formally let Xi,. be the vector of expressions of gene i, and let min(Xi,.) and max(Xi,.) denote respectively the minimal and the maximal gene expression of gene i. EWD computes K equal width bins δ = (max(Xi,.)-min(Xi,.)

)/K. Then ∀j ∈ {1, . . . , J} Xi,j ← k if δ × (k -1) ≤ Xi,j < δ × k for k ∈ {1, . . . , K}.
• Rowkmeans (Row Kmeans Discretization) applies the K-means clustering algorithm [START_REF] Lloyd | Least squares quantization in pcm[END_REF] to partition the gene expression values of gene i into K clusters, and then the values are discretized according to their cluster membership. More formally let Xi,j be the expression of gene i in condition j, ∀j ∈ {1, . . . , J}, the values Xi,j are partitioned in K clusters. Let C k denote the k-th cluster, and let µ k be its centroid, we consider that clusters are sorted according to their centroid location, i.e., µ1 < µ2, < • • • < µK . Then ∀j ∈ {1, . . . , J} Xi,j ← k if Xi,j ∈ C k for k ∈ {1, . . . , K}

• Cokmeans (Column Kmeans Discretization) applies the Kmeans clustering [START_REF] Lloyd | Least squares quantization in pcm[END_REF] to partition the gene expression values of condition j into K clusters, and then the values are discretized according to their cluster membership, similarly as for Rowkmeans.

• Bikmeans (Bidirectional Kmeans Discretization) [START_REF] Li | Comparative study of discretization methods of microarray data for inferring transcriptional regulatory networks[END_REF] computes the Cokmeans and the Rowkmeans discretizations for a given gene expression matrix, combining them in order to consider both genes and conditions at once. More formally, let Xi,j denote the expression of gene i along condition j, let X Cokmeans i,j and X Rowkmeans i,j denote respectively the discretized values of Xi,j using the Cokmeans and the Rowkmeans techniques, with a parameter K + 1. Then the final Bikmeans discretized value of Xi,j is simply

X Rowkmeans i,j × X Colkmeans i,j
, where x = max{z ∈ Z | z ≤ x} denotes the floor function.

The three discretization techniques tested in this work were implemented using scikitlearn Python 3.7 library [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF],

5 Classification-based GRN inference

Definition

Classification and regression oriented GRN inference methods rely on similar principles. Indeed, in both cases the core idea is to train an algorithm to predict the expression level of a TG, from the expressions of a set of TFs. Then, the contribution of each TF to the prediction task is computed, as a feature importance score. Finally such scores can be directly used as proxies to quantify the dependency between the TG and each TF. The difference between the regression and classification oriented approaches is that while a regression algorithm can be trained directly on continuous gene expression data, it is necessary to discretize the TG expressions to form classes to train a classification algorithm. More formally, let tg (i) denote a TG, and let Xi,. be its continuous gene expression vector. Moreover let T F be the set of TFs (different from tg (i) ), and let XT F ,. denote their expression across all experimental conditions. Let y (i) = Disc(Xi,.) denote the discretized vector of gene expressions of tg (i) , where Disc : R J → Y J is simply a discretization function mapping each continuous element of a vector into a set Y = {1, 2, . . . , K}, where K denotes the number of classes (or bins). In practice any method presented in Section 4, can be used to discretize gene expressions. Let Φ θ : R |T F | → Y denote a classification model with parameters θ ∈ Θ, that aims at predicting, the discrete gene expressions y (i) j of tg (i) along each condition j, from the TFs gene expressions XT F ,j the same condition. In this context, each tf ∈ T F is considered as a predictive feature. In practice the classifier Φ θ is trained by finding a set of parameters θ that minimize some kind of average classification error, on the training dataset, i.e., argmin θ j error(y

(i) j ,Φ θ (X T F ,j )) J . Let us consider a function Γ : Θ, R |T F |×J , N J → R |T F | that computes γ (i) = Γ(θ, XT F ,.
, y (i) ) the vector of importance of predictive features (i.e., each tf ∈ T F ) relatively to the classification task achieved by Φ θ , from θ the parameters of the classification model, the training dataset XT F ,. and the target vector (i.e., discretize expressions of TG tg (i) ). Then, the importance of tf (j) in the prediction of y (i) is simply γ (i) j , and this measure can be used to approximate w(X,, tf (j) , tg (i) ) the dependency between tf (j) and tg (i) , as defined in Section 3.

Practical implementation

Target gene expression discretization Since classification algorithms require discrete target variables, we discretized each TG expression vector into K discrete levels (classes) of expression, using the k-means algorithm. The number of classes K was determined using the well-known elbow method: For different values K ranging from 2 to 30, we ran k-means to cluster the gene expressions into K clusters, and then we computed the average sum of square distance (SSE) between the gene expressions and the mean expression of their corresponding clusters. Finally we applied the Kneedle algorithm [START_REF] Satopaa | Finding a" kneedle" in a haystack: Detecting knee points in system behavior[END_REF], to locate the elbow in the plot representing the SSE for different number of clusters K. The elbow is a classical indicator that determines an appropriate number of clusters. Given the location of this indicator, as depicted in Figure 1, in this work the number of discrete levels of expression was set to K = 5. Classification algorithms In this work we decided to test five well-known classification algorithms, which are implemented in the scikit-learn Python 3.7 library [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF] (version 0.20.1). Four of these algorithms are ensemble methods based on decision trees, namely Random Forest (RF) [START_REF] Breiman | Random forests[END_REF], Extremely Randomized Trees (XRT) [START_REF] Geurts | Extremely randomized trees[END_REF], AdaBoost (AB) [START_REF] Freund | A short introduction to boosting[END_REF] and Gradient Boosting (GB) [START_REF] Friedman | Stochastic gradient boosting[END_REF]; and the last algorithm is the One-vs-All linear multi-class SVM classifier [START_REF] Cortes | Support-vector networks[END_REF].

For a decision tree, the importance of a feature d is computed by summing the weighted impurity decrease (e.g., GINI, entropy, variance) for each split defined along d; then for a set of trees, the importance of d is simply the average importance over all trees [START_REF] Breiman | Random forests[END_REF]. In practice, the feature importance was obtained from the dedicated scikit-learn model's feature importances attributes.

In the case of the One-vs-All linear Support Vector Machine (SVM), a linear binary SVM is trained to build a hyperplane that separates each class (i.e., TG's level of expression) form other classes [START_REF] Cortes | Support-vector networks[END_REF]. For a linear SVM, the absolute value of the coordinate of the hyperplane orthogonal vector, along a given feature d, is considered here as the importance of feature d; Since it indicates the participation of feature d, in the class separation, relatively to other features. For each feature (i.e., TF), its final importance is simply the average importance over all binary SVMs. In practice, the hyperplane orthogonal vector coordinates were obtained from the SVM scikit-learn classifier's coef attribute.

6 Experimental Setup

Datasets

All DREAM5 benchmark datasets used in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], contain a gene expression matrix, a list of putative TFs, and a gold standard GRN reporting known regulatory interaction between TFs and their TGs. Table 1 summarizes the major characteristics of this dataset, namely, the number of genes (I), the number of experimental conditions (J), the number of TFs (|T F |); as well as the major caracteristics of the gold standard such as its number of TFs ( E. coli In order to define the gold standard GRN for E. coli, the authors of [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF] collected a set of regulatory links with strong evidence, from the manually curated RegulonDB 6.8 database [START_REF] Gama-Castro | Regulondb version 7.0: transcriptional regulation of escherichia coli k-12 integrated within genetic sensory response units (gensor units)[END_REF]. The gene expression matrix provided in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], is a normalized Microarray dataset constituted of Affymetrix platform chips, downloaded from Gene Expression Omnibus (GEO)2 . According to [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], the Microarray dataset was normalized and filtered using: Robust Multichip Averaging [START_REF] Bolstad | A comparison of normalization methods for high density oligonucleotide array data based on variance and bias[END_REF], background adjustment, quantile normalization, probeset median polishing and, finally, the normalized gene expression values were transformed into a logarithmic scale. Finally, the list of TFs provided in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], is the union of the TFs list provided by RegulonDB 6.8 database [START_REF] Gama-Castro | Regulondb version 7.0: transcriptional regulation of escherichia coli k-12 integrated within genetic sensory response units (gensor units)[END_REF], and a list that the authors derived from Gene Ontology (GO) annotations.

S. cerevisiae

The S. cerevisiae gold standard GRN used in this paper was generated by MacIsaac et al. [START_REF] Macisaac | An improved map of conserved regulatory sites for saccharomyces cerevisiae[END_REF] by analyzing ChIP-chip datasets, and characterizing the presence of conserved TF binding sites motifs. The associated gene expression matrix provided in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], is a Microarray dataset constituted of Affymetrix platform chips downloaded from the GEO website. This Microarray dataset underwent the same normalization and filtering steps described for the E. coli dataset. The list of potential TFs provided by [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], is the union of the comprehensive list presented in [START_REF] Zhu | High-resolution dna-binding specificity analysis of yeast transcription factors[END_REF] and a set of TFs identified using GO terms.

S. aureus Unlike the previous datasets, no experimentally validated gold standard GRN was available for S. aureus. Thus, in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], the authors used RegPrecise [START_REF] Novichkov | Regprecise: a database of curated genomic inferences of transcriptional regulatory interactions in prokaryotes[END_REF], a database of regulatory interactions in prokaryotes, as a gold standard proxy to study the quality of the inference algorithms. As well as for the previous benchmark datasets, the S. aureus gene expression matrix provided in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], is a Microarray dataset constituted of Affymetrix platform chips downloaded from the GEO website. This Microarray dataset underwent the same normalization and filtering steps, that were described for the E. coli dataset. Finally a putative list of TFs was identified in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], considering the GO annotation of S. aureus genes.

In silico Unlike the previous benchmark datasets, this one has been generated by an in silico GRN model, generated using the fourth version of the GeneNetWeaver [START_REF] Schaffter | Genenetweaver: in silico benchmark generation and performance profiling of network inference methods[END_REF] software. According to [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], this synthetic GRN shares the core structure of the RegulonDB E. coli GRN, and incorporates 10% of random regulatory links. Finally the gene expression data was generated numerically, using this in silico GRN as a simulator, and relying on a dedicated dynamical model based on Ordinary Differential Equations, based on a multiplicative regulatory interactions hypothesis.

6.2 Evaluation Procedure

General procedure

In order to assess quality of the inference methods presented in this paper, with respect to gold standard GRNs, we used the evaluation framework presented in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF]. In this framewoork, GRN inference is evaluated as a binary classification task, which aims at predicting the presence of true regulatory links. Gold standard GRNs from real datasets, report lists of experimentally verified regulatory links between studied TFs and TGs. All these links are considered as true interactions, for the binary classification task. At first glance, one could think about considering all other possible links, that are not reported in the gold standards, as false interactions. Nevertheless, according to [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], gold standard links are only an incomplete subset of all the true regulatory interactions involved in the organism. Hence, all missing links should not be considered as false interactions. Indeed, one should avoid penalizing inference methods for detecting true interactions, that have not been experimentally verified yet. In order to choose a suitable set of false interactions, in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF] the authors decided to exclude any link involving a TF or a TG that has not been studied experimentally. And only the links between pairs of experimentally studied TF and TG, that have not been reported in the gold standards, are considered as false interactions, for the binary classification task.

More formally, let T G denote the set of genes of a given organism, and let T F ⊂ T G denote the subset of genes corresponding to TFs. Moreover, let G gold = T F gold ∪ T G gold , E gold be an oriented graph representing a gold standard GRN. Where T F gold ⊆ T F and T G gold ⊆ T G represent respectively the set of TFs and the set of TGs from the gold standard; While E gold represents the set of true regulatory interaction, such that ∀(tf gold , tg gold ) ∈ E gold , tf gold ∈ T F gold and tg gold ∈ T G gold denotes a true regulation of tg gold by tf gold . Let E f ull gold = {(e1, e2) ∈ T F gold ×(T F gold ∪T G gold ) | e1 = e2} be the set of all possible links between TFs and other genes from the gold standard (excluding self-loops). For the evaluation, all edges in E gold are considered as true interactions while edges in E f ull gold \ E gold are considered as false interactions. Notice that any interaction in (T F × T G) \ E f ull gold is not considered in the evaluation process.

Evaluation Measures

Evaluation measures As in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF] we assessed the methods using standard evaluation measures for binary classification, from the machine learning community, namely the Area Under the Receiver Operating Characteristic curve (AUROC) [START_REF] Fawcett | An introduction to roc analysis[END_REF], and the Area Under the Precision Recall curve (AUPR) [START_REF] Davis | The relationship between precision-recall and roc curves[END_REF] values.

p-values In order to make sense of the evaluation results, we proceeded to compare the performance of the inference methods with those that could be achieved by random predictions. To do so we computed empirical p-values for AUROC and AUPR values, following the procedure described in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], as described hereafter. For each gold standard, null empirical distributions of AUROC and AUPR values, termed respectively HAUROC and HAUP R, were computed for 25000 random GRNs. Such random GRNs were built by assigning random scores to all possible regulatory links. Then, the probability density functions of the empirical distributions were modeled by fitting a beta probability density function, in order to extrapolate the probability distributions beyond the empirical values range. Then for xAUROC and xAUP R, an AUROC value and an AUPR value respectively, the corresponding p-values are defined as p AU ROC = Pr(x ≥ xAUROC |HAUROC ) and p AU P R = Pr(x ≥ xAUP R|HAU P R).

Quality scores Finally, quality scores for each method along each dataset were computed by taking the opposite of the decimal logarithm of the p-values: score AU ROC = -log10(p AU ROC ) and score AU P R = -log10(p AU P R ). Notice that higher scores are associated to high evaluation measures, i.e., measures are less likely to be obtained by a random scoring.

Experimental protocol 6.4.1 Comparison with DREAM5 competitors

Gene expression data preprocessing Several methods that competed in the DREAM5 challenge, as well as state-of-the-art approaches (e.g., [START_REF] Irrthum | Inferring regulatory networks from expression data using tree-based methods[END_REF][START_REF] Haury | Tigress: trustful inference of gene regulation using stability selection[END_REF]), apply a simple z-score standardization to each row of the gene expression matrix, in order to center and scale the expression profile of each gene. In order to assess our method, we decided to standardize the benchmark gene expressions using the same preprocessing technique.

Parameter settings Each classification algorithm that has been used to score regulatory links relies on different meta parameters. One of the most important parameters of the methods based on sets of decision trees, is the number of such estimators that should be considered. Using many estimators, require more computational power, but tends to lead to better results. Here, we have set the number of estimators to 100 decision trees for each method. The remaining parameters were left to their default values, as established in their respective implementations on the 0.20.1 version of the scikit-learn Python 3.7 library [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF].

DREAM5 participants results

In [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF] the 35 methods that participated in the DREAM5 challenge were categorized in 6 families, namely 1) MI: data-driven methods based on Mutual Information, 2) Correlation: data-driven methods based on correlation metrics, 3) Regression: data-driven methods based on regression methods, 4) Bayesian: model-based Bayesian networks, 5) Meta: Multi-network methods, 6) Others: Different methods including model-based methods based on Boolean networks, or Genetic Algorithms. Combining the results of all the participants, the authors of [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF] inferred a single robust and high-quality GRN, which was termed Community. The performance measures obtained by each participant on each benchmark dataset, as defined in Section 6.2, have been made available by [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF]. Each one of our inference methods was also executed on the benchmark gene expression datasets, and the inferred GRNs were evaluated against the corresponding gold standards networks, using the procedure presented in Section 6.2.

Impact of preprocessing techniques

In this paper we also assessed the impact of the preprocessing techniques presented in Section 4, on the performance of the methods proposed here. The three standardiza-tion techniques do not need to be parameterized, while the number of bins (discrete classes) is the only parameter of the discretization techniques. In practice, we set the number of bins equal to five, since GRN inference methods tested in [START_REF] Li | Comparative study of discretization methods of microarray data for inferring transcriptional regulatory networks[END_REF], performed better with this number of bins. Then, each inference method was tested on the preprocessed datasets, and the resulting GRNs were evaluated with respect to the gold standards. Control experiments were also conducted by running the inference methods on raw data directly. In order to evaluate the impact of the preprocessing steps, we computed the gains in terms of quality measures, with respect to the control experiment.

All experiments were executed on a 2 GHz Intel Core i7 CPU running MacOS Mojave 10.14.5.

Experimental Results

Comparison with DREAM5 competitors

In order to assess our methods, we have computed the AUROC and the AUPR metrics, as well as the corresponding p-values and scores as described in Section 6.2. Then, we proceeded to compare our results to those reported in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF], for the different families of methods. As depicted in table 2, classification based methods output in average the best results in terms of AUROC and AUPR metrics and scores, even surpassing the average results obtained by the community approach. Moreover, the promising results of classification based methods does not seem to depend on a specific dataset, since their results are among the best for the four benchmark datasets, and they outcompete the community GRN for the three real organisms, as shown in Figures 2 and3.

As depicted in Figure 4, the different methods score differently depending on the benchmark dataset, and there is no ever-winning method among the classification based inference algorithms. For instance, the method based on SVM outputs the best results for E. coli and S. cerevisiae datasets both in terms of AUROC and AUPR, while it performs worse than the other four methods on the In silico dataset. An analogous phenomenon has also been observed in the case of the DREAM5 competitors, as reported in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF].

Impact of preprocessing techniques

As described in the Section 6.4.2, we have studied the impact of eight well-known preprocessing techniques, on the methods performance. The impact of these methods was evaluated by computing the gains in terms of AUROC and AUPR values and scores, with respect to control performance obtained on raw data. In this case, a negative (resp. positive) gain means that the results obtained on preprocessed data are worse (resp. better) than those obtained on raw data. The gains in AUROC range from -0.12 to 0.029, with an average gain equal to -0.013, and a standard deviation equal to 0.018. While the gains for AUPR range from -0.07 to 0.02, with an average gain equal to -0.019 and a standard deviation equal to 0.028. Therefore applying preprocessing techniques to these datasets tends to lead, in average, to a small degradation of the inference quality. Nevertheless, we note that these datasets already underwent a pipeline of sophisticated normalization and filtering steps, and consequently the impact of further preprocessings are likely to be less important. Average gains in AUROC and AUPR values, for each pair of inference method and preprocessing technique, are depicted as cluster-maps in Figure 5. This Figure shows that the vectors of gain of the four decision tree based methods tend to cluster together, showing that these techniques behave similarly to similar changes in the dataset, while the SVM based approach tends to behave rather differently. This last method tends to perform better on datasets that underwent row z-score or polishing standardizations. While decision tree based methods, obtained consistently the highest gains, on datasets that were discretized using the EFD technique. This last result is consistent to the conclusions of some studies (e.g., [START_REF] Jung | Evaluation of data discretization methods to derive platform independent isoform expression signatures for multi-class tumor subtyping[END_REF][START_REF] Kotsiantis | Discretization techniques: A recent survey[END_REF]) that affirm that discretizing continuous variables may lead some machine learning algorithms to produce more accurate models.

In

Conclusion

In this paper we presented a framework of data-driven GRN inference methods based on well-known classification algorithms, that includes five new GRN inference methods. These new techniques have been compared to well-established approaches on bench- mark datasets, using the evaluation procedure described in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF]. The results showed that inference methods based on classification algorithms exhibit satisfactory results, outperforming other families on average. These promising results suggest that a datadriven GRN inference paradigm relying on classification methods is an interesting complementary tool for the community. Future work perspectives include i) Running further analysis on RNAseq datasets, more complex organisms (e.g., [START_REF] Potier | Mapping gene regulatory networks in drosophila eye development by large-scale transcriptome perturbations and motif inference[END_REF]) and running a parameter sensitivity analysis, ii) Incorporating and studying new classification based methods, iii) Combining the results of several classification based techniques in an ensemble learning manner, to get more robust results as suggested in [START_REF] Marbach | Wis-dom of crowds for robust gene network inference[END_REF].

Fig. 1 :

 1 Fig. 1: Average sum of square distance (SSE) between gene expressions and clusters mean expression for different number of clusters K. The plot exhibits an elbow for K = 5 clusters.

Fig. 2 :

 2 Fig. 2: Boxplots representing AUROC and AUPR values for each dataset and each family of GRN inference methods. Classification based approaches are among the best for all datasets, and they out-compete the community GRN for the three real organisms.

Fig. 3 :

 3 Fig. 3: Boxplots representing AUROC and AUPR scores for each dataset and each family of GRN inference methods. Classification based approaches are among the best for all datasets, and they out-compete the Community GRN for the three real organisms.

Fig. 4 :

 4 Fig. 4: Barplots representing the AUROC and the AUPR associated to each classification algorithm, for each benchmark dataset.

  |T F gold |), TGs (|T G gold |) and regulatory links (|E gold |).Hereafter we describe briefly each benchmark dataset.

		Tab. 1: Benchmark datasets summary
	Dataset	Data	J	I	|T F | |T F gold | |T G gold | |E gold |
	In silico	Simulated 805 1,643	195	178	1,565	4,012
	S. aureus	Microarray 160 2,810	99	38	446	515
	E. coli	Microarray 805 4,511	334	141	1,081	2,066
	S. cerevisiae	Microarray 536 5,950	333	114	1,994	3,940

  Tab. 2: Average values and scores for AUROC and AUPR metrics, and number of methods for each GRN inference family.

	AUROC	0.6 0.7 0.8 0.9 1.0					Families Classification Meta MI Regression Correlation Others Bayesian Community Families Classification Meta MI Regression Correlation Others Bayesian Community
		0.5
		0.4	In silico	E. coli S. aureus S. cerevisiae Organism
	AUPR	0.2 0.3 0.4 0.5					Families Classification Meta MI Regression Correlation Others Bayesian Community Families Classification Meta MI Regression Correlation Others Bayesian Community
		0.1
		0.0	In silico	E. coli S. aureus S. cerevisiae Organism
				AUROC	AUPR	Number of
				Values Scores Values Scores	Methods
	Classification	0.67	70.05	0.18	95.60	5
	Community		0.64	35.62	0.13	47.94	1
	Others		0.58	9.97	0.06	8.54	8
	MI			0.60	9.53	0.09	9.09	5
	Meta			0.60	7.69	0.09	13.77	5
	Regression		0.59	7.17	0.09	18.70	8
	Correlation		0.59	4.65	0.08	5.65	3
	Bayesian		0.56	0.71	0.05	3.16	6

https://gitlab.com/speignier/classifiedgrni

http://www.ncbi.nlm.nih.gov/geo

Bikmeans


domain: spi phys math info

text:
Thibault Cantou
email: thibault.cantou@onera.fr

Nicolas Merlinge
email: nicolas.merlinge@onera.fr

Romain Wuilbercq
email: romain.wuilbercq@onera.fr

3DoF simulation model and specific aerodynamic control capabilities for a SpaceX's Starship-like atmospheric reentry vehicle

    3DoF simulation model and specific aerodynamic control capabilities for a SpaceX's Starship-like atmospheric reentry vehicle

Introduction

Since the beginning of the space era, the design of vehicles able to survive atmospheric reentry has been a challenging engineering problem. Reentry at orbital or near-orbital velocity implies very intense heat, aerodynamic pressure, and load factor. These issues are tackled with a combination of careful vehicle design (heat shielding or regenerative cooling, tough structure, etc.), and careful mission planning through reentry trajectory optimization if the vehicle has a maneuvering capability. The reentry trajectory can be divided in two phases: the exo-atmospheric phase, and the atmospheric phase. The latter is the most constraining one due to aerodynamic and aero-thermal effects at high velocity. During the first phase, the trajectory is controlled by exerting thrusts in various directions to track a desired velocity vector (orientation and norm). During the second phase, vehicles will often exploit angle of attack (AoA) and bank angle (BA) control to manoeuver and track a reference constrained trajectory minimizing heating, load factor, and aerodynamic pressure. Atmospheric trajectory optimization relies a lot on the aerodynamic capabilities of the considered vehicle. The first generations of reentry vehicles (i.e. reentry capsules) needed to be simple, robust to a great amount of uncertainties, and small enough to fit in the rockets of that era. The geometry allowing to deal with those constraints with the technology available at that time was not in favor of a high lift-to-drag ratio vehicles. They were nevertheless somewhat maneuverable during the atmospheric part of the flight (AoA and BA) thanks to center-of-mass (CoM) location control and/or cold gas thrusters for attitude control. Later, technological progress and the need for reusability allowed the design of rather safe reentry vehicles with more aerodynamically efficient shapes (i.e. with higher lift-to-drag ratio). The Space Shuttle and the X37-B are great examples of such vehicles. Their actuators (similar to those of an aircraft) allow this category of vehicle to remain stable at high AoA during the atmospheric flight. This then generate lift that can be oriented through BA control and used to shape the trajectory as needed.

The subject of this study, a novel design for a reusable atmospheric reentry vehicle, was proposed in September 2018: Space X's Starship (SXS). It exhibits an innovative configuration of aerodynamic control surfaces granting it unique capabilities. The motivation for this work is to give elements to understand why this configuration is an interesting choice. In the following sections, we will: Thibault CANTOU, Nicolas MERLINGE, Romain WUILBERCQ 2  Propose an OpenVSP [START_REF] Gloudemans | A rapid geometry modeler for conceptual aircraft[END_REF] geometrical model of SpaceX's Starship;  Introduce the methodology used to generate an aerodynamic database for lift, drag, and pitching moment at hypersonic speeds;  Analyze the implications of this design on longitudinal aerodynamic stability and its dependency to the location of the center of mass ;  Discuss the implication of the multi-actuator configuration on aerodynamic forces modulation while holding the vehicle at a fixed angle of attack ;  Provide some atmospheric reentry simulation results to illustrate the performances of the vehicle.

Context

SXS will be the second stage of SpaceX's future heavy rocket. It will be mounted on top of the Super Heavy booster. As it is a reusable second stage, SXS needs propulsion in order to reach the orbit, deorbit itself, and achieve a propulsive landing. But it also needs actuators able to control its attitude during atmospheric reentry. Since this study is focused on the behavior of SXS in the hypersonic regime during its reentry, it does not integrate a propulsive model. Figure 1 below presents the reference axis system and some elements of vocabulary used in this study. The actuated fins will only be used when separated from the booster and in the atmosphere at sufficient speed, i.e. during reentry. As demonstrated in section 4, this aerodynamic control configuration allows the vehicle to achieve aerodynamic stability at AoAs up to 90°, but not for values close to 0° (their rotation axis doesn't allow for torque generation at low AoA). As a matter of comparison, the Space Shuttle was able to control its AoA from 0° to 3DoF simulation model and specific aerodynamic control capabilities for a Space X's Starship-like atmospheric reentry vehicle. 3 approximately 40°. This very high AoA stability provides SXS with different manoeuvering capabilities that need to be explored.

SXS Model creation

3D geometry of the vehicle in OpenVSP

The 3D geometry of the vehicle was created in NASA's OpenVSP [START_REF] Gloudemans | A rapid geometry modeler for conceptual aircraft[END_REF] from photographs. As the vehicle is still in active development, the precise three-dimensional shape varies from a published picture to another. Therefore a representative model of the aerodynamic configuration of SXS has been created rather than a picture perfect model. As stated on SpaceX's website [START_REF] Space | Making Life Multiplanetary[END_REF], the vehicle's diameter is 9 m, and its length is 55 m including the fins. This allowed us to infer the rest of the required dimensions.

Body:

The body is 51 m long. The ogive section ends at 13 m along the X axis from the nose and then blends into a cylindrical shape until the base of the vehicle.

Actuated fins:

The selected profile is a symmetrical NACA 0010 profile. The planform is modelled the following way:

 Front fins: The front point of the front fins is located 2.6 m rearward the nose on the X axis. At the root, their chord is 3.5 m. At the tip, it is 2.2 m. The section tangent to the body forms an angle of 22° around the Z axis. They span 3 m orthogonally to the first section. The taper ratio is 0.63 and the sweep angle is 30°. The dihedral angle is the control angle and is limited to the range [-90°, 90°].  Rear fins: The front point of the rear fins is located 37.8 m rearward the nose on the X axis, at an angle of 30° down the XY plane. At the root, their chord is 11.4 m. At the tip, it is 3.7 m. They span 6.2 m orthogonally to the first section. The sweep angle is 64°. The dihedral angle is the control angle and is limited to [-60°, 120°].

Vertical fin: It is similar to the rear fins but the dihedral angle is kept at -90°. The convention for the fins control angles is chosen negative up and positive down. To account for the dependency of the aerodynamics to the actuated fins angles, a set of 3D models has been generated for every combination of rear and front fins angles in their respective definition domain with a step angle of 10°. From each of these models, a 3D mesh has been created in OpenVSP and then used in SHAMAN, as described over the next section.

Reduced-Order Hypersonic Aerodynamics

Future reusable launchers, such as those of interest in the present study, will face some severe aero-thermodynamic loads during their reentry into the terrestrial atmosphere. They will travel from the near-vacuum conditions of spacei.e. in the free-molecular flow regime -to the denser region of the lower atmosphere -i.e. in the continuum flow regimeand through a mixed-density transitional flow regime in between.

In the context of this study however, since the reentry trajectory of the Starship occurs mostly in regions where the flow behaves as a continuum, it is assumed that methods that are valid in this regime can be applied throughout the flight envelop of the vehicle. An in-house code, dubbed Shaman, is thus used in order to quantify the inviscid aerodynamics of our vehicle when it operates in the realm of the hypersonic flow regime. Shaman embeds a number of simplified methods, the so-called

Local Surface Inclination (LSI) methods [START_REF] Anderson | Hypersonic and High Temperature Gas Dynamics[END_REF]. Those only require, as inputs, a triangulation to describe the geometry of the vehicle being modelled (see Figure 4 below), the local inclination angle of all surface panels and the freestream conditions. The decrease in terms of physical realism that accompanies this panel-based approach, especially when compared to modern numerical methods such as Computation Fluid Dynamics (CFD), is offset almost entirely in many practical applications by its low computational cost and parametric flexibility [START_REF] Wuilbercq | Multi-disciplinary Modelling of Future Re-usable Space-Access Vehicles[END_REF]. Thus, albeit these LSI methods are considered of a low fidelity level, they provide reasonable estimates with regard to the preliminary nature of the analyses presented in this study.

Since the vehicle holds a very high angle-of-attack throughout its reentry trajectory and due to the general bluntness of the Starship, the modified Newtonian theory has been applied to all parts of the configuration geometry. Assuming a so-called Newtonian flow, the pressure exerted on the vehicle is thus considered solely due to the total loss of momentum of the fluid in the direction normal to the vehicle's surface [START_REF] Anderson | Hypersonic and High Temperature Gas Dynamics[END_REF]. The local pressure coefficient 𝐶 𝑝,𝑖 at any panel 𝑖 is then given by:

𝐶 𝑝,𝑖 = 𝐶 𝑝,𝑚𝑎𝑥 sin 2 𝜃 𝑖 (1)
Where 𝜃 𝑖 is the local inclination angle between the freestream flow direction and a tangent to the local surface panel.

As suggested by the work of Lester Lees [START_REF] Lees | Hypersonic Flow[END_REF], the modified Newtonian theory accounts for the total pressure loss across the normal shock that is located ahead of the vehicle. Therefore, 𝐶 𝑝,𝑚𝑎𝑥 is equal to the pressure coefficient at the stagnation point behind a normal shock.

A database of aerodynamic coefficients (i.e. 𝐶 𝐿 , 𝐶 𝐷 and 𝐶 𝑚/𝑟𝑒𝑓 ) has then been created as a function of the freestream Mach number 𝑀 ∞ (from Mach 4 to Mach 30), the vehicle's angle-of-attack 𝛼 (from 0° to 180°), and the front and rear control fins angles (respectively from -90° to 90°, and from -60° to 120°). Some results from the aerodynamic database without any deflection of the control surfaces (𝛿 = 0°) are presented in Figure 5.

The aerodynamic coefficients can then be interpolated in order to provide estimates at any particular point along the trajectory profile through the use of the reference properties associated with our SXS configuration (𝑆 𝑟𝑒𝑓 and 𝐿 𝑟𝑒𝑓 ) and the dynamic pressure. These reference quantities are provided in Table 1. It shall finally be borne in mind that viscous effects have been neglected in this preliminary work.

Reference Point for Moments [𝒎]

[-25, 0, 0] The origin is located at the nose of the vehicle 3DoF simulation model and specific aerodynamic control capabilities for a Space X's Starship-like atmospheric reentry vehicle. The aerodynamic forces and moment in wind axis are to be computed with the following expressions:

𝐹 𝐿𝑖𝑓𝑡 = -𝑄. 𝑆 𝑟𝑒𝑓 . 𝐶 𝐿 (2)
𝐹 𝐷𝑟𝑎𝑔 = -𝑄. 𝑆 𝑟𝑒𝑓 . 𝐶 𝐷

𝑀 𝑃𝑖𝑡𝑐ℎ = 𝑄. 𝑆. 𝑙. 𝐶 𝑚/𝑟𝑒𝑓 (4)

In the body axis system, the vector forces' expressions are:

𝐹 𝐴𝑥𝑖𝑎𝑙 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ = -𝑄. 𝑆 𝑟𝑒𝑓 . 𝐶 𝐴 . 𝑥 (5)
𝐹 𝑁𝑜𝑟𝑚𝑎𝑙 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ = -𝑄. 𝑆 𝑟𝑒𝑓 . 𝐶 𝑁 . 𝑧

With:

𝐶 𝐴 = 𝐶 𝐷 . cos(𝛼) -𝐶 𝐿 . sin (𝛼) (7)
𝐶 𝑧 = 𝐶 𝐷 . sin(𝛼) + 𝐶 𝐿 . cos (𝛼) (8)

Mass and center of mass

The vehicle's dry mass is estimated to be 85 tons [START_REF] Space | Abridged transcript of Elon Musk's presentation -Making Life Multiplanetary[END_REF]. The simulation results displayed in section 6 are based on this dry mass hypothesis, although this value excludes the fuel needed for the landing burn. A study including a propulsive model would be needed to estimate the required fuel mass needed for this part of the flight and then use a more representative mass during reentry simulation.

As demonstrated in the next section, the location of the center of mass has a great influence on the AoA solution for the pitch moment equation equilibrium. We assume that its coordinates in Y and Z are 0, and will only study the influence of variations of its X coordinate. Since the engines have a significant mass, and that most of the fuel has already been burnt before its reentry, the CoM is likely to be closer to the base than to the nose of the vehicle along the X axis. In the next section, we will study aerodynamic stability depending on the CoM's location. In section 6, the CoM will be considered still in the body frame and located 30 m behind the nose of the vehicle on the X axis.

Pitch equilibrium analysis and sensitivity to CoM location

Coefficient of moment at CoM

In order to study the influence of the CoM's location on the pitch equilibrium of the vehicle, the aerodynamic moment is expressed at the center of mass location (˄ being the cross-product operator):

𝑀 𝐶𝑜𝑀 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ = 𝑀 𝐶𝑅𝑀 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ + (𝐶𝑅𝑀 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ -𝐶𝑜𝑀 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) ∧ 𝐹 𝑎𝑒𝑟𝑜 ⃗⃗⃗⃗⃗⃗⃗⃗⃗ (9)
Given that this model is contained in the XZ plan, in the body frame the previous equation becomes:

𝑀 𝐶𝑜𝑀 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ | 𝑏𝑜𝑑𝑦 = 𝑄. 𝑆. 𝑙. 𝐶𝑚 𝐶𝑅𝑀 . 𝑦 + ( 𝑋 𝐶𝑅𝑀 -𝑋 𝐶𝑜𝑀 0 0 ) ∧ ( -𝐶𝐴 0 -𝐶𝑁 ) ( 10
)
With 𝑀𝑌 𝐶𝑜𝑀 being the second component of the total aerodynamic moment in the body frame, we have:

𝑀𝑌 𝐶𝑜𝑀 =

𝑪𝒎 𝑪𝒐𝑴 VS Mach

In this paragraph, the CoM is considered located at -27.5 m on the X axis. The following curves show the evolution of 𝐶𝑚 𝐶𝑜𝑀 against AoA for several Mach values, and for two sets of fins angles : 0° both (fins parallel to the XY plan), and -90° front and -60° rear (maximum fins angle above the XY plan, cf. middle picture of Figure 2):

𝑪𝒎 𝑪𝒐𝑴 VS CoM location

In this paragraph, we focus on the sensitivity of the pitching moment to the location of the CoM. In order to reduce the complexity of the analysis, we will only consider the two following extreme configurations of the fins angles (i.e. negative values of the angles):

 The first configuration is when the rear fins are at 0° (maximally exposed to the airflow) and the front fins are at 90° up (minimally exposed to the airflow). This configuration corresponds to the case of lowest equilibrium AoA. The rear fins act like the pitch stabilizers of an airplane's tail (cf. blue curve on Figure 7).

 The second configuration is the opposite: the rear fins are at 60° up (minimally exposed to the airflow) and the front fins are at 0° (maximally exposed to the airflow). This corresponds to the case of highest equilibrium AoA (cf. orange curve on Figure 7).

We deliberately limit the control fins angles to the negative part of their domain. This rather restrictive hypothesis is justified by the fact that high positive values of fins angles (i.e. down) when the vehicle is at high AoA implies a high risk of roll instability, while on the contrary high negative values (i.e. up) favor roll stability. As stated in paragraph 4.2, we can draw general conclusions with acceptable accuracy from the study of one point of Mach. We then arbitrarily choose Mach 15 in the following results. For the two considered fins angles configurations, the following figures display: the 𝐶𝑚 𝐶𝑜𝑀 versus AoA for several positions of CoM on the left, and the obtained equilibrium AoA versus the CoM's location on the right. Propulsion systems protection: As it is not desirable to frontally expose the nozzles to the hypersonic airflow during reentry, it is necessary for the vehicle's AoA to be in a stable equilibrium between 0° and 90° during the flight. It is also required to minimize the thermal flow. From this standpoint, it is optimal to fly with an AoA comprised between 0° and 90°.

Global heating minimization:

Since the thermal flow highly depends on the curvature of the exposed surface, it reaches a maximum value at 0° AoA, and a minimum value at 90°. From this standpoint, the optimum value would be a high AoA during reentry, likely around 90°.

Aerodynamic manoeuver capability for trajectory shaping:

The vehicle needs to manoeuver aerodynamically in the atmosphere by generating lift for several reasons: it allows for trajectory optimization, and it also helps to ensure flight safety thanks to closed-loop trajectory control during reentry. The maximum lift generation occurs with an AoA near 50°, independently from the fins angles. The minimum lift occurs at 0° and 90°. Then, in order to fully modulate the lift, the controllable AoA domain must span from 50° to 90°. Those constraints are satisfied with a minimum AoA domain between 50° and 90°, which is ensured for a CoM location between -27 m and -30 m. As stated earlier, the heaviest systems (the propulsion systems) will likely be closer to the base than the nose. For this reason, we will now consider the CoM to be still and located at -30 m.

3DoF simulation model and specific aerodynamic control capabilities for a Space X's Starship-like atmospheric reentry vehicle. 9

Lift and drag modulation capacity at a constant equilibrium AoA

Available control domain for fins angles at iso-equilibrium AoA

Since this vehicle has two sets of actuators, there are several combinations of front and rear fins angles that cancel out aerodynamic moment. Therefore, an equilibrium AoA can be obtained for several combinations of controls. Figure 8 displays the functions 𝛿 𝑟 = 𝑓(𝛿 𝑓 , 𝐴𝑜𝐴 0 ) with the constraint of null aerodynamic moment at the CoM (at Mach 15 and altitude 60 km). It corresponds to a family of all the couples (𝛿 𝑟 , 𝛿 𝑓 ) that stabilize the aerodynamic moment for a given AoA 𝛼 0 .

Figure 8 : Front and rear fins angles domain for several equilibrium AoAs

For exemple, at AoA = 90°, the rear fins control angle can vary from -41° to -27° for front fins control angles ranging from -90° to 0°.

Lift and drag modulation at an equilibrium AoA

Since the fins have a coupled effect on rotational and translational dynamics of the vehicle, the available controls combinations domain at iso-equilibrium AoA allows for modulation of lift and drag without changing the vehicle's attitude. Figure 9 displays the limits of lift and drag acceleration generated by the vehicle at 60 km altitude and Mach 15, versus the desired equilibrium AoA.

As it is visible on these curves, the modulation of drag has greater effects for an AoA near 90°. Concerning lift, the modulation is greater at an AoA close to 65°. The modulation capability is about 12% at best for the lift and drag at their respective preferred AoA. While a 12% modulation does not yield great variations of the aerodynamic forces, it can have a noticeable impact on the final ground distance, as it will be demonstrated on the simulation results presented in the next section.

It is of importance to note that this model does not include friction drag. As a consequence, the results presented here depend on altitude only because of the dynamic pressure in the forces and moment computation. Further work will account for the effect of friction drag on the aerodynamic database and system performance.

Figure 9 : Limits of lift and drag modulation capability versus stable equilibrium AoA

Simulation of translational dynamics

The translational dynamics part of the model has been extracted and used to simulate reentry trajectories of the vehicle in the vertical plan. As it is common for this kind of study, rotational dynamics has been neglected in those simulations for two main reasons: the lack of information about the inertia tensor and pitch damping coefficient 𝐶𝑚 𝑞 , and that accurate rotational dynamics simulation often require a much shorter time step than the one use for translational dynamics. However, the result showed here are still accurate because they can be viewed as simulations of a nominally controlled vehicle in terms of AoA. We use the iso-equilibrium AoA curves of the previous section to select the fins angles to apply depending on the desired AoA. The AoA is considered constant in those simulations.

The reference frame is the local geographic trihedral (North -East -Down, NED). The gravitational acceleration is simply computed as follows:

𝑔 = 𝜇 (𝑅 𝑇 + ℎ) 2 (14)
With 𝜇 the Earth's standard gravitational parameter (398 600 km 3 .s -2 ), 𝑅 𝑇 the Earth's radius (6 378.139 km), and ℎ the current altitude of the vehicle obtained as follows:

ℎ = -𝑧 (1)
With z the coordinate of the vehicle along the Z axis of the North-East-Down frame.

To account for orbital effects of speed, the vehicle is subject to an apparent gravitational acceleration in the reference frame equal to the summation of 𝑔 and the centrifugal acceleration. This centrifugal acceleration is computed as follows:

𝐴𝑐𝑐 𝑐𝑡𝑓𝑔 = - (𝑉. cos(𝛾)) 2 𝑅 𝑇 + ℎ (2)
The atmosphere model used is the US76 model. The Earth is considered still, which means the aerodynamic speed is the same as the inertial speed. Because the aerodynamic characterization method is limited to hypersonic regime, the simulations are interrupted when the vehicle's speed becomes inferior to Mach 4 (the implemented speed limit corresponds to the Mach number reached at 35 km altitude).

3DoF simulation model and specific aerodynamic control capabilities for a Space X's Starship-like atmospheric reentry vehicle.

11

The state derivatives equations are the following:

( 𝑥ż V γ ̇) = ( 𝑉.
With 𝑥 and 𝑧 the coordinates of the vehicle in the reference North-East-Down frame, 𝑉 its speed, 𝛾 its flight path angle, 𝑚 its mass, and 𝐹 𝐷𝑟𝑎𝑔 and 𝐹 𝐿𝑖𝑓𝑡 respectively the drag force and lift force as explained in paragraph 3.2. The state is initialized with the following values:

Reentry with different AoAs

To illustrate the diverse aerodynamic capabilities of the vehicle, we computed the trajectories for two constant AoAs during reentry. The first case corresponds to an AoA of 90° (negligible lift), and the second one to an AoA 50° (high lift). The control fins angles have been selected on the curves displayed in Figure 8. Their value is selected on the corresponding AoA curve at the point corresponding to the middle of the 𝛿 𝑓 available domain: in the first case we choose 𝛿 𝑓 = -45° and 𝛿 𝑟 = -36°, and in the second case we choose 𝛿 𝑓 = -75° and 𝛿 𝑟 = -8°. The Figure 10 displays the obtained results for those two simulations.

As it is demonstrated on the Altitude VS 𝑥 plot, the vehicle's ability to fly with an AoA between 90° and 50° will allow for very different reentry trajectories. The load factor is drastically lower on the AoA 50° trajectory (about 1.6 g) than on the AoA 90° trajectory (almost 9 g). It is explained by the fact that aerodynamic lift allows the vehicle to stay longer in the thin high atmosphere which allows it as a consequence to reduce less fast its kinetic energy than in the case without lift. Manned transportation applications and structure design requirements urge to select a reentry trajectory with an AoA allowing for lift generation.

Note that these conclusions heavily rely on initial reentry conditions. With the selected set of initial conditions, in the AoA 50° case, the vehicle generates enough lift force to engage phugoid motion, which is not desirable as we look for a smooth and steady reentry. Simulations with more representative initial conditions will be part of future studies to ensure that phugoid motion still happens with an AoA of 50°. An ideal constant AoA would then be found between 90° and 50° at the lowest AoA without phugoid motion. Future work will also include trajectory optimization using the AoA as time-dependent control parameter.

Thibault CANTOU, Nicolas MERLINGE, Romain WUILBERCQ

Reentry at a constant AoA with drag modulation

To illustrate the modulation of aerodynamic force presented in the previous section, we compute the two extreme feasible trajectories for an AoA of 90°. As in the previous paragraph, the control fins angles have been selected on the curves displayed in Figure 8. The first corresponds to fins angles of 𝛿 𝑓 = -90° and 𝛿 𝑟 = -41° (minimum aerodynamic force), and in the second to 𝛿 𝑓 = 0° and 𝛿 𝑟 = -27° (maximum aerodynamic force). Those two control configurations produce sensible difference only in terms of final 𝑥 distance. The first trajectory stops 120 km further than the second one.

Conclusion

This study's goal was to render available some data and insights about the aerodynamics of SpaceX's Starship vehicle class. The article presents an estimation of the vehicle's non-viscous aerodynamic coefficients in the hypersonic flow regime dependent on Mach number, angle of attack, and the actuated fins angles. This estimation relies on 3D meshes generated using the OpenVSP tool, and the modified Newtonian theory for local pressure computation. The aerodynamics coefficients reveal little dependence on the Mach number.

An analysis of the longitudinal stability depending on the center of mass' location is also presented. It is shown that this parameter is of great influence on the limits of the AoA domain stabilizable by the actuated fins. As a noticeable conclusion, we observe that a displacement of the center of mass from 30 m to 32 m behind the nose implies a drift of the lower stabilizable AoA boundary from 50° to 80°, which deteriorates drastically the aerodynamic maneuvering capacity of the vehicle.

With the selected set of initial reentry parameters, simulation results at constant AoAs of 90° (non-lifting trajectory) and 50° (lifting trajectory) demonstrates a reduction of the peak load factor from 9 g to 1.6 g. This advocates for lifting reentry trajectories, but will require a more extensive study of reentry trajectories with a more representative set of initial velocity and flight path angle.

On a fixed angle of attack, the multi-actuator configuration (front and rear fins) allowing for aerodynamic force modulation of 12% at best, simulations result demonstrates on the case of study little influence on reentry parameters. We observe a final distance difference of 4% between the two extreme cases.

This preliminary work has arisen numerous questions and needs for improvement of the 3DoF model. Thus, future work will include:

 Reentry simulations with more representative initial conditions ;  Adding a thermal flux computation model ;  Trajectory optimization using the AoA as a time-dependent control parameter in order to minimize thermal flux and load factor.

 Estimation of the inertia tensor to allow for rotational dynamics simulation and accurate attitude algorithm design ;  Creation of a more analytic model fitted on the generated database in order to reduce the computational cost of interpolating the aerodynamic coefficients. Identifying functional dependency of the aerodynamic coefficients to the flight parameters will allow for easier control studies, particularly in the non-linear control field.

Figure 1 :

 1 Figure 1 : Body frame and parts of the SXS model

Figure 2 :

 2 Figure 2 : Degrees of freedom of the actuated fins

Figure 3 :

 3 Figure 3 : Up view (left) and right-side view (right) of the OpenVSP model

Figure 4 :

 4 Figure 4 : Unstructured triangular mesh of the Starship (made up of 16800 panels)

Figure 5 :

 5 Figure 5 : Aerodynamic Coefficients at 𝛿 𝑟 = 0° and 𝛿 𝑓 = 0°

Figure 6 :

 6 Figure 6 : Variation of 𝐶𝑚 𝐶𝑜𝑀 versus Mach number for a CoM located at -27.5 m on the X axis

Figure 7 :

 7 Figure 7 : Variation of 𝐶𝑚 𝐶𝑜𝑀 versus CoM location along the X axis for the two extreme configurations of fins angles

Figure 10 :

 10 Figure 10 : Illustration of reentries with constant AoAs of 90° (blue curve) and 50° (red curve)

Figure 11 :

 11 Figure 11 : Illustration of the aerodynamic force modulation capability using two control fins angle configurations for a stable AoA of 90°: 𝛿 𝑓 = -90° and 𝛿 𝑟 = -41° (blue curve), and 𝛿 𝑓 = 0° and 𝛿 𝑟 = -27° (red curve).

Table 1 :

 1 Reference Properties for the Aerodynamics of the Starship

	Reference	Value	Remark
	Length, 𝑳 𝒓𝒆𝒇 [𝒎]	9	Fuselage Diameter
	Area, 𝑺 𝒓𝒆𝒇 [𝒎 𝟐 ]	64	Fuselage Circular Cross-Section

  𝐶𝑚 𝐶𝑜𝑀 is a function of Mach, AoA, CoM location, 𝛿 𝑟 , and 𝛿 𝑓 . Solving for the AoA, the zeros of 𝐶𝑚 𝐶𝑜𝑀 are the equilibrium points of the pitching moment equation. Moreover, an equilibrium point is open-loop locally stable if the slope of 𝐶𝑚 𝐶𝑜𝑀 is negative. For safety reasons, open-loop rotational stability of the vehicle during atmospheric reentry is a desirable behavior.

	3DoF simulation model and specific aerodynamic control capabilities for a Space X's Starship-like atmospheric
	reentry vehicle.
	𝑀𝑌 𝐶𝑜𝑀 = 𝑄. 𝑆. 𝑙. 𝐶𝑚 𝐶𝑜𝑀		(12)
	With 𝐶𝑚 𝐶𝑜𝑀 :
	𝐶𝑚 𝐶𝑜𝑀 = 𝐶𝑚 𝐶𝑅𝑀 +	𝑋 𝐶𝑅𝑀 -𝑋 𝐶𝑜𝑀 𝑙	. 𝐶𝑁	(13)
				(11)
	7

𝑄. 𝑆. (𝑙. 𝐶𝑚 𝐶𝑅𝑀 + (𝑋 𝐶𝑅𝑀 -𝑋 𝐶𝑜𝑀 ). 𝐶𝑁)

Table 2 :

 2 Initial state values used in simulation

	State parameter	Initial value
	𝑥	0 m
	𝑧	-90 000 m
	𝑉	7 850 m/s
	𝛾	0 °

Acknowledgements

We would like to thank Mathieu Balesdent and Loïc Brevault for their constructive remarks that have enriched this work, and ONERA's Space Programs Direction to provide us the financial means to attend EUCASS 2019.