Sadaf Ul Zuhra
email: sadaf-ul.zuhra@inria.fr

Samir M Perlaza
email: samir.perlaza@inria.fr

Eitan Altman
email: eitan.altman@inria.fr

In this paper, the fundamental limits on the rates at which information and energy can be simultaneously transmitted over an additive white Gaussian noise channel are studied under the following assumptions: (a) the channel is memoryless; (b) the number of channel input symbols (constellation size) and block length are finite; and (c) the decoding error probability (DEP) and the energy outage probability (EOP) are bounded away from zero. In particular, it is shown that the limits on the maximum information and energy transmission rates; and the minimum DEP and EOP, are essentially set by the type induced by the code used to perform the transmission. That is, the empirical frequency with which each channel input symbol appears in the codewords. Using this observation, guidelines for optimal constellation design for simultaneous energy and information transmission are presented.

I. INTRODUCTION

. In [4], a capacity-energy function is defined in order to determine the fundamental limit on the information transmission rate subject to the fact that the average energy at the channel output is not smaller than a given threshold. Therein, the underlying assumption is that the communication duration in channel uses is infinitely long. This guarantees that the decoding error probability (DEP) and the energy outage probability (EOP) can be made arbitrarily close to zero, and thus, the focus is only on the information transmission rate and the energy transmission rate. This analysis has been extended to multi-user channels. In this case, the notion of information-energy capacity region generalizes to the set of all information and energy rate

tuples that can be simultaneously achieved in the asymptotic block length regime [START_REF] Amor | Fundamental limits of simultaneous energy and information transmission[END_REF]. For instance, the information-energy capacity region of the Gaussian multiple access channel is characterized in [START_REF] Amor | Feedback enhances simultaneous energy and information transmission in multiple access channels[END_REF], whereas the information-energy capacity region of the Gaussian interference channel is approximated in [START_REF] Khalfet | Simultaneous information and energy transmission in the two-user Gaussian interference channel[END_REF].

A first attempt to study the fundamental limits of SIET under the assumption of finite transmission duration with finite constellation sizes is presented in [START_REF] Perlaza | Simultaneous information and energy transmission: A finite block-length analysis[END_REF]. Therein, the study is restricted to discrete memoryless channels and it is shown that finite transmission duration implies DEPs and EOPs that are bounded away from zero. This work contributes in this direction and considers the problem of SIET in additive white Gaussian noise channels considering finite transmission duration and finite constellation sizes. The main results in this work highlight the intuition that codes that uniformly use all channel input symbols are associated with high information rates, whereas, codes that exclusively use the channel input symbols that carry the largest amount of energy are associated with high energy rates. More specifically, a characterization of the maximum information and energy rates and minimum DEP and EOP that can be simultaneously achieved is formulated in terms of the type the code induces on the set of channel input symbols. In this work, a type is understood in the sense of the empirical frequency with which each channel input symbol appears in the codewords [START_REF] Csiszár | The method of types[END_REF].

II. SYSTEM MODEL Consider a communication system formed by a transmitter, an information receiver (IR), and an energy harvester (EH). The objective of the transmitter is to simultaneously send information to the IR at a rate of R bits per second; and energy to the EH at a rate of B Joules per second over an additive white Gaussian noise (AWGN) channel. That is, given a channel input x = (x 1 , x 2 , . . . , x n ) T ∈ C n , with n ∈ N, the outputs of the channel are the random vectors

Y =x + N 1 , and (1a) Z=x + N 2 , ( 1b
)
where n is the duration of the transmission in channel uses; and the vectors

Y = (Y 1 , Y 2 , . . . , Y n ) T ∈ C n and Z = (Z 1 , Z 2 , . . . , Z n ) T ∈ C n
are the inputs of the IR and the EH, respectively. The components of the random vectors

N 1 = (N 1,1 , N 1,2 , . . . , N 1,n ) T ∈ C n and N 2 = (N 2,1 , N
f Z|X (z|x)= n t=1 f Z|X (z t |x t ), (3)
and for all t ∈ {1, 2, . . . , n},

f Y |X (y t |x t ) = 1 πσ 2 exp Å - ( (y t ) -(x t )) 2 + ( (y t ) -(x t )) 2 σ 2 ã , (4a)
f Z|X (z t |x t ) = 1 πσ 2 exp Å - ( (z t ) -(x t )) 2 + ( (z t ) -(x t )) 2 σ 2 ã . (4b)
Within this framework, two tasks must be accomplished: information transmission and energy transmission.

A. Information Transmission

Assume that the information transmission takes place using a modulation scheme that uses L symbols. That is, there is a set X {x (1) , x (2) , . . . , x (L) } ⊂ C

that contains all possible channel input symbols, and L |X |.

Let M be the number of message indices to be transmitted within n channel uses. That is, M 2 n log 2 L To reliably transmit a message index, the transmitter uses an (n, M )-code defined as follows.

Definition 1. (n, M )-code: An (n, M )-code for the random transformation in (1) is a system:

{(u(1), D 1 ), (u(2), D 2 ), . . . (u(M ), D M )}, (6)
where, for all (i, j) ∈ {1, 2, . . . , M } 2 , i = j,

u(i) = (u 1 (i), u 2 (i), . . . , u n (i)) ∈ X n , (7a) D i ∩ D j = φ, (7b)
M i=1 D i ⊆ C n , and (7c)
|u t (i)| P, ( 7d
)
where P is the peak-power constraint. Assume that the transmitter uses the (n, M )-code

C {(u(1), D 1 ), (u(2), D 2 ), . . . , (u(M ), D M )}, (8)
that satisfies [START_REF] Perera | Simultaneous wireless information and power transfer (SWIPT): Recent advances and future challenges[END_REF]. The information rate of any (n, M )-code is given by

R = log 2 M n (9)
in bits per channel use. To transmit the message index i, with i ∈ {1, 2, . . . , M }, the transmitter uses the codeword u(i) = (u 1 (i), u 2 (i), . . . , u n (i)). That is, at channel use t, with t ∈ {1, 2, . . . , n}, the transmitter inputs the symbol u t (i) into the channel. At the end of n channel uses, the IR observes a realization of the random vector Y = (Y 1 , Y 2 , . . . , Y n ) T in (1a). The IR decides that message index j, with j ∈ {1, 2, . . . , M }, was transmitted, if the event Y ∈ D j takes place, with D j in [START_REF] Huang | Simultaneous wireless information and power transfer: Technologies, applications, and research challenges[END_REF]. That is, the set D j ∈ C n is the region of correct detection for message index j. Therefore, the DEP associated with the transmission of message index i is

γ i (C ) 1 - Di f Y |X (y|u(i))dy, (10)
and the average DEP is

γ(C ) 1 M M i=1 γ i (C ). (11)
Using this notation, Definition 1 can be refined as follows.

Definition 2 ((n, M, )-codes). An (n, M )-code for the random transformation in (1), denoted by C , is said to be an

(n, M, )-code if γ(C ) < . (12)

B. Energy Transmission

Let g : C → [0, +∞] be a positive function such that given a channel output z ∈ C, the value g(z) is the energy harvested from such channel output. The energy transmission task must ensure that a minimum average energy B is harvested at the EH at the end of n channel uses. Let ḡ : C n → [0, +∞] be a positive function such that given n channel outputs z = (z 1 , z 2 , . . . , z n ), the average energy is

ḡ(z) = 1 n n t=1 g(z t ), (13)
in energy units per channel use. Assume that the transmitter uses the code C in [START_REF] Huang | Simultaneous wireless information and power transfer: Technologies, applications, and research challenges[END_REF]. Then, the EOP associated with the transmission of message index i, with i ∈ {1, 2, . . . , M }, is

θ i (C , B) Pr[ḡ(Z) < B|X = u(i)], (14)
where the probability is with respect to the probability density function f Z|X in (3); and the average EOP is given by

θ(C , B) 1 M M i=1 θ i (C , B). ( 15
)
This leads to the following refinement of Definition 2.

Definition 3 ((n, M, , B, δ)-code). An (n, M, )-code for the random transformation in (1), denoted by C , is said to be an

(n, M, , B, δ)-code if θ(C , B) < δ. (16)
III. MAIN RESULTS The results in this section are presented in terms of the types induced by the codewords of a given code. Given an (n, M, , B, δ)-code denoted by C of the form in [START_REF] Huang | Simultaneous wireless information and power transfer: Technologies, applications, and research challenges[END_REF], the type induced by the codeword u(i), with i ∈ {1, 2, . . . , M }, is a probability mass function (pmf) whose support is X in [START_REF] Amor | Fundamental limits of simultaneous energy and information transmission[END_REF]. Such pmf is denoted by P u(i) and for all ∈ {1, 2, . . . , L},

P u(i) (x ( ) ) 1 n n t=1 1 {ut(i)=x ( ) } , (17)
where x ( ) is an element of X in [START_REF] Amor | Fundamental limits of simultaneous energy and information transmission[END_REF]. The type induced by all the codewords in C is also a pmf on the set X in [START_REF] Amor | Fundamental limits of simultaneous energy and information transmission[END_REF]. Such pmf is denoted by P C and for all ∈ {1, 2, . . . , L},

P C (x ( ) ) 1 M M i=1 P u(i) (x ( ) ). ( 18
)
Using this notation, the main results are essentially upper bounds on the information rate R in ( 9) and energy rate B in ( 14); as well as, lower bounds on the average DEP and average EOP δ for all possible (n, M, , B, δ)-codes. These bounds are provided assuming that the L channel input symbols in X and the block length n are finite.

A. Energy Transmission Rate and Average EOP

The following lemma introduces an upper bound on the energy transmission rate that holds for all possible (n, M, , B, δ)codes.

Lemma 1. Given an (n, M, , B, δ)-code denoted by C for the random transformation in (1) of the form in [START_REF] Huang | Simultaneous wireless information and power transfer: Technologies, applications, and research challenges[END_REF], the energy transmission rate B satisfies,

B ≤ 1 1 -δ L =1 P C (x ( ) )E î g Ä x ( ) + W äó , ( 19
)
where P C is defined in (18) and the expectation is with respect to W , which is a complex circularly symmetric Gaussian random variable whose real and imaginary parts have zero means and variances 1 2 σ 2 . Proof: The proof of Lemma 1 is presented in [START_REF] Zuhra | Simultaneous information and energy transmission with finite constellations[END_REF].

An observation from Lemma 1 is that the limit on the energy rate B does not depend on the individual types P u(1) , P u(2) , . . ., P u(M ) but on the type induced by all codewords, i.e., P C . The inequality in (19) implies the following corollary.

Corollary 1. Given an (n, M, , B, δ)-code denoted by C for the random transformation in (1) of the form in (8), the average EOP δ satisfies,

δ ≥ 1 - 1 B L =1 P C (x ( ) )E î g(x ( ) + W ó + . ( 20
)
where P C is defined in (18) and the expectation is with respect to W , which is a complex circularly symmetric Gaussian random variable whose real and imaginary parts have zero means and variances 1 2 σ 2 .

B. Average Decoding Error Probability

The analysis of the information transmission rate and the average DEP of a given code C of the form in (8) depends on the choice of the decoding sets D 1 , D 2 , . . ., D M . Without any loss of generality, assume that for all i ∈ {1, 2, . . . , M }, the decoding set D i is written in the form

D i = D i,1 × D i,2 × . . . × D i,n , (21)
where for all t ∈ {1, 2, . . . , n}, the set D i,t is a subset of C.

The following lemma introduces a lower bound on the average DEP that holds for all possible (n, M, , B, δ)-codes.

Lemma 2. Given an (n, M, , B, δ)-code C for the random transformation in (1) of the form in (8), the average DEP satisfies,

≥1 - 1 M M i=1 exp -nH(P u(i) ) -nD(P u(i) ||Q) +n log L j=1 Ej f Y |X (y|x (j) )dy , (22)
where, P u(i) is the type defined in (17); the function Q is a pmf on X in (5) such that for all i ∈ {1, 2, . . . , L},

Q Ä x (i) ä = Ei f Y |X (y|x (i) )dy L j=1 Ej f Y |X (y|x (j) )dy ; (23)
and for all ∈ {1, 2, . . . , L}, the set E is

E = D i ,t , (24)
where (i , t ) satisfies (i , t ) arg min (i,t)∈{1,2,...,M }×{1,2,...,n} Di,t f Y |X (y|x ( ) )dy.

(25)

Proof: The proof of Lemma 2 is presented in [START_REF] Zuhra | Simultaneous information and energy transmission with finite constellations[END_REF].

A class of codes that is of particular interest in this study is that of homogeneous codes, which are defined hereunder.

Definition 4 (Homogeneous Codes

). An (n, M, , B, δ)-code denoted by C for the random transformation in (1) of the form in (8) is said to be homogeneous if for all i ∈ {1, 2, . . . , M } and for all ∈ {1, 2, . . . , L}, it holds that

P u(i) (x ( ) ) = P C (x ( ) ), (26)
and for all (i, t) ∈ {1, 2, . . . , M } × {1, 2, . . . , n} for which u t (i) = x ( ) , for some ∈ {1, 2, . . . , L}, it holds that,

D i,t = E , (27)
where, P u(i) and P C are the types defined in (17) and (18), respectively; the sets D i,t are defined in (21); and the sets E are defined in (25).

Homogeneous codes are essentially (n, M, , B, δ)-codes that satisfy two conditions. First, a given channel input symbol is used the same number of times in all codewords; and second, every channel input symbol is decoded with the same decoding set independently of the codeword and/or the position in the codeword. Lemma 2 simplifies for the case of homogeneous codes (Defintion 4) as follows.

Corollary 2. Given a homogeneous (n, M, , B, δ)-code denoted by C for the random transformation in (1) of the form in (8), the DEP satisfies that

≥1 -exp -nH(P C ) -nD(P C ||Q) +n log L j=1 Ej f Y |X (y|x (j) )dy , (28)
where, P C is the type defined in (18); the pmf Q is defined in (23) and for all j ∈ {1, 2, . . . , L}, the sets E j are defined in (25).

Note that the right-hand side of ( 28) is minimized when the following condition is met,

P C = Q. ( 29
)
This observation leads to the following corollary.

Corollary 3. Given a homogeneous (n, M, , B, δ)-code denoted by C for the random transformation in (1) of the form in [START_REF] Huang | Simultaneous wireless information and power transfer: Technologies, applications, and research challenges[END_REF], such that (29) is satisfied, the DEP satisfies that

1 -exp Ñ -nH(Q) + n log Ñ L j=1 Ej f Y |X (y|x (j) )dy é é , (30)
where the pmf Q is defined in (23) and for all j ∈ {1, 2, . . . , L}, the sets E j are defined in (25).

Note that the equality in (29) reveals the optimal choice of P C with respect to the average DEP, which is using the channel input symbols in X with a probability that is proportional to the probability of correct decoding.

C. Information Transmission Rate

A first upper bound on the information rate is obtained by upper bounding the number of codewords that a code might possess given the particular types P u(1) , P u(2) , . . ., P u(n) ; or the average type P C in (18). The following lemma introduces such a bound for the case of a homogeneous code.

Lemma 3. Given a homogeneous (n, M, , B, δ)-code denoted by C for the random transformation in (1) of the form in (8), the information transmission rate R in (9) is such that

R ≤ 1 n log 2 Ç n! L =1 (nP C (x ( ) ))! å ≤ log 2 L, (31)
where, P C is the type defined in (18).

Proof: The proof of Lemma 3 is presented in [START_REF] Zuhra | Simultaneous information and energy transmission with finite constellations[END_REF].

Lemma 3 can be written in terms of the entropy of the type P C as shown by the following corollary.

Corollary 4. Given a homogeneous (n, M, , B, δ)-code denoted by C for the random transformation in (1) of the form in [START_REF] Huang | Simultaneous wireless information and power transfer: Technologies, applications, and research challenges[END_REF], the information transmission rate R in ( 9) is such that

R≤H (P C ) + 1 n 2 1 12 - L =1 1 12P C (x ( ) ) + 1 + 1 n log Ä√ 2π ä - L =1 log » 2πP C (x ( ) ) - log n n Å L -1 2 ã , (32)
where, P C is the type defined in (18).

Proof: The proof of Corollary 4 is presented in [START_REF] Zuhra | Simultaneous information and energy transmission with finite constellations[END_REF].

Note that all terms in (32), except the entropy H (P C ), vanish with the block length n. This implies that the information rate is essentially constrained by the entropy of the channel input symbols. In particular, note that H (P C ) log 2 L.

D. Information and Energy Trade-Off

The results presented above lead to the following result.

Theorem 1. Given a homogeneous (n, M, , B, δ)-code denoted by C for the random transformation in (1) of the form in (8), the following holds

R≤H (P C ) + 1 n 2 1 12 - L =1 1 12P C (x ( ) ) + 1 + 1 n log Ä√ 2π ä - L =1 log » 2πP C (x ( ) ) - log n n Å L -1 2 ã ; (33)
B≤ 1 1 -δ L =1 P C (x ( ) )E î g Ä x ( ) + W äó ; and (34)
≥1 -exp -nH(P C ) -nD(P C ||Q) +n log L j=1 Ej f Y |X (y|x (j) )dy . (35)
where, P C is the type defined in (18); the pmf Q is defined in (23); for all j ∈ {1, 2, . . . , L}, the sets E j are defined in (25); and the expectation in (34) is with respect to W , which is a complex circularly symmetric Gaussian random variable whose real and imaginary parts have zero means and variances 1 2 σ 2 .

Proof. The proof of Theorem 1 follows from Corollary 4, Lemma 1; and Corollary 2.

Theorem 1 quantifies the trade-off between all parameters of a homogeneous (n, M, , B, δ)-code for the random transformation in [START_REF] Tesla | Apparatus for transmitting electrical energy[END_REF]. Let such a code be denoted by C of the form in [START_REF] Huang | Simultaneous wireless information and power transfer: Technologies, applications, and research challenges[END_REF]. Hence, both the upper bound on the information rate in (33) and the upper bound on the energy rate in (34) depend on the type P C . The information rate is essentially upper bounded by the entropy of P C . Hence, codes whose codewords are such that every channel input symbol is used the same number of times are less constrained in terms of information rate. This is the case in which P C is a uniform distribution. Alternatively, using a uniform type P C might dramatically constrain the energy transmission rate. For instance, if the constellation is such that for at least one pair

(x 1 , x 2 ) ∈ X 2 , it holds that E [g(x 1 + W )] < E [g(x 2 + W )]
where the random variable W is defined in Lemma 1, then using the symbol x 1 equally often as x 2 certainly constraints the energy rate. Codes than might potentially exhibit the largest energy rates are those in which the symbols x that maximize E [g(x + W )] are used more often. Clearly, this deviates from the uniform distribution and thus, constraints the information rate.

Another interesting trade-off appears between the DEP and EOP. From Corollary 1, it follows that codes whose codewords contain mainly channel input symbols x that maximize E [g(x + W )] are more reliable from the perspective of energy transmission. Alternatively, from Corollary 3, it follows that codes that are more reliable in terms of information transmission are those whose codewords contain more channel input symbols with the smallest DEP. That is, when P C approaches the pmf Q in (23).

IV. EXAMPLES

Consider a homogeneous (n, M, , B, δ)-code C (Definition 4) for the random transformation in (1) of the form in ( 8) built upon a 16-QAM constellation. More specifically, the set of channel input symbols X in ( 5) is: X = {(5 + i5), (-5 + i5), (-5 -i5), (5 -i5), (15 + i5), (5 + i15), (-5 + i15), (-15 + i5), (-15 -i5), (-5 -i15), (5 -i15), (15 -i5), (15 + i15), (-15 + i15), (-15 -i15), (15 -i15)}. The decoding sets D 1 , D 2 , . . ., D 16 in (27) are those obtained from the maximum-à-posteriori (MAP) estimator [START_REF] Proakis | Digital Communications[END_REF] given the type P C in (18). The block-length is fixed at n = 80 channel uses and the noise variance σ 2 = 2 in (4). The energy harvested at the EH is calculated taking into account the non-linearities at the receiver described in [START_REF] Varasteh | Wireless information and power transfer over an AWGN channel: Nonlinearity and asymmetric Gaussian signaling[END_REF]. In particular, given a type P C , the expression In Figure 1, given a type P C , the upper bound on the information transmission rate and the upper bound on the energy transmission rate in (33) and (34), respectively, are depicted. Each point in the curves corresponds to a different type P C . Therein, it becomes evident that smaller EOPs imply smaller information and energy transmission rates. The maximum information rate achievable with a 16-QAM constellation is four bits per channel use. Nonetheless, when codewords are forced to exhibit the same type (homogeneous code), the information transmission rate does not exceed 3.51 bits per channel use. This maximum holds when the energy transmission rate is smaller than 2 × 10 4 energy units, which corresponds to the maximum energy that might be delivered by codes whose type P C approaches the uniform distribution on X . Beyond this value, increasing the energy transmission rate implies deviating from the uniform distribution, which penalises the information transmission rate. The maximum energy rate implies types P C that concentrate on the symbols (15+i15), (-15+i15), (-15-i15), (15-i15). A type P C that approaches the uniform distribution over these four channel inputs induces a lower bound on the information transmission rate of 1.9 bits per channel use. A type P C that concentrates on one of these channel input symbols induces a zero information rate.

L =1 P C (x ( ) )E g x ( ) + W in (34)

Fig. 1 :

 1 Fig.1: Upper bounds on the information rate in (33) and energy rate in (34) for a homogeneous code built upon a 16-QAM constellation as a function of the EOP δ.

Fig. 2 :

 2 Fig. 2: Upper bound on the energy rate in (34) for a homogeneous code built upon a 16-QAM constellation as a function of the lower bound on the DEP in (35).

Figure 2

 2 Figure2illustrates the upper bound on the energy transmission rate in (34) as a function of the lower bound on the DEP in (35) when the EOP is kept constant at δ = 10 -4 and n = 80. More specifically, each point on the curve in Figure2for a given value of the DEP is associated with a type P C . This type is used to determine the upper bound on the energy transmission rate, which leads to the curve in Figure2. Higher DEPs require types P C with less constraints, which increases the bound on the energy rate.

This research was supported in part by the European Commission through the H2020-MSCA-RISE-2019 program under grant 872172; in part by the Agence Nationale de la Recherche (ANR) through the project MAESTRO-5G (ANR-18-CE25-0012); and in part by the Fondation Mathématique Jacques Hadamard (FMJH) through the Programme Gaspard Monge.

Simultaneous Information and Energy Transmission with Finite Constellations


Rola El Osta
email: rola.elosta@ul.edu.lb

Maryline Chetto
email: maryline.chetto@univ-nantes.fr

Online Admission of Hard Aperiodic Tasks in Real-Time Energy Harvesting Systems

Keywords: Real-time system, Energy harvesting, Schedulability analysis, Online scheduling, Aperiodic task, Time laxity, Energy laxity

It is expected that energy harvesting technology will increasingly be deployed to realize autonomous embedded realtime systems such as wireless sensors. In this paper, we consider the problem of jointly scheduling periodic tasks and dynamically arriving hard deadline aperiodic tasks in a monoprocessor system powered with regenerative energy. Our result is an exact admission algorithm in that sense that it accepts every schedulable task based on the amount of available energy and computing capacity. It is assumed that all the accepted tasks are preemptively scheduled according to the optimal algorithm, ED-H which is an energy aware extension of Earliest Deadline First. We also discuss issues concerning the implementation of the admission test.

I. INTRODUCTION

Nowadays, batteries are the dominant source of energy for small electronic devices such as wireless sensor nodes. However, besides their impact on the environment, their limited energy storage capacity and their finite lifetime make their use restricted. For these reasons, alternative energy sources have been sought. Energy harvesting (EH), or energy scavenging, is a technology that captures unused ambient energy and convert it into usable electrical energy that can be used immediately or later thanks to a storage unit [START_REF] Tanaka | Energy Harvesting Technology for Maintenance-free Sensors[END_REF]. Generally, we use this terminology for low power and small autonomous devices such as wireless sensor networks and portable electronic equipments. Hence, ambient energy scavenging has engendered a lot of interest for researchers. It permits that autonomous embedded systems be supplied perpetually. In comparison with energy stored in classical storage unit as batteries, the environment represents an infinite source of available energy. Many environmental sources can be exploited, including solar energy, electromagnetic waves, thermal energy, mechanical, etc. and must be selected based on the application characteristics. An EH system consists of three components (see Fig. 1): a single processing unit, an energy harvester and a rechargeable energy storage. The energy harvester scavenges energy from the environment and converts it to usable electrical power. The rechargeable energy storage (e.g battery or capacitor) stores energy. And after that, the processor unit uses this energy to execute the programs. In EH systems, we have to exploit perpetually the available ambient energy which is strongly dependent on the environment. The consumed energy of the system should be adapted to maximize the performance instead of minimizing it as in the classical battery powered systems. So, the operating system has to manage properly the activity of the processing unit such that there is sufficient energy in the storage unit to satisfy all the constraints at every time. Most environmental energy sources like the electromagnetic waves, do not deliver a constant power over time. The energy generated may arrive in bursts and has to be stored so that the device can still be operated at a later moment. The main challenge for an EH system lies in the optimization of its performance while respecting the time-varying amount of energy without wasting or exhausting the energy stored in the battery or capacitor. We then say that the system operates in an energy neutral mode by consuming only as much energy as harvested [START_REF] Kansal | Harvesting aware power management for sensor networks[END_REF].

Most of real-time applications consist of a set of hard periodic tasks. Hard aperiodic tasks can arrive at a certain time from alert conditions or from failures of hard periodic tasks which fail to check and validate their results and must be retried and completed before the original deadline [START_REF] Thuel | Enhancing Fault Tolerance of Real-Time Systems through Time Redundancy[END_REF]. Unlike the problem of scheduling periodic tasks only, the joint scheduling problem is difficult because the scheduler has to choose which tasks to accept for processing when not all of the timing requirements can be met due to processing overload or energy shortage. The operating system should be provided with a procedure called acceptance test which dynamically decides whether a newly arriving aperiodic task could be accepted [START_REF] Tia | A linear-time optimal acceptance test for scheduling of hard real-time tasks[END_REF]. As aperiodic tasks are not known before run-time, the acceptance test is performed online and at unpredictable times so as to guarantee the deadlines of all periodic tasks and any already accepted. If a hard aperiodic task cannot be guaranteed, it is rejected. The mechanism of the acceptance test was also presented in [START_REF] Kim | Scheduling of Hard-Aperiodic Requests in Dynamic Priority Systems[END_REF] based on the the APS (Alternative Priority Algorithm). In this paper, we specifically focus on hard deadline periodic tasks with hard deadline aperiodic tasks that execute on a uniprocessor platform which is powered by regenerative energy. We present the theoretical basis of our approach to solve the acceptance problem. This method is based on the assumption that all periodic and aperiodic tasks are scheduled according to the optimal algorithm ED-H [START_REF] Chetto | Optimal Scheduling for Real-Time Jobs in Energy Harvesting Computing Systems[END_REF] which is an energy aware extension of Earliest Deadline First [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF]. We will answer to the acceptance question by calculating, for each newly arriving aperiodic task, two key values called time laxity and energy laxity. The time laxity is a measure of processing time surplus. The energy laxity is a measure of energy surplus. Both will be assessed and available after accepting the task. Clearly, these two values should be non negative to admit the newly arriving task. The remainder of the paper is organized as follows. Section II presents background materials about Earliest Deadline First scheduling and online admission of aperiodic tasks with no energy considerations. Section III introduces the model and terminology relative to the scheduling issue in real time energy harvesting systems.Section IV describes the ED-H scheduler. Our admission test is presented in section V. In Section VI, we give some suggestions for an efficient implementation of the test. Section VII summarizes the paper and describes future work.

II. BACKGROUND MATERIAL

Hereafter, we give necessary background about scheduling real-time tasks considering that energy is unlimited. We only focus on preemptive scheduling for deadline constrained and independent tasks.

A. Scheduling periodic tasks with EDF

The most common dynamic priority scheduling algorithm for real-time systems is the Earliest Deadline First (EDF) which was introduced by Liu and Layland [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF] in 1973. EDF is clearly a dynamic priority driven scheduler since closer is the deadline, higher is the priority. The EDF algorithm has been proven by Dertouzos to be optimal among all scheduling algorithms on a uniprocessor in the sense that if a task set cannot be feasibly scheduled by EDF, then this task set cannot be feasibly scheduled by any other algorithm [START_REF] Dertouzos | Control Robotics: The Procedural Control of Physical Processes[END_REF]. The classical implementation of EDF consists in executing the tasks under the ASAP (As Soon As Possible) version. We then say that EDF is greedy, consuming the available energy whenever at least one task is pending for execution and ignoring the future energy requirements. Let us notice that the ALAP (As Late As Possible) version of EDF is no more appropriate to schedule deadline constrained tasks in energy harvesting systems.

B. Admission test for aperiodic tasks and EDF

Most of real-time software consist of a set of hard periodic tasks with aperiodic tasks that may be authorized for execution as long as they do not cause any periodic task to violate its deadline. A hard aperiodic task consists of one computation that responds to an internal or external event. The admission test considers the processing capacity required by future periodic invocations and aperiodic tasks. The optimal admission test is based on the exact computation of the largest amount of processor idle time available during the interval between the arrival time and the deadline to complete the processing, while still ensuring that all accepted tasks meet their deadline. The fact that the utilization of the periodic tasks is less than 1 means that over the hyperperiod given by the least common multiple of the periods they are a positive number of slack time units remaining to be used by aperiodic tasks. However, to guarantee the deadlines of the periodic tasks while executing the aperiodic tasks as soon as possible before deadlines, this slack should be distributed adequately. In that objective, we have to exploit the central property of the earliest-deadlinelate (EDL) algorithm which schedules periodic tasks as late as possible. Whenever a hard aperiodic task occurs, the idle times of the EDL schedule are computed and compared to the execution requirements of the aperiodic tasks not yet completed, checking that the aperiodic tasks will not use more idle time than available. The optimality of this test was stated in [START_REF] Chetto | Scheduling Periodic and Sporadic Tasks in a Real-Time System[END_REF]. Nonetheless, such a test considers that there is no limitation on energy availability and should be revisited to take into account, in one part the energy required by any task in execution and in the other part both the profile of the energy source and the size of the energy storage unit.

III. MODEL UNDER STUDY

The Energy Harvesting model (hereafter, referred to as RTEH) consists of a computing element, a set of tasks, an energy storage unit, an energy harvesting unit and an energy source. We consider a set of real-time tasks that are executed on a uniprocessor system that supports only one operating frequency and that consumes negligible energy in the idle state when it does not execute any task. The system contains a set of n independent periodic tasks. We refer to this set by τ = {τ i (C i , D i , T i , E i ); i = 1, ..., n}, where task τ i has a worst case execution time of C i time units, a relative deadline D i , a period T i (with D i ≤ T i ) and a Worst Case Energy Consumption of E i energy units. We assume that E i is not necessarily proportional to C i [START_REF] Jayaseelan | Estimating the Worst-Case Energy Consumption of Embedded Software[END_REF]. Each task gives rise to an infinite sequence of invocation requests that starts at the time origin, here equal to zero. Let H, the hyperperiod, be equal to the least common multiple of the periods T 1 , T 2 ,. . . ,T n .

Let t c denote the current time. We next introduce the aperiodic tasks. It is assumed that there is no a priori knowledge about which set of aperiodic task requests will be encountered. Furthermore, it is supposed that any aperiodic task is ready to execute as soon as it arrives and does not interact with periodic tasks in the sense that it does not share resources. Let t c the current time that coincides with the arrival of a new aperiodic task. Let J = {J i |1 ≤ i ≤ m} the stream of m hard aperiodic tasks already admitted at t c . Every uncompleted aperiodic task at t c is characterized by c i , e i and d i that respectively denote its remaining execution time, remaining energy consumption and absolute deadline. In what follows, J c (c, e, d) denotes the arriving aperiodic task to be tested.

Let P p (t) be the instantaneous charging rate produced by the power source that incorporates all losses. The energy produced on [t 1 , t 2 ) is given as E p (t 1 , t 2 ) = t2 t1 P p (t)dt. We assume that the energy produced in any unit of time can overlap with the energy consumed in one unit of time. We also assume that a short term prediction of the energy harvesting rate is possible. Our system uses an ideal energy storage unit with a nominal capacity C of units of energy. The energy level at time t is denoted E(t). The stored energy may be used at any time later and does not leak any energy over time. Finally, we also suppose that the available average power from the energy source is sufficient to support the long term power demand of the periodic task set τ .

IV. THE ENERGY AWARE ED-H SCHEDULING ALGORITHM

A. Principles

In [START_REF] Chetto | Optimal Scheduling for Real-Time Jobs in Energy Harvesting Computing Systems[END_REF], the author has introduced a novel energy-aware scheduling algorithm, namely ED-H. ED-H is an energy harvesting aware variant of EDF, proved to be optimal for real time tasks scheduling under energy harvesting settings. ED-H answers the question: when should the processing unit use energy to execute a deadline constrained task, and when should the processing unit idle and recharge the energy storage unit. ED-H allows a task to execute only if its execution cannot cause energy starvation for a future occurring task. In that sense, ED-H is a clairvoyant algorithm since its decision necessarily lies on the short term prediction of the incoming energy and the pattern of future task arrivals.

B. Implementation

The implementation of ED-H needs the online computation of the following values:

• The slack time of the task set τ at current time t, denoted ST τ (t) and defined as the maximum continuous processor idle time that could be made available from time t while still guaranteeing the feasibility of all the tasks in the set τ . Whenever the energy storage unit is completely deplenished, the system enters the idle state so as to recharge the energy storage during at most ST τ (t) unit of time.

• The slack energy of the task set τ at current time t, denoted SE τ (t) and defined as the maximum energy that could be made available and consumed from time t while still guaranteeing the feasibility of all the task set τ . Clearly, no energy should be wasted so as to guarantee that SE τ (t) be always non negative. • The preemption slack energy of the task set τ at current time t, denoted P SE τ (t) and defined as the maximum energy that could be consumed from time t by the current active task while still guaranteeing the feasibility of all the future tasks with a higher priority i.e. the tasks that may preempt the current active one.

C. Illustrative Example

Consider a periodic task set τ that is composed of three periodic tasks. [START_REF] Kim | Scheduling of Hard-Aperiodic Requests in Dynamic Priority Systems[END_REF][START_REF] Chetto | Optimal Scheduling for Real-Time Jobs in Energy Harvesting Computing Systems[END_REF]12), τ 2 = (2, 8, 10, 15) and τ 3 = (4, 11, 15, 22). We assume that the energy storage capacity is E = 40. For simplicity, the rechargeable power, P p , is constant along time and equals 5. First of all, we have to schedule the periodic task set τ according to ED-H. We verify that τ is schedulable since all the tasks are executed before their deadline and without depleting the energy reservoir (Fig. 2). At time t = 0, all tasks are released. τ 1 is the highest priority task and is executed until t = 1 where E(1) = 33 energy units. At time t = 1, τ 2 is the highest priority task and is executed until t = 3 where E(3) = 28 energy units. τ 3 is now the highest priority task and is executed until t = 7 where the energy storage capacity becomes 26 energy units. At time t = 7, τ 1 as the highest priority task, ready to be processed, runs and finishes at t = 8. E(8) = 19 energy units. From time t = 8 up to t = 10, the processor remains idle because there are no pending tasks. During that time interval, the energy storage will recharge and the energy level at t = 10 is given by E(10) = 29 energy units. At t = 10, τ 2 the highest priority task, is ready to be processed. Here, computation of the slack energy is necessary because task τ 1 with deadline less than deadline of τ 2 that will be released after time t = 10. SE τ1 (t = 10) = E [START_REF] Jayaseelan | Estimating the Worst-Case Energy Consumption of Embedded Software[END_REF]

τ = {τ i | 1 ≤ i ≤ 3 and τ i = (C i , D i , T i , E i )}. Let τ 1 = (1,
+ 17 10 P p dt -E 1 -E 2 = 37 > 0.
As E(10) = 29 and as the slack energy is greater than zero, τ 2 is authorized to execute immediately. At time t = 12, E(12) = 24. Now, τ 1 has the highest priority. It is executed till t = 13, where E(13) = 17 energy units. Again, the processor is idle from t = 13 up to t = 15 where the energy reservoir will recharge leading to E(15) = 27 energy units. At t = 15, τ 3 is ready and has the highest priority. Slack energy is requires to be computed:

SE τ3 (t = 15) = E(15) + 23 15 P p dt -E 1 -E 3 = 33 > 0.
Consequently, the slack energy is 33 energy units. τ 3 is executed till t = 18, where preemption occurs since τ 1 has a higher priority. At t = 18, we can evaluate the maximum energy that remains to be consumed by τ 3 as follows: we memorize the energy level of the storage Fig. 2. ED-H scheduling on a periodic task set unit at every scheduling point i.e. whenever a preemption occurs or a task begins execution. Thus, we deduce that:

E rem 3 = E 3 -(E(15) -E(18)) = 15.
At t = 18, τ 1 is the highest priority task ready to be processed, runs and finishes at t = 19. E(19) = 13 energy units. τ 3 will now resume its execution until t = 20. At t = 20, τ 2 is the highest priority task ready to be processed and E(20) = 3. But there is no sufficient energy in the energy storage unit for execution. So, we have to insert an idle time to let the processor inactive as long as the energy storage has not filled completely and the latest start time of the next periodic task has not been attained. The slack time is equal to 6. Hence, the processor is let idle till t = 26 permitting the energy reservoir to recharge so as E(26) = 33 energy units. τ 2 can be executed until t = 28, where E(28) = 28 energy units. At t = 28, τ 1 is the highest priority task ready to be processed. It executes and finishes at t = 29 with E(29) = 21 energy units. All invocations of periodic tasks are feasibly scheduled until the end of the hyperperiod where the energy reservoir contains 26 energy units.

V. ACCEPTANCE TEST

This section is concerned with the dynamic acceptance test that answers the following decision question relative to a given arriving aperiodic task J c : "Can J c be accepted?". The approach is based on the assumption that all periodic and aperiodic tasks are executed according to the optimal scheduling algorithm ED-H. The main principle of the admission test is to authorize aperiodic job executions as long as it does not involve a deadline violation for all the jobs generated by the periodic task set τ and the previously accepted aperiodic tasks. Let us recall that a deadline violation occurs either because of processing time starvation (lack of time to complete a task before deadline) or energy starvation (lack of energy to complete a task before deadline). The basic idea of our acceptance test is to determine both processing time and energy that could be stolen for the newly arriving task. In order to provide an exact admission test, we will separate the constraints in timing and energy domains.

A. Time schedulability analysis

Let J c (c, e, d) be the arriving aperiodic task. The acceptance test of J c is executed once at the time of its arrival, t c . The idea is to postpone the execution of periodic activities, making any spare processing time available as soon as possible for the aperiodic activities. Our approach is based on the exact computation of the largest amount of processor idle time available during the interval between the arrival time and the deadline to complete the processing, while still ensuring that all guaranteed tasks meet their deadline. Consequently, when J c occurs, the idle times of the EDL schedule are computed and compared to the execution requirements of the aperiodic tasks not yet completed [START_REF] Chetto | Scheduling Periodic and Sporadic Tasks in a Real-Time System[END_REF].

Let us introduce δ t i (t c ) called the dynamic time laxity of task J i at time t c . By definition, δ t i (t c ) gives the highest quantity of processing time that could be consumed between t c and d i by aperiodic tasks with priority higher than or equal to J i while guaranteeing the deadlines of all the tasks. Let Ω EDF τ ( tc) (t c , d i ) be the processing time which may be stolen from the periodic task set between t c and d i . As proved in [START_REF] Chetto | Scheduling Periodic and Sporadic Tasks in a Real-Time System[END_REF] it is derived from the EDL schedule produced at current time t c . Consequently, we have

δ t i (t c ) = Ω EDF τ (t c , d i ) -Σ i j=1 c j .
We may draw the following necessary condition for the admission test [START_REF] Osta | An exact admission test for hard deadline aperiodic tasks with energy harvesting considerations[END_REF]:

Lemma 1: Aperiodic task J c is accepted only if for every aperiodic task J i such that d i ≥ d, δ t i (t c ) ≥ 0 (1)
Note that Σ i j=1 c j gives the total remaining execution time of all the aperiodic tasks with a priority higher than J i .

The value of Ω EDF τ (t c , d i ) is dynamically computed from two arrays which respectively give the start times and the durations of all idle time intervals in the EDL schedule starting at current time t c and finishing at the end of the current hyperperiod.

B. Energy schedulability analysis

Let us introduce δ e i (t c ) the energy laxity of task J i at time t c . δ e i (t c ) gives the maximal energy which may be stolen from the periodic task set between t c and d i and consumed by aperiodic tasks with priority higher than or equal to J i . Let Π EDF τ ( tc) (t c , d i ) give the maximal energy which remain after executing the periodic tasks between t c and d i . Consequently,

δ e i (t c ) = Π EDF τ (t c , d i ) -Σ i j=1 e j . Π EDF τ ( tc) (t c , d i ) is computed through the following equation: Π EDF τ (t c , d i ) = E(t c ) + E s (t c , d i ) -g(t c , d i ) (2)
where

• E(t c
) is the residual capacity of the energy reservoir at t c ,

• E s (t c , d i ) is the total energy produced by the environmental source between t c and d i ,

• g(t c , d i ) = n j=1 di-tc Tj E j + n k=1(d k >di) E k (t)
and refers to the energy demand issued from the periodic tasks between t c and d i . Here, E k (t) denotes the remaining amount of energy of the current invocation of the periodic task τ k at t c . We may draw the following necessary condition for the admission test [START_REF] Osta | An exact admission test for hard deadline aperiodic tasks with energy harvesting considerations[END_REF]:

Lemma 2: Aperiodic task J c is accepted only if for every aperiodic task J i such that d i ≥ d, δ e i (t c ) ≥ 0 (3)

C. Exact admission test

From the previous lemmas, we prove the following exact admission test [START_REF] Osta | An exact admission test for hard deadline aperiodic tasks with energy harvesting considerations[END_REF]:

Theorem 1: Aperiodic task J c is accepted if and only if for every aperiodic task J i such that

d i ≥ d, δ t i (t c ) ≥ 0 and δ e i (t c ) ≥ 0 (4)
The purpose of this admission test is to accept or reject any arriving task based on first, the amount of processing time and energy required for its execution and second, the amount of processing time and energy available. Theorem 1 tells us that energy and time constraints can be decoupled if we have to decide whether a newly arriving task is admitted or not.

D. Illustrative example

Consider the previous periodic task set τ and suppose now that some aperiodic tasks arrive to be jointly executed with τ . Let us consider J = {J j | J j = (a i , c i , d j , e j )} the stream of 3 aperiodic tasks where J 1 = [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF][START_REF] Tia | A linear-time optimal acceptance test for scheduling of hard real-time tasks[END_REF][START_REF] Osta | An exact admission test for hard deadline aperiodic tasks with energy harvesting considerations[END_REF]17), J 2 = (18, 3, 26, 20) and J 3 = (20, 2, 28, 10). Periodic task set τ is scheduled according to ED-H till t = 7 where E(7) = 26 energy units (Fig. 3). At time 7, J 1 is released. The admission test requires to compute energy laxity and time laxity of J 1 at t = 7. δ e 1 (7) = E(7) +

11 7 P p dt -E 1 -E 2 -e 1 = 7 > 0 and δ t 1 (7) = 3 -4 = -1 < 0.
As the slack time is less than zero, from lemma 1, J 1 is rejected and not authorized for execution.

J 2 arrives at t = 18 (Fig. 4). J 2 is rejected since δ e 2 (18) = -2 < 0.

At t = 20, J 3 arrives (Fig. 5). δ e 3 (20) = 6 and δ t 3 (20) = 4. As both the time laxity and energy laxity of J 3 are positive, the test leads to admit J 3 . As the absolute deadline of J 3 is 28, J 3 is the highest priority task in the ready list. At t = 20, we have E(t) = 3 and P SE(t) = 18. J 3 may execute immediately to be finished at t = 22. τ 2 is the highest priority task ready for execution at t = 22. Nonetheless, there is no sufficient energy in the energy storage unit for execution since E(t) = 3. The energy storage unit should be recharged. The slack time has to be computed in computation of its energy laxity requires to estimate future incoming energy. The better the approximation, the better the acceptance algorithm performs in terms of exactness. During short time intervals, the produced power can be considered as constant. Let us consider an energy harvester composed of a solar panel in order to draw energy from outdoor light. A sudden change of the light intensity is improbable during a time interval with a duration in the order of the second.

Let us reminder that a real time periodic task generally consists in sampling so as to control a physical process. As a consequence, we may consider that the incoming power is constant over time. This simplifies online computations since sophisticated prediction methods may not be necessary.

There are mainly two techniques to store energy in a small device such as a wireless sensor node. The most common one uses rechargeable batteries. The second technique uses supercapacitors or ultra-capacitors which support higher power flows than batteries. In addition, they offer very long lifetimes with very high number of charge/recharge cycles without performance degradation. An important feature of energy harvesting systems is the capability to estimate the remaining energy in the storage unit. Unlike batteries, the energy of supercapacitors can be measured in a straightforward way. So, we may assume that at every time instant, the operating system may determine precisely the energy level in the storage unit and the energy harvested on a short term time interval in order to compute the energy laxity with accuracy according to equation 3. It follows that the implementation complexity of the acceptance test is O(k.n + m) in the worst case where n is the number of hard periodic tasks, m is the number of previously guaranteed but yet uncompleted hard aperiodic tasks and k is an integer which depends on periods and deadlines.

If the technique used to store energy considers a non-ideal energy storage unit, consuming the stored energy as soon as possible is recommended to avoid energy leakage and task execution failure. However, when the leakage rate of the ultracapacitor is very high during charging and discharging, we may calculate the amount of energy leakage of ultra-capacitor for each periodic task and each arriving aperiodic task at charging and discharging phases. We take into account energy leakage by modifying the formula for the computation of the energy laxity.

VII. CONCLUSION AND FUTURE WORKS

This paper has presented an exact admission test to prevent processing overload and energy shortage in energy harvesting systems such as sensor nodes in the IoT. Such systems consist of a baseload of hard deadline periodic tasks where aperiodic tasks may be accepted for execution as long as they do not cause any periodic task to miss its deadline. We assume that all the deadline constrained tasks are scheduled according to the optimal scheduling algorithm ED-H [START_REF] Chetto | Optimal Scheduling for Real-Time Jobs in Energy Harvesting Computing Systems[END_REF]. Upon its arrival, any aperiodic task is tested for admission calculating two key values: the time laxity and the energy laxity. Provided that precise information about the future energy generation is available, the proposed online admission test is exact: it enables us to continuously optimize the capacity of the system in the two dimensions: processing time and energy. Work is in progress to extend our analysis to more general applicative software composed of deadline constrained tasks that may share resources and exhibit precedence constraints.

Fig. 1 .

 1 Fig. 1. A typical Real-Time Energy Harvesting System

Fig. 5 .

 5 Fig. 5. Admission of J3; δ e 3 (20) = 6 > 0 and δ t 3 (20) = 4 > 0


David Simoncini Lisbp

Sophie Barbe
email: sbarbe@insa-toulouse.fr

Thomas Schiex
email: tschiex@toulouse.inra.fr

Sébastien Verel
email: verel@uni-littoral.fr

Fitness Landscape Analysis around the Optimum in Computational Protein Design

The geometry and properties of the fitness landscapes of Computational Protein Design (CPD) are not well understood, due to the difficulty for sampling methods to access the NP-hard optima and explore their neighborhoods. In this paper, we use a state-of-theart AI complete algorithm to enumerate all solutions within a 2 kcal/mol energy interval of the optimum of two CPD problems. We compute the number of local minima, the size of the attraction basins, and the local optima network. We provide various features in order to characterize the fitness landscapes, in particular the multimodality, and the ruggedness of the fitness landscape. Results show some key differences in the fitness landscapes and help to understand the successes and failures of metaheuristics on CPD problems. Our analysis gives some previously inaccessible and valuable information on the problem structure related to the optima of the CPD instances (multi-funnel structure), and could lead to the development of more efficient metaheuristic methods.

INTRODUCTION

Present in all living organisms, proteins are polymeric chains of amino acids that play a central role in cellular processes such as gene expression, catalysis, communication, regulation, transport and signaling. The succession of amino acids in a protein sequence defines how the protein folds into a given three-dimensional (3D) structure, and thus its biological function. By changing the aminoacid sequence, protein design seeks to produce proteins with new structures and functions. Besides the ever-growing interest to produce tailor-made proteins for applications for the industry and medicine, protein design is also motivated by the desire to understand the relationships between the sequence, evolution, structure and function of biomolecules.

Experimental protein design methods consists in synthesizing various amino-acid sequences using genetic engineering and testing them. However, with a choice among 20 natural amino acids at each position, the size of the combinatorial sequence space (20 𝑛 sequences for a protein of length 𝑛) is out of reach of current experimental methods. Therefore, computational methods have been proposed for rationalizing protein design process by focusing experimental tests on a computationally selected small library of sequences of particular relevance for the targeted structure/function.

To achieve this, Computational Protein Design (CPD) tries to identify, in this exponential amount of amino acid sequences, those that are compatible with a targeted 3D structure, chosen for its functional and/or structural properties. To achieve this, CPD relies on an energy function to assess the stability of the amino-acid sequence in the desired structure and a search algorithm to identify a sequence that can organize itself in space in such a way that it reaches optimal stability. An optimally organized sequence defines a Global Minimum Energy Conformation or GMEC.

Any amino acid can be decomposed in a central invariable part that participates in the formation of the relatively rigid polymeric chain (or backbone) and a variable and highly flexible lateral chemical component that defines the amino acid and that can be chosen among the 20 possible natural variants. For these reasons, and following an approach introduced by Ponder and Richards [START_REF] Jay | Tertiary templates for proteins: use of packing criteria in the enumeration of allowed sequences for different structural classes[END_REF], most CPD approaches assume a rigid backbone structure while amino acid side-chains are allowed to move among a finite set of compositions and preferred organization conformations, called rotamers. While CPD is still a young and rapidly evolving field, success stories of computationally designed proteins highlight its ability to adequately capture fundamental rudiments of molecular recognition and interactions enabling the design of several kinds of proteins for different purposes [START_REF] Sagar D Khare | Computational redesign of a mononuclear zinc metalloenzyme for organophosphate hydrolysis[END_REF][START_REF] Verges | Computer-Aided Engineering of a Transglycosylase for the Glucosylation of an Unnatural Disaccharide of Relevance for Bacterial Antigen Synthesis[END_REF][START_REF] Arnout | Computational design of a self-assembling symmetrical beta-propeller protein[END_REF].

The CPD problem can be formulated as an optimization problem which consists in searching for combinations of rotamers at designable positions that will lead to a stable 3D-fold. Despite its apparent simplicity, the rigid backbone/discrete rotamer CPD problem has been proven NP-hard [START_REF] Niles | Protein design is NP-hard[END_REF] and hard to approximate [START_REF] Chazelle | A semidefinite programming approach to side chain positioning with new rounding strategies[END_REF]. Consequently, most CPD approaches rely on metaheuristic search methods based on simulated annealing [START_REF] Kuhlman | Native protein sequences are close to optimal for their structures[END_REF][START_REF] Voigt | Trading accuracy for speed: A quantitative comparison of search algorithms in protein sequence design[END_REF] (implemented into the Rosetta modeling software [START_REF] Leaver-Fay | ROSETTA3: an object-oriented software suite for the simulation and design of macromolecules[END_REF]), genetic algorithms [START_REF] Raha | Prediction of amino acid sequence from structure[END_REF] or other local search algorithms [START_REF] Desmet | Fast and accurate side-chain topology and energy refinement (FASTER) as a new method for protein structure optimization[END_REF][START_REF] Kyle Lassila | Combinatorial methods for small-molecule placement in computational enzyme design[END_REF]. Although metaheuristics have the advantage of providing a solution at any time, they neither guarantee finding the global optimal solution (or GMEC) in finite time nor a bounded energetic distance to the optimal solution. To try to circumvent this limitation, multiple independent runs are performed (each with a predefined number of steps) in order to cover, as well as possible, a rugged energy landscape. However, the accuracy of metaheuristic methods drastically degrades as problem size increases [START_REF] Simoncini | Guaranteed Discrete Energy Optimization on Large Protein Design Problems[END_REF][START_REF] Voigt | Trading accuracy for speed: A quantitative comparison of search algorithms in protein sequence design[END_REF] and the probability of finding the GMEC drops very quickly close to 0 as problems get more difficult. Additionally, the average energy gap to optimality tends to increase with the number of designable positions, putting a limit on the size of systems for which a reasonably good solution can be found with confidence. The lack of knowledge on the features of the energy landscape of CPD problems makes it difficult to understand why metaheuristics methods fail on given instances.

Alternative CPD approaches rely on exact deterministic algorithms guaranteeing that the solution returned at the completion of the algorithm is the GMEC. Unfortunately, these methods, historically based on the Dead-End Elimination theorem combined with the A* algorithm [START_REF] Leach | Exploring the conformational space of protein side chains using dead-end elimination and the A* algorithm[END_REF] are often rapidly outstripped by the complexity of the search space and do not provide any solution in reasonable time. Recent artificial intelligence methods have allowed to push back these limitations [START_REF] Allouche | Computational protein design as an optimization problem[END_REF][START_REF] Simoncini | Guaranteed Discrete Energy Optimization on Large Protein Design Problems[END_REF][START_REF] Traoré | A New Framework for Computational Protein Design through Cost Function Network Optimization[END_REF][START_REF] Traoré | Fast search algorithms for computational protein design[END_REF]. Based on graphical models and more specifically Cost Function Networks (CFN) [START_REF] Cooper | Soft arc consistency revisited[END_REF] methods, they can handle complex search spaces that were previously unsolvable by state-of-the-art provable CPD methods. The CFN-based approaches speed-up search by several orders of magnitude and can usually provide a guaranteed GMEC in reasonable time. In addition to the optimal solution, these methods can provide an exhaustive list of sub-optimal solutions in a given energy threshold around the optimum. This new ability opens new opportunities, notably to exhaustively explore the energy landscape in the neighborhood of global optima. In this paper, we use CFN-based methods to access the guaranteed NP-hard optima of CPD problems and exhaustively explore their neighborhoods. This study provides previously inaccessible and highly valuable data on the CPD energy landscape. Analysis highlights some key characteristics of the two CPD problems explaining the distinct behaviors of a metaheuristic, the simulated annealing algorithm implemented in Rosetta.

RESULTS

For each of two CPD instances, 2ckx and 2gkt, using the Toulbar2 CFN solver, the optimal solution (GMEC), was identified and an exhaustive enumeration of sub-optimal solutions within an energy interval of 2 𝑘𝑐𝑎𝑙 .𝑚𝑜𝑙 -1 above the GMEC energy was performed. The number of solutions within this energy interval is 497, 282 for 2gkt and 6, 209, 729 for 2ckx. These two CPD instances were selected because of the difference of performance achieved with simulated annealing. In a previous study, a simulated annealing algorithm was able to reach the GMEC 250 times out of 1000 trajectories on the instance 2gkt and failed to reach the GMEC on 2ckx, even after 1 million trajectories [START_REF] Simoncini | Guaranteed Discrete Energy Optimization on Large Protein Design Problems[END_REF]. These results were obtained using the Rosetta software with talaris14 energy function. In order to confirm this difference in performance, we ran 1, 000 trajectories of simulated annealing on both instances using Rosetta software with the recent beta_nov16 energy function (see details in Section 3). Our result confirms that 2gkt has several features indicating that it is easier than 2ckx. For 2gkt, the simulated annealing was able to find the GMEC 24 times, and found 98 distinct sequences in the enumerated landscape (269 trajectories found a solution within 2kcal/mol of the GMEC, several trajectories finding the same sequence). In the case of 2ckx, the simulated annealing never found the GMEC and only 13 sequences (out of the 6, 209, 729 possibilities) appeared in the enumerated landscape. In this section, we present a fitness landscape analysis of the two CPD problems.

For a given solution, its fitness is the difference of energy of the optimal solution and this solution (maximum at zero).

Fitness distance correlation

The density of states according to the fitness value decreases exponentially for both instances, with a value of 𝑅 2 of the regression over 0.98 (Figure 1). The density of states decreases faster for 2ckx, which means fewer solutions near the optimum and suggests that this instance is more difficult. The distributions of the Hamming distances to the optimal sequence have different profiles (Figure 2). The distribution is unimodal for 2gkt, with a mean distance to the optimal sequence of 6. The distribution is bimodal for 2ckx, with one mode at distance 6 and the other one at distance 13 of the optimal sequence. This profile is more problematic for sampling methods because of the risk of being trapped far away from the attraction basin of the global optimum. The fitness of the solutions is better correlated to their distance to the global minimum for 2gkt (Figure 3). The linear correlation of 0.2 for 2gkt puts this instance in a class of easy problems according to the scale of Jones [START_REF] Jones | Evolutionary Algorithms, Fitness Landscapes and Search[END_REF]. With a linear correlation of 0.14, 2ckx is a more difficult instance according to this measure. The distance to the optimum decreases with the fitness for both instances, which explains the performance of simulated annealing to find good approximations of the optimal solution on both instances. However, the figure shows a recess at distance 8 of the optimum for 2ckx: this shape suggests that solutions have to break through an energy barrier in order to get closer to the optimal sequence. For 2gkt, the slope is smoother and it seems easier to reach the optimal sequence without having to cross high energy valleys.

Evolvability

The distribution of the number of improving neighbors (i.e., the number of sequences in the neighborhood with a lower energy) is unimodal for both instances (Figure 4). The distribution is narrower for 2ckx, and has a longer tail for 2gkt, which indicates that a solution of the instance 2gkt is easier to improve. Near-optimal solutions are easier to improve for 2gkt. Figure 5 shows the 2D density between fitness values and the number of improving neighbors. The number of improving neighbors decreases with the fitness for both instances. Although there is no linear relation, the correlation is higher for 2gkt: 0.28, whereas it is 0.23 for 2ckx. Figure 5 shows differences in the bottom left quarter: here the number of improving neighbors decreases faster for 2ckx.

Local optima and basins of attraction

The number of local optima is higher for 2ckx (459) than for 2gkt (151). The number of local optima is generally directly linked to the difficulty of the problems. However, in this study, the number of solutions is different for both instances and a direct comparison is difficult. If we only look at the best 497, 282 solutions for 2ckx, a number that fits the sample size of 2gkt, the number of local minima is 220, which is still higher than for 2gkt.

The fitness distribution of local optima has a larger tail for 2gkt (Figure 6), and local optima are more evenly distributed. In contrast, the number of local minima drops down quickly with the quality of the fitness for 2ckx. Figure 7 shows the scatter plot of the basin size normalized by the number of solutions (in log scale) as a function of the fitness value of the corresponding local optima. The logarithm of the basins size is linearly correlated to the fitness of the local optima. The correlation is lost for high energy levels, but this may be an effect of the enumeration threshold: some parts of the attraction basins may not belong to the enumerated ensemble of solutions. The attraction basin sizes tend to be smaller for 2ckx away from the global optimum. This trend is inverted for near-optimal basins. For 2ckx, the cumulative basin size of local optima with high fitness values is larger than those of 2gkt: even taking into account the basin sizes, the near-optimal solutions are more difficult to reach for 2ckx. The length of an adaptive walk (number of steps to reach a local optimum with a steepest-descent algorithm), is positively highly correlated with the logarithm of the corresponding basin size and its fitness value (Figure 8). This property could potentially be used to design a restart strategy. According to the length of steepestdescent, and the fitness value of the local optima found, a restart distance could be estimated in order to maximize the probability to escape from the local optimum, and to find a better local optimum in a next steepest-descent. Moreover, notice that the length of adaptive walks of 2ckx are shorter than the length of 2gkt, and indicates a more difficult multi-modal landscape.

Local optima network

Fig. 9 shows different views of the local optima networks. The graph representing the network of 2gkt is more densely connected than the one of 2ckx. Even though the set of local optima is split in two clusters, the edges show that there is a probability to escape from the basins of attraction of sub-optimal solutions and to reach the cluster leading to the global optimum. For 2ckx, there are several clusters of basins of attraction as well, but no escape edge is visible. Furthermore, the basins of attraction are larger for 2ckx than for 2gkt which confirms that the probability to be trapped into a suboptimal local optimum basin is higher for 2ckx.

The nodes strength of 2gkt is higher than the nodes strength of 2ckx (Figure 10-left): on average 0.142 for 2gkt, and 0.093 for 2ckx. The strength is the sum of out-going weights, and then the opposite of self-loop weights which is the probability to remain in the same basin of attraction. So, the probability to escape from local optima is 1.5 times higher for 2gkt than for 2ckx: it is easier to escape the basin of attraction of 2gkt, in particular for nearoptimal fitness values. This confirms the 2D/3D representation of LON. The local density of the network measured by the weighted clustering coefficients consolidates the picture. 2ckx is locally more clustered (𝑤𝐶𝐶 = 0.103) than 2gkt (𝑤𝐶𝐶 = 0.034), and it should be even more difficult to escape from a cluster of local optima for 2ckx. As a consequence of the local structure of LON, the average path length from one random node to the global minimum node is longer for 2ckx (349.0) than for 2gkt (110.0). The structure of the LON indicates that 2ckx is a more difficult problem than 2gkt: the network of 2ckx is more clustered, with few edges toward better nodes, and a higher probability to stay in the basin of attraction.

Ruggedness

2ckx is slightly more rugged than 2gkt. Fig. 11 (top) shows the number of improving neighbors during random walks computed on the sample with the solution with lowest energy. The autocorrelation length of 2gkt (20) is a little bit longer than 2ckx [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF]. The autocorrelation function of fitness computed from random walks starting from random solution and with an energy bounded to the energy of the initial solution (Figure 12) also confirms this trend. 2gkt seems to be slightly smoother than 2ckx. However, during the random walks, the fitness value suddenly increases or decreases. As a consequence, the sequence of fitness values does not fit an autoregressive model, and another fitness landscape tool should be created to analyze such landscapes. However, this first difference of autocorrelation opens a way to infer the structural difference between the two fitness landscapes from a random solutions, and computationally not expensive sampling. The set 𝑆 ⊂ 𝐼 of a cost function 𝑤 𝑆 is called the scope of the cost function. We denote by 𝐷 𝑆 the Cartesian product of the domains of all variables indexed in 𝑆: 𝐷 𝑆 = 𝑖 ∈𝑆 𝐷 𝑖 . Given a tuple 𝑡 ∈ 𝐷 𝑆 , and 𝑆 ′ ⊂ 𝑆 we denote by 𝑡 [𝑆 ′ ] the projection of 𝑡 on 𝐷 𝑆 ′ . The cost of an assignment 𝑡 of all variables is defined as the sum 𝑤 𝑆 ∈𝑊 𝑤 𝑆 (𝑡 [𝑆]) of all cost functions. If it is strictly less than 𝑘, the assignment is said to be a solution. The weighted constraint satisfaction problem (WCSP) is to identify a solution of guaranteed minimum cost over all 𝑡 ∈ 𝐷 𝑋 . Because of the non negativity of all cost functions, the cost function 𝑤 ∅ ∈ 𝑊 , a constant cost function with no parameters, defines a lower bound on this minimum cost.

CPD modeling

We model CPD with a fixed protein backbone, a set of statistically preferred side chain orientations, or rotamers, and a pairwise decomposable energy function taking the form:

𝐸 (𝑐) = 𝐸 𝑡 + 𝑖 𝐸 (𝑖 𝑟 ) + 𝑖< 𝑗 𝐸 (𝑖 𝑟 , 𝑗 𝑠 )
𝐸 𝑡 is a constant energy terms that captures the internal energy of the protein backbone and the energy of the interactions between the protein backbone and the environment. 𝐸 (𝑖 𝑟 ) captures internal sidechain energies and rotamer-backbone interactions for rotamer 𝑟 at position 𝑖. 𝐸 (𝑖 𝑟 , 𝑗 𝑠 ) captures pairwise interactions between rotamers 𝑟 and 𝑠 at positions 𝑖 and 𝑗 respectively. Because the energy function is pairwise decomposable, all terms can be precomputed and stored in an energy matrix. The objective is to find the combination of rotamers which minimizes the energy of the protein. The side chains differentiate the amino acid types. So, selected rotamers determine the sequence of the protein and solve the CPD problem. With this representation, modeling CPD as a CFN is straightforward and has be done several times in the literature [START_REF] Allouche | Computational protein design as an optimization problem[END_REF][START_REF] Simoncini | Guaranteed Discrete Energy Optimization on Large Protein Design Problems[END_REF][START_REF] Traoré | A New Framework for Computational Protein Design through Cost Function Network Optimization[END_REF]. Each position 𝑖 in the protein sequence defines a variable 𝑥 𝑖 which takes its values in the set of possible rotamers. The precomputed energy terms are converted into costs by shifting them by a constant in order to make them non negative, and by rounding after multiplication by a large constant in order to make them integers. The constant energy term 𝐸 𝑡 is represented as the constant cost function 𝑤 ∅ , rotamer-backbone interactions as unary cost functions and pairwise interactions between rotamers as binary cost functions. The backbone dependent CPD representation and an example of CPD modeling as a CFN are illustrated on Figure 13.

We use ToulBar2 1 to compute the optima and enumerate suboptimal solutions. ToulBar2 explores a search tree with a Branch and Bound algorithm. At the root of the search tree, all variables are unassigned. At each node, the algorithm either assigns a value to a variable (left branch) or removes this value (right branch). Each leaf of the tree is an assignment of all variables. The optimal solution can be obtained by examining the leaves and identifying the assignment that gives the best cost. In order to avoid exploring the whole exponentially sized tree, ToulBar2 maintains an upper bound and a fast incremental lower bound on the cost of solutions to prune branches that cannot lead to better solutions [START_REF] Cooper | Soft arc consistency revisited[END_REF]. In a minimization problem, the cost of the best solution found so far defines the upper bound, and the lower bound provides an underestimation of the best cost below an unexplored branch. If the lower bound reaches the upper bound, the branch is pruned. When every branch has been either explored or pruned, the algorithm returns the optimal solution. In order to perform enumerations, we run our algorithm with a fixed upper bound 𝑘 = 𝑥 * + 𝑟 where 𝑥 * is the optimum and 𝑟 the radius of the enumeration around the optimum. This way, a branch is cut only if it contains solutions whose costs fall outside of the enumeration radius.

1 http://www.inra.fr/mia/T/toulbar2

Models preparation

The structures of two proteins used as targets were downloaded from the Protein Data Bank (PDB ID : 2ckx and 2gkt) and relaxed in all-atom representation with PyRosetta, using the energy function beta_nov16. From the relaxed models, pairwise energy matrices were computed and converted into wcsp format using PyRosetta, the beta_nov16 energy function, the Dunbrack backbone dependent rotamer library [START_REF] Maxim | A smoothed backbonedependent rotamer library for proteins derived from adaptive kernel density estimates and regressions[END_REF] and an in-house Python script. 2gkt is 51 amino acid residues long, 2ckx is 84 amino acid residues long. The global optima were computed from the wcsp files using Toul-Bar2 with the following options: -dee: -O=-3 -B=1 -A -s -cpd -scpbranch. The enumerations were performed by running Toul-Bar2 with the options: -dee: -A -s -hbfs: -cpd -scpbranch -bestconf -a -ub=threshold, where threshold is the cost of the global optimum augmented by the cost of the desired enumeration radius.

For each instance, 1, 000 simulated annealing trajectories were performed using the fixbb program from the Rosetta software with the same energy function and the relaxed protein models as input.

Fitness landscape analysis

In evolutionary computation, fitness landscapes (FL) is one of the fundamental approaches to understand the geometry of the search space from the point of view of local search techniques such as evolutionary algorithm, simulated annealing, etc. On one side, FL depicts the search space with metaphorical pictures such as mountains, peaks, valleys, plateaus, etc. that could help one researcher to design better algorithms; On the other side, it brings a set of numerical metrics that can be used to compare problem difficulty, or be used as features for machine learning techniques that predict algorithm performance. One goal of a fitness landscape analysis is then to contrast optimization difficulty of different problem instances, representations, or neighborhood relations, etc.

Formally, a fitness landscape [START_REF] Stadler | Fitness Landscapes[END_REF] is a triplet (X, N, 𝑓 ) where X is the set of potential solutions of the optimization problem, 𝑓 : X → I R is the fitness function, and N : X → 2 X is the neighborhood relation between solutions: for each solution 𝑥 ∈ X, N (𝑥) is the set of solutions so-called neighbors of 𝑥. Several features of fitness has been proposed and are detailed in the following paragraphs.

The Density Of States (DOS) [START_REF] Rosé | The density of states-a measure of the difficulty of optimisation problems[END_REF] is the distribution of fitness values of solutions from the search space. It has been shown that the speed of decrease of the distribution tail towards good solutions is an indicator of the problem difficulty. The faster the decrease, the higher the difficulty. Indeed, this distribution corresponds to the probability density of a random search to reach a fitness value. In our experiments, in order to compare distributions, the fitness value is shifted in such a way that the global minimum of energy defines an optimal zero fitness.

One of the oldest fitness landscape feature is the Fitness Distance Correlation (FDC) [START_REF] Jones | Evolutionary Algorithms, Fitness Landscapes and Search[END_REF]. The idea of FDC is to measure how the fitness function can guide search towards the global optimum. The FDC is the scatter plot of the fitness vs. the distance of solutions to optimum. The problem becomes easier when the correlation is higher (higher than 0.15 according to the scale proposed by Jones [START_REF] Jones | Evolutionary Algorithms, Fitness Landscapes and Search[END_REF]). The global optimum is required to compute the FDC, which is usually a drawback of the feature. But, for the CPD problem, the global minimum is known. The distance used in this work is the Hamming distance i.e., the number of amino-acid substitutions between two sequences. The distance is clearly related to the local search operator used by existing Simulated Annealing implementations.

The set of FL metrics is related to evolvability [START_REF] Altenberg | The evolution of evolvability in genetic programming[END_REF] which is defined in this context by the ability to evolve towards better solutions. One tool to measure the evolvability is the Fitness Cloud (FC) [START_REF] Verel | Where are bottlenecks in NK fitness landscapes?[END_REF] which is the bi-variate distribution of fitness of solution before and after applying a local search operator 𝑜𝑝: (𝑓 (𝑥), 𝑓 (𝑜𝑝 (𝑥))). Several statistics can be computed from the FC. For example, the number of improving solutions in the neighborhood: 𝑛 + (𝑥) = ♯{𝑥 ′ ∈ N (𝑥) : 𝑓 (𝑥 ′ ) < 𝑓 (𝑥)}. A high number means that the probability to improve the current solution is higher. So, the correlation between the fitness of the solution, and the number of improving neighbors is an indicator of the difficulty of the problem: the probability decreases more slowly for an easier problem.

The main geometries of fitness landscape are the neutral ones dominated by large plateaus of solutions with the same fitness values, and the multi-modal ones with many local optima. The ruggedness of fitness landscapes often impacts the multi-modality of the fitness landscapes. Indeed, the ruggedness is the local continuity of the fitness function thanks to the neighborhood relation. More rugged FL tends to be more multi-modal. The autocorrelation of fitness during a random walks [START_REF] Weinberger | Correlated and uncorrelatated fitness landscapes and how to tell the difference[END_REF] is one of the metric of ruggedness. A random walk is a sequence of solutions (𝑥 1 , 𝑥 2 , 𝑥 3 , . . . , 𝑥 𝑛 ) such that 𝑥 𝑡 +1 is a neighbor of 𝑥 𝑡 . The first solution 𝑥 1 is a random solution from the search space, and the next solution 𝑥 𝑡 +1 of the walk is selected at random in the neighborhood of 𝑥 𝑡 . The autocorrelation function 𝜌 (𝑠) of the fitness is the correlation coefficient of {(𝑥 𝑡 , 𝑥 𝑡 +𝑠 ) : 1 ≤ 𝑡 ≤ 𝑛 -𝑠}. The autocorrelation length ℓ is defined by the smallest integer such that |𝜌 (ℓ)| < 2/ √ 𝑛. It measures the degree of ruggedness of the landscape. A larger autocorrelation length indicates a smoother and easier landscape.

One major feature of multi-modal fitness landscapes is the number of local optima: a high number of local optima usually leads to a difficult problem. However, a number of research works have shown that the link between the number of local optima and difficulty must be contrasted by the size of the basins of attraction [START_REF] Garnier | Efficiency of Local Search with Multiple Local Optima[END_REF]. The basin of attraction of a local optimum 𝑥 ★ is the set of solutions which converge to 𝑥 ★ by steepest-descent search. A steepest-descent is a walk such that 𝑥 𝑡 +1 is the best neighbor of 𝑥 𝑡 according to the fitness function. Intuitively, if the size of the basin of attraction of the global minimum is "small", then the computation effort to reach it is "important". The size of a basin of attraction can be estimated using the average length of steepest-descent towards the local optimum [START_REF] Verel | Pareto local optima of multiobjective NK-landscapes with correlated objectives[END_REF]. It has also been observed on several classes of combinatorial problems (nk-landscapes, Quadratic Assignment Problem, Flow Shop Scheduling Problem) that the fitness of the local optima is correlated with the logarithm of the size of basins of attraction: a deeper local optimum has more chance to have a larger basin of attraction. Notice that the size of global minimum usually remains exponentially small compared to the search size.

Recently, the Local Optima Network (LON) [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF] have been proposed to compress the information of the whole search space into a directed weighted graph where the nodes are the local optima, and the edges are the possible weighted transitions between them. Then, the graph can be analyzed using features from complex networks science. Two definitions of edges have been proposed [START_REF] Ochoa | Local optima networks: A new model of combinatorial fitness landscapes[END_REF]. An edge from the escape-edges type between local optima 𝑖 and 𝑗 is defined if the distance between 𝑖 and a solution of the attraction basin of 𝑗 is below a constant. An edge from the basin-edges type between 𝑖 and 𝑗 exists if there are two solutions 𝑥 𝑖 and 𝑥 𝑗 respectively in the basin of 𝑖 and 𝑗 such 𝑥 𝑗 is a neighbor of 𝑥 𝑖 . The weight is defined by

𝑤 𝑖 𝑗 = 1 ♯𝑏 𝑖 𝑥 ∈𝑏 𝑖 𝑥 ′ ∈𝑏 𝑗 𝑝 (𝑥 → 𝑥 ′ )
where 𝑏 𝑖 and 𝑏 𝑗 are the basins of attraction of 𝑖 and 𝑗, and 𝑝 (𝑥 → 𝑥 ′ ) is the probability to pass from the solution 𝑥 to 𝑥 ′ with the given neighborhood structure. The main metrics of the network which are correlated with the performance of local search algorithm are the average strength, the average path length to the global optimum, and the weighted clustering coefficient [START_REF] Ochoa | Local optima networks: A new model of combinatorial fitness landscapes[END_REF]. The strength of a node is the sum of outgoing weights: 𝑠 𝑖 = 𝑗≠𝑖 𝑤 𝑖 𝑗 . A higher strength means a higher probability to escape from the basin of attraction. The distance 𝑑 𝑖 𝑗 between two local optima 𝑖 and 𝑗 on the network is defined by the inverse of the weight: 𝑑 𝑖 𝑗 = 1 𝑤 𝑖 𝑗 . Intuitively, the weight is the probability to pass from a local optimum to another, so, the distance can be interpreted as the expected number of mutations needed to pass from a local optimum to another. The path length between one local optimum in the network and the global optimum is the length of the shortest path in the graph weighted by the distances 𝑑 𝑖 𝑗 . The local optima network may have more than one connected component. So, to compute the average distance, the harmonic mean is used. The weighted clustering coefficient (wCC) evaluates the local density of the graph. It measures how two neighbors of a node can be also neighbors with each other. A high clustering coefficient shows a high density network which may indicate either a smooth path to global optimum, or, conversely, point out a hard region of the search from which it is difficult to escape. The wCC remains, however, a useful metric to understand the structure of the network. As complementary information, in addition to numerical metrics, it may be useful to represent the local optima using a 3D representation [START_REF] William B Langdon | Visualising the Search Landscape of the Triangle Program[END_REF]. The layout is obtained with Fruchterman-Reingold layout algorithm based on the weighted edges. The width of the edge is proportional to the weight, the size of nodes is proportional to the logarithm of the basin size, and the color and height of nodes is given by the fitness value of the local optima.

CONCLUSION

CPD methods play an important role in protein engineering and have many practical applications in Biotechnology and Synthetic Biology. Metaheuristics are suitable methods for this task, given the size of the search space and the NP-hardness of the problem. However, without any knowledge of the fitness landscape near the global optimum, and without any guarantee on the quality of the solutions found, it is hard to improve the sampling of the search space and to increase the success rate of CPD-based protein engineering. Using some powerful CFN-based exact optimization methods, we were able to exhaustively enumerate all solutions in the vicinity of the global optimum for a CPD problem on which the simulated annealing algorithm implemented in Rosetta, the most famous and commonly used molecular modeling software, behaves poorly. We analyzed the fitness landscape of this CPD problem, and compared it to the fitness landscape of an easier CPD problem on which the simulated annealing algorithm behaves correctly.

Our analysis highlights some crucial differences in the fitness landscapes, explaining the failure of the simulated annealing on the difficult instance. The structure of the landscape, with the presence of several sub-optimal clusters of local optima disconnected from the basin of attraction of the global minimum, is likely the main reason for the bad performance of the simulated annealing. With this type of landscapes, a method able to periodically perturb the solutions in order to help them escape the basin of attraction of a sub-optimal solution would be more suitable. Statistics on random walks also showed some differences between the two instances, in agreement with the fitness landscape analysis, and could be used to assess the difficulty of CPD problems and adapt the sampling strategy before or during sampling.

Figure 1 :

 1 Figure 1: Density of states over fitness.

Figure 2 :

 2 Figure 2: Distribution of distance to the optimum sequence.

Figure 3 :

 3 Figure 3: Correlation between fitness and distance to the optimum sequence (FDC) with color in log scale. 2ckx (left), and 2gkt (right).

Figure 4 :

 4 Figure 4: Distribution of the # of improving neighbors.

Figure 5 :

 5 Figure 5: Correlation between fitness and the number of improving neighbors with color in log scale. 2ckx (left), and 2gkt (right).

Figure 6 :

 6 Figure 6: Fitness distribution of local optima.

Figure 7 :

 7 Figure 7: Scatter plot of the fitness of local optima and their corresponding normalized basin of attraction size. Notice the log scale for the basin size.

Figure 8 :

 8 Figure 8: Scatter plot between the basin size, and fitness of local optima and the average length of steepest descents to reach the local optima. Notice the log scale for the basin size.

3 METHODS 3 . 1

 331 Cost function networks. Cost function networks (CFNs) are deterministic Graphical Models derived from Constraint Satisfaction Problems [20] Definition 3.1. A CFN (𝑋,𝑊 , 𝑘) is defined by: • a set 𝑋 of discrete variables 𝑥 𝑖 ∈ 𝑋 indexed by 𝐼 = {1, . . . , 𝑛}, each variable 𝑥 𝑖 takes its values in a finite domain 𝐷 𝑖 of maximum cardinality 𝑑. • a set of cost functions 𝑤 𝑆 ∈ 𝑊 each involving a subset {𝑥 𝑖 ∈ 𝑋 | 𝑖 ∈ 𝑆 } of all variables and taking non negative values in [0, 𝑘]. • The value 𝑘 is a finite or infinite cost representing an upper bound on costs: a cost of 𝑘 or above is considered as forbidden.

Figure 9 :

 9 Figure 9: Local optima network with basin edges for 2ckx (left side) and 2gkt (right side) in 2D and 3D representations. The size is the log of the basin of attraction size. The color is the fitness: the better, the warmer. The thickness of the edges is proportional to the probability of passing from one basin to the other.

Figure 10 :

 10 Figure 10: Strength and distance to global optimum on the network as a function of fitness.

Figure 11 :

 11 Figure 11: Example of random walks. At each step the number of improving neighbors and autocorrelation function of the random walks are computed.

Figure 12 :

 12 Figure 12: Example of random walks from a random solution and energy bounded to the energy of the initial point. At each step of the walk, the fitness of the solution is computed.

<Figure 13 :

 13 Figure 13: Top: the backbone of the protein is shown in red. The side chains are shown in blue. In the fixed backbone representation, only side chains can move. Statistically preferred side-chain orientations are stored in rotamer libraries. Bottom: toy CFN example representing the interactions between the side chains of two amino acids: a Valine (variable 𝑋 ) and a Leucine (variable 𝑌 ). Values 𝑎 and 𝑏 represent rotamers. Assigning 𝑎 to 𝑋 has a unary cost of 1. The assignments 𝑋 = 𝑎, 𝑌 = 𝑎, 𝑋 = 𝑏, 𝑌 = 𝑎, and 𝑋 = 𝑏, 𝑌 = 𝑏 have binary costs of 1. All other assignments have costs of 0.

ACKNOWLEDGMENTS

This work relied on the HPC resources of the Computing Center of Region Midi-Pyrénées (CALMIP, Toulouse, France). This work has been partially funded by the french Agence Nationale de la Recherche, under Grant ANR-16-C40-0028 and the program EMER-GENCE IDEX Toulouse (project E-CODE).