Keywords: Ad hoc, DTN, vehicular network, information, propagation speed

This research report analyses information propagation speed in bidirectional vehicular delay tolerant networks such as roads or highways. The analysis provided in this report shows that phase transition occurs concerning information propagation speed, with respect to vehicle density in each direction of the highway. Under a certain threshold, information propagates on average at vehicle speed, while above this threshold, information propagates faster. This report gives the exact expression of this threshold and of the average propagation speed in this context, around the threshold. Simulations results, obtained on several platforms, are presented and confirm the analytical results.

Vitesse de Propagation de l'Information dans les Réseaux DTN Véhiculaires Bidirectionnels

Résumé : Ce rapport étudie la vitesse de propagation de l'information dans les réseaux DTN véhiculaires bidirectionnels, tels des automobiles sur une une autoroute. L'analyse présentée dans ce rapport montre qu'une transition de phase se produit concernant la vitesse de propagation de l'information, par rapport à la densité des véhicules dans chaque direction. Sous une certaine borne, l'information se propage à la vitesse des véhicules, alors qu'au dessus de cette borne, l'information se propage plus vite. Ce rapport donne l'expression exacte de cette borne et de la vitesse moyenne de propagation de l'information dans ce contexte, autour de cette borne. Des résultats de simulations, obtenus sur différentes plateformes, sont présentés et confirment ces résultats analytiques.

Mots-clés : Ad hoc, DTN, réseau véhiculaire, information, vitesse de propagation

Introduction

The limits of the performance of multi-hop packet radio networks have been studied for more than a decade, yielding fundamental results such as those of Gupta and Kumar [START_REF] Gupta | The capacity of wireless networks[END_REF] on the capacity of fixed ad hoc networks. These studies assume that either end-to-end paths are available or packets are dropped on the spot. Following seminal works such as [START_REF] Grossglauser | Mobility increases the capacity of ad hoc wireless networks[END_REF] evaluating the potential of mobility to increase capacity, recent research studies focussed on the limits of the performance beyond the end-to-end hypothesis, i.e., when end-to-end paths may not exist and communication routes may only be available through time and mobility. In this context nodes may carry packets for a while until advancing further towards the destination is possible. Such networks are generally referred as Intermittently Connected Networks (ICNs) or Delay Tolerant Networks (DTNs).

Interest in DTN modeling and analysis has risen as novel network protocols and architectures are being elaborated to accommodate various forms of new, intermittently connected networks, which include vehicular ad hoc networks (VANETs), power-saving sensor networks, or even Interplanetary Internet [START_REF] Cerf | Delay-Tolerant Networking Architecture[END_REF].

In this paper, we study the information propagation speed in the typical case of bidirectional vehicular DTNs, such as roads or highways (e.g., about 75% of the total statute miles in the USA [18]). Our objective consists in determining the maximum speed at which a packet (or beacon) of information can propagate in such a bidirectional vehicular network. Our analysis shows that a phase transition occurs concerning information propagation speed, with respect to the vehicle density. We prove that under a certain threshold, information propagates on average at vehicle speed, while above this threshold, information propagates much faster. We provide the exact expressions of the threshold and of the average propagation speed near the threshold.

With applications such as safety, ad hoc vehicular networks are receiving increasing attention (see recent surveys [START_REF] Casteigts | Communication protocols for vehicular ad hoc networks[END_REF][START_REF] Toor | Vehicle ad hoc networks: Applications and related technical issues[END_REF]). Delay tolerant architectures have thus been considered in this context in recent studies, and various analytical models have been proposed. In [START_REF] Wisitpongphan | On the Routing Problem in Disconnected Vehicular Ad-hoc Networks[END_REF], the authors study vehicle traces and conclude that vehicles are very close to being exponentially distributed on highways. The authors of [START_REF] Burgess | MaxProp: Routing for Vehicle-Based Disruption-Tolerant Networks[END_REF] also base themselves on traces gathered in DieselNET (the experimental vehicular network deployed by UMass) to elaborate and evaluate a novel DTN routing algorithm. In [START_REF] Fracchia | Analysis and Design of Warning Delivery Service in Inter-vehicular Networks[END_REF], the authors provide a model for critical message dissemination in vehicular networks and derive results on the average delay in delivery of messages with respect to vehicle density. The authors of [START_REF] Yousefi | Analytical Model for Connectivity in Vehicular Ad Hoc Networks[END_REF] propose an alternative model for vehicular DTNs and derived results on node connectivity, under the hypothesis that vehicles are exponentially distributed. The study is based on queuing theory techniques and characterizes the relationship between node connectivity and several parameters including speed distribution and traffic flow. In [START_REF] Wu | Analytical Models for Information Propagation in Vehicle-to-Vehicle Networks[END_REF], the authors model vehicles on a highway, and study message propagation among vehicles in the same direction, taking into account speed differences between vehicles, while in [START_REF] Nadeem | A Comparative Study of Data Dissemination Models for VANETs[END_REF] authors study message dissemination among vehicles in opposing directions and conclude that using both directions increases dissemination significantly.

Several studies focus on characterizing the packet propagation delay in DTNs: [START_REF] Pellegrini | A Graph-based model for disconnected ad hoc networks[END_REF] which models DTNs as Erdös-Rényi random graphs to derive results concerning packet propagation delay, [START_REF] Zhang | Performance modeling of epidemic routing[END_REF] which uses fluid limit techniques to derive relationships between buffer space, packet duplication and dissemination delay.

Other studies focus on information propagation speed in DTNs. In [START_REF] Kong | On the latency for information dissemination in Mobile Wireless Networks[END_REF] the authors show that when a two-dimensional network is not percolated, the latency scales linearly with the Euclidean distance between the sender and the receiver, while in [START_REF] Jacquet | Information propagation speed in mobile and delay tolerant networks[END_REF], the authors obtained analytical estimates of the constant bounds on the speed at which information can propagate in two-dimensional DTNs. Studies such as [START_REF] Agarwal | Phase Transition Behavior of Message Propagation in Delay Tolerant Vehicular Ad Hoc Networks[END_REF][START_REF] Agarwal | Analytical Model for Message Propagation in Delay Tolerant Vehicular Networks[END_REF][START_REF] Agarwal | Exploiting Mobility to Achieve Fast Upstream Propagation[END_REF] are the closest related work, also focusing on information propagation speed in one-dimensional DTNs. These studies introduce a model based on space discretization to derive upper and lower bounds in the highway model under the assumption that the radio propagation speed is finite. Their bounds, although not converging, clearly indicates the existence of a phase transition phenomenon for the information propagation speed. Comparatively, we introduce a model based on Poisson point process on continuous space, that allows both infinite and finite radio propagation speed, and derive more fine-grained results above and below the threshold (some of the work described in the following was presented in [START_REF] Baccelli | Information Propagation Speed in Bidirectional Vehicular Delay Tolerant Networks[END_REF]). Using our model, we prove and explicitly characterize the phase transition.

In this context, our contributions are as follows: (1) we develop a new vehicule-to-vehicule model for information propagation in bidirectional vehicular DTNs in Section 2; (2) we show the existence of a threshold (with respect to the vehicle density), above which the information speed increases dramatically over the vehicle speed, and below which the information propagation speed is on average equal to the vehicle speed, and (3) we give the exact expression of this threshold, in Section 3; (4) in Section 4, we prove that, under the threshold, even though the average propagation speed equals the vehicle speed, DTN routing using cars moving on both directions provides a gain in the propagation distance, and this gain follows a sub-linear power law with respect to the elapsed time, in the referential of the moving cars; (5) we characterize information propagation speed as increasing quasi-exponentially with the vehicle density when the latter becomes large above the threshold, in Section 5; [START_REF] Casteigts | Communication protocols for vehicular ad hoc networks[END_REF] we cover both infinite radio propagation speed cases, then finite radio propagation speed cases in Section 6, and [START_REF] Cerf | Delay-Tolerant Networking Architecture[END_REF] we validate the provided analysis with simulations in multiple environments (The One and Maple), in Section 7. In the following, we consider a bidirectional vehicular network, such as a road or a highway, where vehicles move in two opposite directions (say east and west, respectively) at speed v, as depicted in Figure 1. Let us consider eastbound vehicle density as Poisson with intensity λ e , while westbound vehicle density is

Model and Results

Westbound Westbound cluster

Eastbound cluster Eastbound

INRIA Figure 2:
Information propagation threshold with respect to (λ e , λ w ). Below the curve, the average information propagation speed is limited to the vehicle speed (i.e., the propagation speed is 0 in the referential of the eastbound cars), while above the curve, information propagates faster on average.

Poisson with intensity λ w . We note that the Poisson distribution is indeed a reasonable approximation of vehicles moving on non-congested highways [START_REF] Wisitpongphan | On the Routing Problem in Disconnected Vehicular Ad-hoc Networks[END_REF]. Furthermore, we consider that the radio propagation speed (including store and forward processing time) is infinite, and that the radio range of each transmission in each direction is of length R.

The main result presented in this paper is that, concerning the information propagation speed in such an environment, a phase transition occurs when λ e R and λ w R coincide on the curve y = xe -x , i.e., λ e Re -λeR = λ w Re -λwR .

(

) 1
Figure 2 shows the corresponding threshold curve for R = 1 (in the following, we will always consider the case R = 1, without loss of generality). We show that below this threshold, the average information propagation speed is limited to the vehicle speed, while above the curve, information propagates faster on average.

We focus on the propagation of information in the eastbound lane. Our aim is the evaluation of the maximum speed at which a packet (or beacon) of information can propagate. An information beacon propagates in the following manner, illustrated in Figure 3: it moves toward the east jumping from car to car until it stops because the next car is beyond radio range. The propagation is instantaneous, since we assume that radio routing speed is infinite. The beacon waits on the last eastbound car until the gap is filled by westbound cars, so that the beacon can move again to the next eastbound car.

We denote T i the duration the beacon waits when blocked for the ith time and D i the distance traveled by the beacon just after. The random variables T i and D i are dependent but, due to the Poisson nature of vehicle traffic, the tuples in the sequence (T i , D i ) are i.i.d.. From now on, we denote (T, D) the independent random variable.

We denote L(t) the distance traveled by the beacon during a time t on the eastbound lane. We consider the distance traveled with respect to the referential of the eastbound cars. We also define the average information propagation speed v p as:

v p = lim t→∞ E(L(t)) t . (2)
By virtue of the renewal processes, we have

v p = E(D) E(T) . (3)
For the remainder of the paper, for x > 0, we denote x * the conjugate of x with respect to the function xe -x : x * is the alternate solution of the equation x * e -x * = xe -x . Notice that x * * = x and 1 * = 1.

We prove the following theorems, concerning the information propagation phase transition threshold (Theorem 1), the distribution of the waiting time spent by information packets in the phase depicted in Figure 3a (Theorem 2), and the total propagation distance achieved (under the phase transition threshold) because of multi-hop bridging using nodes in both traffic directions, as depicted in Figure 3b (Theorem 3).

Theorem 1. For all (λ e , λ w ), the information propagation speed v p with respect to the referential of the eastbound cars is v p < ∞, and,

λ e < λ * w ⇒ v p = 0, (4)
λ e > λ * w ⇒ v p > 0. ( 5
) INRIA Theorem 2. When t → ∞, P (T > t) = A(λ e , λ w )(2vt) -λe λ * w (1 + o(1)) , (6)
for some A(λ e , λ w ), function of (λ e , λ w ).

Notice that Theorem 1 is in fact a corollary of Theorem 2.

Theorem 3. When λ * w > λ e (case v p = 0), when t → ∞, E(L(t)) = B(λ e , λ w )(2vt) λe λ * w + O(t 2 λe λ * w -1 ) . ( 7
)
for some B(λ e , λ w ), function of (λ e , λ w ).

3 Phase Transition: Proof of Theorem 1

Proof Outline

We call cluster a maximal sequence of cars such that two consecutive cars are within radio range. A westbound (respectively, eastbound) cluster is a cluster made exclusively of westbound (respectively, eastbound) cars. A full cluster is made of westbound and eastbound cars.

We define the length of the cluster as the distance between the first and last cars augmented by a radio range. We denote L w a westbound cluster length. We start by proving in Section 3.2 that the Laplace transform of L w : f w (θ) = E(e -θLw ) equals:

f w (θ) = (λ w + θ)e -λw-θ θ + λ w e -λw-θ , (8)
and the exponential tail of the distribution of L w is given by

P (L w > x) = Θ(e -λ * w x ). (9)
To evaluate how information will propagate according to Figure 3, we compute the distribution of the gap length G e between the cluster of eastbound cars on which the beacon is blocked and the next cluster of eastbound cars. We show in Section 3.4 that P (G e > x) = O(e -λex ). Now, let T(x) be the time needed to meet a westbound cluster long enough to fill a gap of length x (i.e., a westbound cluster of length larger than x). We show in Section 3.3 that:

E(T(x)) = Θ( 1 vP (L w > x) ) = Θ(e λ * w x ) . (10)
Therefore, the average time T to get a bridge over all possible gaps is

E(T) = ∞ 1 E(T(x))e -xλe dx = 1 2v ∞ 1 Θ(exp((λ * w -λ e )x))dx . ( 11
)
As a result, the threshold with respect to (λ w , λ e ) where E(T) diverges is clearly when we have:

λ * w = λ e , (12)
or, in other words, since λ * w e -λ * w = λ w e -λw , when we have:

λ w e -λw = λ e e -λe . (13)

Cluster Length Distribution

Lemma 1. The Laplace transform of the westbound cluster length f w (θ) = E(e -θLw ) satisfies:

f w (θ) = (λ w + θ)e -λw-θ θ + λ w e -λw-θ . ( 14
)
Proof. The length of the cluster is counted from the first car. The random variable L w satisfies:

• L w = 1, with probability e -λw , when the first car has no car behind within the radio range;

• L w = g w + L w , with probability 1 -e -λw , where g w is the distance to the next car and g w < 1.

Translating this in terms of Laplace transforms yields:

f w (θ) = e -λw-θ + f w (θ) 1 0 λ w e -(λw+θ)x dx (15)
= e -λw-θ + f w (θ)

λ w λ w + θ (1 -e -λw-θ ) . (16)
In passing, we get E(L w ) = -f w (0) = e λw -1 λw . Lemma 2. We have the asymptotic formula:

P (L w > x) = (λ w -λ * w )e λ * w -λw (1 -λ * w )λ * w e -λ * w x (1 + o(1)) (17)
Proof. The asymptotics on P (L w > x) are given by inverse Laplace transform:

P (L w > x) = - 1 2iπ -ε+i∞ -ε-i∞ f w (θ) θ e θx dθ = - 1 2iπ -ε+i∞ -ε-i∞ (λ w + θ)e -λw (θe θ + λ w e -λw )θ e θx dθ ,
for some ε > 0 small enough. For (θ) < 0 the denominator (θe θ + λ w e -λw ) has two simple roots at θ = -λ * w and θ = -λ w and is absolutely integrable elsewhere. The root -λ w does not lead to a singularity since it is canceled by the numerator λ w + θ. The residues theorem neutralizes the pole at -λ w , therefore for some ε > 0:

P (L w > x) = (λ w -λ * w )e λ * w -λw (1 -λ * w )λ * w e -λ * w x + O(e -(λ * w +ε)x ) . INRIA

Road Length to Bridge a Gap

Now, let us assume that we want to fill a gap of length x. We want to know the average length of westbound road until the first cluster that has a length greater than x -1. Figure 4 depicts a gap of length x, and the length of westbound road until a cluster is encountered which can bridge the gap. Let f w (θ, x) be the Laplace transform of the cluster length, under the condition that it is smaller than x: f w (θ, x) = E(1 (Lw<x) e -θLw ).

(starting from arbitrary cluster)

v v Lw1 < x -1 Lw2 < x -1 Unbridged gap length x Lw3 > x -1 R = 1
Road length to bridge gap Bw(x) Lemma 3. The Laplace transform of the road length B w (x) to bridge a gap of length x, starting from the beginning of an arbitrary cluster, is:

(a) R = 1 v v Lw3 > x -1 (b)
β w (θ, x) = E(e -θBw(x) ) = P (L w > x -1) 1 -λw λw+θ f w (θ, x -1) , (18)
and,

E(B w (x)) = 1 + O(e -εx ) e λw λ w (1 -λ * w )λ * w (λ w -λ * w )e λ * w -λw e (x-1)λ * w . ( 19
)
Proof. Before a cluster of length greater than x appears there is a succession of clusters L 1 , L 2 , . . . each of length smaller than x (see Figure 4b). Several cases are possible:

• the first cluster to come is a cluster of length greater than x, with probability P (L w > x -1); the road length to this cluster is equal to 0.

• the first cluster greater than x is the second cluster; in this case, the road length Laplace transform is equal to f w (θ, x) λw λw+θ , i.e., the Laplace transform of a cluster multiplied the Laplace transform of the exponentially distributed inter-cluster distance, namely λw λw+θ (as Laplace transform multiplication represents the addition of the corresponding independent random variables).

• or, in general, the first cluster greater than x is the kth cluster; in this case, the road length Laplace transform is equal to

f w (θ, x) λw λw+θ k
Therefore, the Laplace transform of the road length to the cluster of length greater than x is equal to the sum of the Laplace transforms of the previous cases, i.e., P

(L w > x -1) ∞ k=0 f w (θ, x) λw λw+θ k .
Thus, the average is (combining with Lemma 2):

E(B w (x)) = - ∂ ∂θ β w (0, x) = - ∂ ∂θ f w (0, x -1) - 1 λ w f w (0, x -1)
× 1 P (L w > x -1) = e λw λ w + O(e -(x-1)λ * w ) 1 P (L w > x -1) = 1 + O(e -εx ) e λw (1 -λ * w )λ * w λ w (λ w -λ * w )e λ *
w -λw e (x-1)λ * w

Gap Distribution

Let us call G e an eastbound gap which is not bridged (see Figure 5). As illustrated in Figure 6, G e can be decomposed into a westbound cluster length L * w without eastbound cars, plus a random exponentially distributed distance I e to the next eastbound car.

Unbridged eastbound gap Ge

Eastbound Westbound Bridged eastbound gap Ge which is defined for all (θ) > -λ e , and

E(G e ) = - f w (λ e ) f w (λ e ) + 1 λ e . ( 21
)
Proof. Let p w (x) be the probability density of a westbound cluster length L w .

The probability that a westbound cluster has no eastbound cars is:

∞ 0 p w (x)e -λex dx = f w (λ e ). (22)
The Laplace transform of the westbound cluster length without eastbound cars, defined for all (θ) > λ * w + λ e , equals:

E(e -θL * w ) = ∞ 0 p w (x)e -λex e -θx dx P r(no eastbound) = f w (θ + λ e ) f w (λ e ) . ( 23
)
and the Laplace transform of G e = L * w + I e follows, as well as its average estimate.

Lemma 5. The probability density p e (x) of G e is:

p e (x) = λ e f w (λ e ) e -λex (1 + O(e -εx )) . (24)
Proof. The proof comes from a straightforward singularity analysis on the inverse Laplace transform.

Distribution of Waiting Time T

Lemma 6. We have 2vT = L * w + I w + B w -1 (where I w is a random exponentially distributed distance to the next westbound car), and, therefore,

2vE(T) = E(L * w ) -1 + 1 λ w + ∞ 1 E(B w (x))p e (x)dx . ( 25
)
Proof. The total distance to bridge a gap, as depicted in Figure 7, equals the distance to the beginning of the first westbound cluster (L * w + I w ) plus the road length to bridge a gap starting from an arbitrary cluster (B w ) minus 1, since communication can start at exactly one radio range. Since the distance is covered by cars moving in opposite directions, we have 2vT = L * w + I w + B w -1. We complete the proof by taking the expectations, and averaging on all possible gap lengths x.

Corollary 1. The quantity E(T) converges when λ e > λ * w and diverges when λ e < λ * w . Proof. The proof comes from the leading terms of E(B w (x)) and p e (x).

Distance D Traveled after Waiting Time T

We denote C e the distance traveled beyond the first gap. As depicted in Figure 8 Proof. The random variable C e is smaller in probability than a full cluster. Lemma 8. The average value of C e satisfies:

E(C e ) = 1 λ e 1 -f w (λ e ) f w (λ e ) + f w (λ e ) f w (λ e ) . ( 26
)
Proof. From [START_REF] Zhang | Performance modeling of epidemic routing[END_REF], the probability that an eastbound car is not bridged to the next eastbound car equals f w (λ e ). The unconditional gap length is 1 λe . We define Ḡe as an eastbound gap, under the condition that the gap is bridged (see Figure 5). Therefore, the average gap length Ḡe satisfies the following:

f w (λ e )E(G e ) + (1 -f w (λ e ))E( Ḡe ) = 1 λ e , ( 27
) INRIA which gives E( Ḡe ) = 1 λe + f w (λe)
1-fw(λe) . The distance C e traveled in bridging (beyond the first gap, and extended to the next cluster, which is eventually bridged) satisfies:

E(C e ) = (1 -f w (λ e )) E( Ḡe ) + E(C e ) (28) = 1 λ e 1 -f w (λ e ) f w (λ e ) + f w (λ e ) f w (λ e ) . ( 29
)
Corollary 2. The total distance D e traveled including the first gap satisfies

E(D e ) = E(G e ) + E(C e ) =
1 λefw(λe) , which remains finite for all vehicle densities.

Since E(D e ) is finite (Corollary 2) and E(T) converges when λ e > λ * w , and diverges when λ e < λ * w (Corollary 1), we obtain the proof of Theorem 1.

4 Power Laws: Proof of Theorems 2 and 3

Waiting Time Distribution

In this section, we are interested in finding an evaluation of the waiting time distribution P (T > y), when y → ∞, in case λ e < λ * w , i.e., when the information propagation speed is 0 on average. Proof. See appendix.

Lemma 10. When t → ∞:

P (T > t) = A(λ e , λ w ) (t2v) -λe λ * w (1 + o(1)) . (30)
Proof. We have the relation

T 2v = G ew + B w , with G ew = L * w + I w -1.
We know that G ew , in analogy with G e , has an exponential tail, i.e., E(e -θGew ) < ∞ for all θ > -λ w . In other words P (G ew > y) = O(exp(-θy)). The other particularity is that G ew and B w are dependent. First we have the inequality for all y P (G ew + B w > y) ≥ P (B w > y), therefore, we have P (T > t) ≥ A(λ e , λ w ) (t2v)

-λe λ * w (1 + o(1)
). Second, we have the other inequality for all (y, z):

P (G ew + B w > y) ≤ P (G ew > z) + P (B w > y -z) .

Thus, by selecting

z = O(log t) such that P (G ew > z) = o(t -λe λ * w ), we get P (T > t) ≤ A(λ e , λ w ) (t -O(log t)2v) -λe λ * w + o(t -λe λ * w )
Therefore, we conclude the proof of Theorem 2.

RR n°7266

Traveled Distance Distribution

Now, we focus on the renewal process made of the various waiting time intervals T, experienced by the information beacon. Considering the sequence of waiting phases, T 1 , T 2 , . . . , T n , . . .: the beacon moves at time T 1 , then at time T 1 +T 2 , T 1 + T 2 + T 3 , etc. We denote n(t) the number of phases achieved before time t:

i=n(t) i=1 T i ≤ t < i=n(t)+1 i=1 T i .
Since the T i are independent, this is a renewal process, and we have the identity:

P (n(t) ≤ n) = P (T 1 + • • • + T n ≥ t) .
We can get a precise estimate of the average number of renewals E(n(t)) during time t. Lemma 11. There exists b > 0 such that, when t → ∞, the following estimate is valid:

E(n(t)) = sin 2 (π λe λ * w ) bπ 2 Γ( λ e λ * w )t λe λ * w + O(t 2 λe λ * w -1 ) . (31)
Proof. See appendix.

In parallel to the sequence of waiting intervals {T i } i≥1 , we have the sequence {D i } i≥1 the distances traveled by the beacon after every waiting interval T i . Now we denote L(t) = i=n(t) i=1 D i which is the total distance traveled by the beacon until time t.

Lemma 12. E(L(t)) = E(n(t))E(D) .

Proof. We have the identity:

E(L(t)) = n>0 E(1 n(t)≥n D n ), ( 32
)
where 1 n(t)≥n is the indicative function of the event n(t) ≥ n. Since, from the definition of n(t), we have the equivalence:

n(t) ≥ n ⇔ T 1 +• • •+T n-1 < t, then 1 n≤n(t)
and D n are independent random variables and, therefore,

E(1 n(t)≥n D n ) = P (n(t) ≥ n)E(D).

Quantity E(D) has a closed expression (E(D) =

1 λefw(λe) , from Corollary 2. Substituting and using Lemma 11, the power law for E(L(t)) in Theorem 3 is shown.

Asymptotic Estimates

Near the Threshold

v p ∼ 2v (λ w -λ * w )λ w λ 2 e (1 -λ * w )λ * w (λ e -λ * w )e λ * w +λe-2λw . ( 33
)
Proof. Using Lemma 6, we have:

2vE(T) = - f w (λ e ) f w (λ e ) + ∞ 1 E(B w (x))p e (x)dx ∼ ∞ 1 e λw λ w (1 -λ * w )λ * w λ w -λ * w e λw-λ * w e (x-1)λ * w × λ e f w (λ e )
e -λex dx .

Integrating, we obtain v p = E(D) E(T) .

Large Densities

Corollary 4. When the vehicle densities become large, i.e., λ e , λ w → ∞:

v p ∼ 2v e λe+λw 1 + λw λe + λe λw . (34)
Proof. According to Lemma 4, we have:

E(L * w ) = 1 + λ w λ e (λ w + λ e ) , ( 35
)
and the expected gap length tends to 1. The average road length to bridge such a gap tends to 1 λw . From Lemma 6:

2vE(T) ∼ E(L * w ) -1 + 1 λ w . (36)
From corollary 2, the average distance traveled in bridging is:

E(D) = 1 λ e f w (λ e ) ∼ e λe+λw λ e + λ w , (37)
and we obtain v p = E(D) E(T) . Note that the information propagation speed grows quasi-exponentially with respect to the total vehicle density.

Finite radio propagation speed

If the radio propagation speed (including store and forward timings) is finite and constant, equal to v r , then the average information propagation speed becomes:

s = E(T w )v + E(D e ) E(T w ) + 1 vr E(D e ) . (38)
Figure 9: Threshold zones for (λ e , λ w ) for various γ = v vr , γ = 0 red, γ = 0.1 blue, γ = 0.5 green, γ = 0.8 yellow.

But the speed v r impacts the random variables G e and T . The main impact is that, to fill an eastbound gap of length x, one needs a westbound cluster of length at least x(1+γ) 1-γ with γ = v vr , otherwise the message will be off the gap when arriving at the end of the cluster. Thus, E(T w (x)) = O(e (1+γ)θw 1-γ

x ). The gap length is also modified, since we must consider 1-γ 1+γ L * w :

E(e -θGe ) = f w ( 1-γ 1+γ (θ + λ e )) f w ( 1-γ 1+γ λ e ) λ e λ e + θ . (39)
But, this does not change the exponential term e -λex in the asymptotic expression of p e (x). Therefore, the threshold condition becomes:

λ w e -λw = 1 -γ 1 + γ λ e e -1-γ 1+γ λe . (40)
This corresponds to a dilatation by a factor 1+γ 1-γ of the horizontal axis in the diagram of Figure 2. Notice how the diagram then loses its symmetry with respect to λ e versus λ w , which can be observed in Figure 9.

Notice also that, when v r decreases to v, the threshold limit tends to infinity.

Similarly, in correspondance to Theorem 3, we have E(L(t)) = Ω(t

1-γ 1+γ λe λ *
w ), due to the dilatation by the factor 1+γ 1-γ .

INRIA

In this section, we present simulation results obtained with Maple on one hand, and the Opportunistic Network Environment (ONE [START_REF] Keranen | The ONE Simulator for DTN Protocol Evaluation[END_REF]) simulator on the other hand.

Maple Simulations

We first compare the theoretical analysis with measurements performed using Maple. In this case, the simulations follow precisely the bidirectional highway model described in Section 2: we generate Poisson traffic eastbound and westbound traffic on two opposite lanes moving at constant speed, which is set to v = 1m/s. The radio propagation range is R = 1m, and radio transmissions are instantaneous; the length of the highway is sufficiently large to provide a large number of bridging operations (of order at least 10 3 ) for all considered traffic densities. We measure the information propagation speed which is achieved using optimal DTN routing, by selecting a source and destination pairs at large distances, taking the ratio of the propagation distance over the corresponding delay, and averaging over multiple iterations of randomly generated traffic. We vary the total traffic density, and we plot the resulting information propagation speed. Figures 10 and11 show the evolution of the information propagation speed near the threshold versus the total vehicle density, when λ e = λ w , in linear and semilogarithmic plots, respectively. We can observe the threshold at λ e +λ w = 2 in Figure 10, which confirms the analysis presented previously in Section 3, and corresponds to λ e = λ w = 1 in Figure 2. In semilogarithmic scale (Figure 11), we observe that the simulation measurements quickly approach a straight line, and therefore are close to the theoretically predicted exponential growth above the phase transition threshold, in Section 5. In Figure 12, we present a 3dimensional plot of the eastbound information propagation speed v p by varying the vehicle densities in both eastbound (λ e ) and westbound (λ w ) traffic.

Finally, we perform detailed measurements of the waiting time T that each packet of information spends in the buffer of an eastbound car until it encounters a westbound cluster which allows it to propagate faster to the next eastbound car. For the measurements, we set the vehicle densities to λ e = λ w = 0.9; thus the conjugate λ * w = 1.107 . . . . In Figure 13, (i.e., below the phase transition threshold), we plot the distribution of the waiting time T, and we compare it to the predicted power law in Theorem 2:

P (T > t) = t -λe λ * w .

ONE Simulations

In this section, we depart from the exact Poisson model simulations in Maple, and we present simulation results obtained with the Opportunistic Network Environment (ONE [START_REF] Keranen | The ONE Simulator for DTN Protocol Evaluation[END_REF]) simulator. Vehicles are distributed uniformly on the length of both lanes of a road, and move at a constant unit speed. The total number of vehicles varies from 1000 to 5000. Similarly to the previous section, we measure the fastest possible information propagation speed which is achieved using epidemic broadcast, assuming that radio transmissions are instantaneous and that there are no buffering or congestion delays, with a radio range R = 10m. Again, we vary the vehicle densities λ e and λ w , which are given in vehicles per radio range, and we perform several simulation iterations of randomly generated traffic. As shown in Figure 14, similarly to the Maple simulation results and to the analysis, we observe the threshold phenomenon at λ e = λ w = 1: the information propagation speed remains almost constant below the threshold but increases dramatically beyond it. We also observe an exponential growth above the threshold in Figure 15.

We remark that measurements below the phase transition threshold yield an average information propagation speed which is slightly larger than the vehicle speed. This is due to the finite duration of the simulations and in the computations of the expectations. This phenomenon can also be explained from the theoretical analysis in Section 4: even below the threshold, DTN routing using cars moving on both directions still provides a gain in the propagation distance, which follows a sub-linear power law with respect to the elapsed time (in the referential of the moving cars). Figure 16 confirms the predicted power law t λe/λ * w for E(L(t)), the average distance traveled by information with respect to time (in Theorem 3), which is shown in Figure 16 for λ e = 1 and λ w = 0.9, therefore, a growth of order t 0.903 . vehicle speed, and below which information propagation speed is on average equal to vehicle speed (in Theorem 1), and we computed the exact expression of this threshold. We exactly characterized the information speed near the threshold (Corollary 3), and we showed that, above the threshold, the information propagation speed increases quasi-exponentially with vehicle density (Corollary 4). We also analyzed in detail the way information propagates under the threshold, and we showed that DTN routing using bidirectional traffic provides a gain in the propagation distance, which follows a sub-linear power law with respect to the elapsed time (Theorems 2 and 3). Combining all these different situations, we obtain a complete image of the way information propagates in vehicular networks on roads and highways, which is useful in determining the performance limits and designing appropriate routing protocols for VANETs. All our theoretical results were validated with simulations in several environments (The One and Maple).

Our analysis can be extended to investigate other models of vehicle traffic and radio propagation. In future works, we intend to provide a detailed expression of this threshold in specific VANET models (e.g., intersections). Finally, an interesting direction for further research consists in collecting large traces of real traffic on roads and highways, and evaluating the information propagation properties in this context. Proof. We want to find an evaluation of P (T > y) when y → ∞. For this we will evaluate for x given P (B w (x) > y). We know that

INRIA

E(e θBw(x+1) ) = β w (θ, x + 1) = λ w P (L w > x) θ -λ w (f w (θ, x + 1) -1) Since P (L w > x) = α λ * w e -λ * w x (1 + o(e -εx )), with α = (λw-λ * w )e λ * w -λw 1-λ * w , we have f w (θ, x) = f w (θ) - αe -λ * w x θ + λ * w (1 + o(e -εx )) (41)
From here we drop the o(e -εx ) for simplicity, since it will just bring an exponentially small factor. Therefore we have which is in power law as claimed.

Proof of Lemma 11

We first prove the property for an hypothetic renewal process based on the B i s on the real line. Let b(x) be this renewal process at time t: when θ → 0. This is also true for complex θ.

We can now prove Lemma 11 for the original renewal process.

Proof. From the previous lemma, we have that: N B (θ) = sin(πa) bπ θ -a-1 +O(θ -1 + θ -2a ), when θ → 0.

The inversion of the Laplace transform yields

E(b(y)) = 1 2iπ c+i∞ c-i∞ N B (θ)e θy dθ .
for any c > 0. we bend the integration path so that it resemble the path described in figure 17 with c < 0. Since far from its singularities on the line -θ(x) for x > 0 (which corresponds to the negative real axis), the functions β(θ, x) is uniformly in The cranted part of the integral, using the Flajolet Odlyzko theorem on continuous functions [START_REF] Flajolet | Singularity analysis of generating functions[END_REF], gives a contribution which is sin 2 (πa) bπ 2 Γ(-a)t a + O(1 + t 2a-1 ).

Figure 1 :

 1 Figure 1: Model of a bidirectional vehicular network on a highway.

Figure 3 :

 3 Figure 3: Eastbound information propagation: the beacon waits on the last eastbound car (a), until the gap is bridged by westbound cars so that the beacon can move again (b).

Figure 4 :

 4 Figure 4: Illustration of the road length B w (x) until a gap x is bridged: (a) smaller clusters cannot bridge the gap, (b) until a westbound cluster of length at least x -1 is encountered.

Figure 5 :Lemma 4 .

 54 Figure 5: Illustration of a bridged gap Ḡe , and an unbridged gap G e .

Figure 6 :

 6 Figure 6: Unbridged gap G e model; L * w corresponds to a westbound cluster length without eastbound cars.

Figure 7 :

 7 Figure 7: Waiting time T: the total distance to bridge a gap is L * w +I w +B w -1.

Figure 8 :Lemma 7 .

 87 Figure 8: Total distance D traveled when a bridge is created D = G e + C e .

Lemma 9 .

 9 When y tends to infinity, P (B w > y) = A(λ e , λ w )y -λe λ * w (1 + o(1)), with A(λ e , λ w ) = λ e e -λe λ * w f w (λ e ) .) is the Euler "Gamma" function, β = λw(λw-λ * w )e λ * w -2λw 1-λ * w .

First 3 .

 3 , we investigate the case where (λ e , λ w ) is close to the threshold boundary.INRIACorollary When λ e → (λ * w ) + :

Figure 10 :

 10 Figure 10: Maple simulations. Information propagation speed v p for λ e = λ w , versus the total vehicle density λ e + λ w .

Figure 11 :

 11 Figure 11: Maple simulations. Information propagation speed v p for λ e = λ w , versus the total vehicle density λ e + λ w , in semi-log scale, compared to the theoretically predicted asymptotic exponential growth.

Figure 12 :

 12 Figure 12: Maple simulations. Eastbound information propagation speed v p for different values of the vehicule densities, λ e and λ w .

Figure 13 :

 13 Figure 13: Maple simulations. Cumulative probability distribution P (T > t) of the waiting time T, compared to the power law t -λe λ * w , for λ e = λ w = 0.9.

Figure 14 :

 14 Figure 14: ONE simulations. Information propagation speed for λ e = λ w , with respect to λ e + λ w .

Figure 15 :

 15 Figure 15: ONE simulations. Information propagation speed for λ e = λ w , with respect to λ e + λ w , in semi-log scale, compared to the theoretically predicted asymptotic exponential growth.

Figure 16 : 8 . 1

 1681 Figure16: ONE simulations. Average distance traveled by information with respect to elapsed time E(L(t)), divided by the predicted power law t λe/λ * w for λ e = 1 and λ w = 0.9, compared to the constant value 1.80 (dash).

θ

  λ w (f w (θ) -1αe -λ * w (x + 1) = y) = 1 2iπβ w (θ, x + 1)e yθ dθ .Let θ(x) be the root of θλ w (f w (θ) -1αe -λ * analysis gives θ(x) = -βe -λ * w x + O(e -2λ * w x ), with β =P (B w (x + 1) = y) = λ w αe -λ * w x λ * w × e θ(x)y 1 -λ w ∂ ∂θ f w (θ(x))αe -λ * w x (θ(x)+λ * w ) 2+O(e (θ(x)-ε)y )Omitting the O( ) terms we getP (B w (x + 1) = y) = λ w αe -λ * = βe -θx exp(-βe -λ * w x y) , Or P (B w (x + 1) > y) = exp(-βe -λ * w x y) .(43) Therefore stating P (B w > y) = p e (x)P (B w (x) > y)dx we get, omitting O( ) terms P (B w > y) = λ e f w (λ e ) e -(x+1)λe exp(-βe -λ * w x y)dx . (44) RR n°7266 with the change of variable u = e -λ * w x we get P (B w > y) = λ e e -λe λ * w f w (λ e )

ELemma 13 .

 13 (b(y)) = P (B 1 < y) + P (B 1 + B 2 < y)+ • • • + P (B 1 + • • • + B n < y) + • • • Let us define β = 1 -E(e -θBw ). We also define N B (θ) = ∞ 0 E(b(y))e -yθ the Laplace transform of E(b(y)). We have:N B (θ) = 1 -B(θ) θB(θ) . When θ → 0, then B(θ) = b πθ a sin(πa) for some b and a = λe λ * w θ θ + y g(θ, θ -1 (y))p e (θ -1 (y))(θ -1 (y)) dy .We use together the following estimates θ(x) = -βe -λ * w x + O(e -2λ * w x ), and p e (x) = λe fw(λe) e -λex (1 + O(e -εx )) to statep e (θ -1 (y))(θ -1 (y)) = λ e f w (λ e ) (yβ ) λe λ * w -1 (1 + O(y ε )) , INRIA for some ε > 0. Then, we use the fact that ∞ Therefore, B(θ) = b πθ a sin(πa) + O(θ) for some b and a = λe λ * w

1 θFigure 17 :

 117 Figure 17: The Flajolet-Odlyzko cranted integral path.

RR n°7266

Concluding RemarksIn this paper, we provided a detailed analysis for information propagation in bidirectional vehicular DTNs. We proved the existence of a threshold, concerning vehicle density, above which information speed increases dramatically over INRIA

Therefore, we obtain that, for y → ∞: E(b(y)) = sin 2 (πa) bπ 2 Γ(-a)y a + O(1 + y 2a-1 ).

To terminate the argument we will prove that E(n(t)) = E(b(2vt)) + o(1). We have 2vT = B w +H w , where H w = L w +I w -1, and H w ≥ 0 has distribution with an exponentially decreasing tail. However, we stress the fact that B w and H w are not independent. We denote H i the successive values of the H w . We have therefore 2vT i = B i + H i .

Since T i ≥ 2vB i , we already have E(n(t)) ≤ E(b(2vt)). We also have for any t ∈ [0, t] and for any integer i:

Using Chernoff bounds, we have for all θ > 0 such that E(e θH ) exists, ∀i ≤ k and ∀x > 0:

On the other side, we also have ∀y > 0:

We have

We use the fact that P (B > y) ≤ Ay -a for a = λe λ * w and for some A > 0 to state that

and, finally,

Collecting all results yields:

If we take t = t a with 1 > a > a, k = 2vt γ+E(H) = Ω(t a , then we have

which is exponentially small when t → ∞. Therefore,


Maria Suarez

Pablo Tortosa

Maria M Garcia-Mira

David Rodríguez-Larrea

Raquel Godoy-Ruiz

Beatriz Molero

Alfonso Jaramillo
email: alfonso.jaramillo@polytechnique.fr

Jose M Sanchez-Ruiz

David Rodríguez- Larrea

Beatriz Ibarra-Molero

Using multi-objective computational design to extend protein promiscuity

Keywords: protein design, protein promiscuity, computational chemistry, multipurpose catalysts, Pareto Set

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T

INTRODUCTION

Enzymes are natural protein catalysts that can enhance reaction rates up to 23 orders of magnitude [START_REF] Kraut | Challenges in enzyme mechanism and energetics[END_REF] with an impressive affinity and/or specificity. Their applications are enormous, and enzyme design holds the promise to provide important impact in society areas like medicine (i.e. treatment of neurodegenerative diseases [START_REF] Lecanu | Cutting-edge patents in Alzheimer's disease drug discovery: anticipation of potential future treatments[END_REF], design of therapeutically antibodies to bind tumor-associated antigenic determinants while maintaining a small immunogenicity, development of peptide-based vaccines), biotechnology (biosensors [START_REF] Marvin | The rational design of allosteric interactions in a monomeric protein and its applications to the construction of biosensors[END_REF][START_REF] Looger | Computational design of receptor and sensor proteins with novel functions[END_REF], biocatalysts with activity in non-natural environments)

and bioremediation (design of enzymes that would reduce waste by-products and toxicity [START_REF] Anga | Recent advances in the bioremediation of persistent organic pollutants via biomolecular engineering[END_REF]).

In the last few years, studies on enzyme evolvability have given convincing evidence on the mechanism of evolution of protein functions from cross-reactive proteins [START_REF] Jensen | Enzyme recruitment in evolution of new function[END_REF][START_REF] O'brien | Catalytic Promiscuity And The Evolution Of New Enzymatic Activities Chemistry and Biology[END_REF][START_REF] Aharoni | The 'evolvability' of promiscuous protein functions[END_REF][START_REF] Copley | Enzymes with extra talents: moonlighting functions and catalytic promiscuity[END_REF][START_REF] Khersonsky | Enzyme promiscuity: evolutionary and mechanistic aspects[END_REF],

starting from the observation that most proteins possess, besides their native function, additional promiscuous functions with specificities in the range of k cat /K m ~10 -2 -10 6 [START_REF] Khersonsky | Enzyme promiscuity: evolutionary and mechanistic aspects[END_REF]. In the proposed mechanism, a weak promiscuous function arises due to neutral evolution, protein robustness (the ability of proteins to tolerate mutations without compromising fitness) and plasticity (the ability of gaining new functions by a reduced number of mutations). Under the right selection pressure, natural selection can improve the new function once it has arisen until, at some point, the protein may become specialized for the new function.

In the classic view, selection constraints on the native function were believed to be determinant in the process (if a given function appeared by natural selection, there must be a large penalty associated to the loss of that function) reducing the allowed evolution scenario to gene duplication at an early stage of the specialization. However, modern studies have shown that in most cases the coupling of the nascent and the original

A C C E P T E D M A N U S C R I P T

functions is smaller than expected (see [START_REF] Khersonsky | Enzyme promiscuity: evolutionary and mechanistic aspects[END_REF] and references therein). This allows for several generations of "generalist" proteins [START_REF] Rothman | How Does an Enzyme Evolved in vitro Compare to Naturally Occurring Homologs Possessing the Targeted Function?. Tyrosine Aminotransferase from Aspartate Aminotransferase[END_REF] that are able to perform both functions, and suggests that gene duplication acts after the new function has appeared [START_REF] Aharoni | The 'evolvability' of promiscuous protein functions[END_REF] and not before.

Proteins with several activities (multipurpose enzymes) may have a wide range of biotechnological applications related (but not limited) to industrial organic synthesis and metabolic engineering [START_REF] Jürgens | Directed evolution of a (βα)8-barrel enzyme to catalyze related reactions in two different metabolic pathways[END_REF][START_REF] Kazlauskas | Enhancing catalytic promiscuity for biocatalysis[END_REF][START_REF] Canada | Directed evolution of toluene ortho-monooxygenase for enhanced 1-naphthol synthesis and chlorinated ethene degradation[END_REF][START_REF] Zheng | Mechanistic studies on the alkyltransferase activity of serotonin N-acetyl transferase[END_REF][START_REF] Spreitzer | RUBISCO: structure, regulatory interactions, and possibilities for a better enzyme[END_REF][START_REF] Bornscheuer | Catalytic Promiscuity in Biocatalysis: Using Old Enzymes to Form New Bonds and Follow New Pathways[END_REF]. However, natural protein promiscuity (as described in the two preceding paragraphs) may perhaps be of limited use in the development of these multipurpose catalysts. First of all, natural promiscuous activities are often related to the evolved ones, sharing the same active sites and even the basic chemical mechanisms and, in general, bearing a significant resemblance to the original function [START_REF] Khersonsky | Enzyme promiscuity: evolutionary and mechanistic aspects[END_REF]. Secondly, development of the promiscuous activity upon suitable mutations usually brings about a decrease in the evolved activity (see Table 2 in [START_REF] Khersonsky | Enzyme promiscuity: evolutionary and mechanistic aspects[END_REF]). We explore in this work the use of computational design to overcome these limitations

We have selected thioredoxin from E. coli (PDB code 2trx, 1.5 Å resolution [START_REF] Katti | Crystal structure of thioredoxin from Escherichia coli at 1.68 Å resolution[END_REF]) as scaffold for our studies. It is a small (108 residues) general disulphide oxidoreductase found in all the kingdoms of living organisms; it is a common model for protein design studies because of its high stability and good expression properties [START_REF] Ladbury | Stability of oxidized Escherichia coli thioredoxin and its dependence on protonation of the aspartic acid residue in the 26 position[END_REF]. We aim at introducing an esterase activity [the nucleophilic hydrolysis of the p-nitrophenyl acetate (PNPA) into p-nitrophenol (PNP) and acetate] in E. coli thioredoxin, following an approach similar to that used in ref [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF]. Unlike these previous studies [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF], however, we intend to preserve the natural thioredoxin activity (see below for details).

Natural promiscuous activities appear to be shaped by residues at the "wall and perimeter" of the native active site [START_REF] Khersonsky | Enzyme promiscuity: evolutionary and mechanistic aspects[END_REF]. These residues show high-plasticity, likely because they do not belong to the protein scaffold or the native catalytic machinery [START_REF] Khersonsky | Enzyme promiscuity: evolutionary and mechanistic aspects[END_REF] A C C E P T E D M A N U S C R I P T and provide a suitable target for the introduction of new, non-natural promiscuous activities [START_REF] Yoshikumi | Designed divergent evolution of enzyme function[END_REF][START_REF] Bone | Structural plasticity broadens the specificity of an engineered protease[END_REF]. However, designing a new active-site implies introduces "unsatisfied" destabilizing interactions (which will hopefully be satisfied upon ligand and/or transition-state binding). Designing a new active-site in close proximity to the native one poses the additional problem that the introduced destabilizing interactions may disrupt the native active site and affect the original activity. Our computational design approach, therefore, is based on a multi-objective optimization. Both a measure of protein stability and de-novo catalytic activity are simultaneously optimized by using two competing score functions for folding free energy and binding free energy of the protein-ligand (i.e. transition-state-model) complex , obtaining the Pareto Set [START_REF] Papalambros | Principles of Optimal Design[END_REF] of optimal stability/promiscuous-function solutions. The goal of this procedure is the development of the new activity with the lowest possible stability cost, which has two advantages. First of all, the need to introduce additional stabilizing mutations (to compensate the destabilizing effect of the new active-site mutations) is minimized. This is an important point, since the original (native) active-site would be one obvious target for stabilizing mutations (active-site residues are optimized for function, not for stability) and, in this case, we aim at preserving the original activity. Indeed, in a previous work [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF], the catalytic D26 residue was mutated to isoleucine, a change which enhances stability but which will impair the natural thioredoxin oxido-reductase activity. Secondly and most important, since our designed active-site is spatially close to the original one, the low-stability-cost strategy guarantees the minimal perturbation required for maintaining the original activity .

Our computational design approach is based upon the DESIGNER software [START_REF] Jaramillo | Folding free energy function selects native-like protein sequences in the core but not on the surface[END_REF][START_REF] Wernisch | Automatic Protein Design with All Atom Force-fields by Exact and Heuristic Optimization[END_REF], which optimizes protein sequence for a given target structure. This procedure uses atomic models and rotamer libraries to represent side-chain conformations. The

A C C E P T E D M A N U S C R I P T

free energies of the different models are calculated with the CHARMM force-field [START_REF] Mackerell | All-atom empirical potential for molecular modeling and dynamics studies of proteins[END_REF] and a free-energy solvation term proportional to the surface area [START_REF] Ooi | Accessible Surface Areas as a Measure of the Thermodynamic Parameters of Hydration of Peptides[END_REF]. Calculations on a reference state (taken to represent the unfolded state) allow a quantity akin to the unfolding free energy to be computed for each model.

The DESIGNER program was originally developed to address the inverse folding problem. However, as we show in this work, it can also be used to design new activesites. This requires that the design is approached as an optimization in sequence-space of a protein structure which includes a model of the transition state of the chemical reaction (the tetrahedral PNPA intermediate in the case of interest here) with reference to the structure without the transition-state model bound, thus yielding a quantity akin to the free energy barrier of the reaction. Furthermore, since DESIGNER also leads to a quantity akin to the unfolding free-energy, the multi-objective stability-function optimization is indeed feasible.

MATERIALS AND METHODS

Molecular modeling of the new active site

General approach. We aim at introducing a promiscuous esterase activity in E.

coli thioredoxin: the nucleophilic hydrolysis of the p-nitrophenyl acetate (PNPA) into pnitrophenol (PNP) and acetate. The general approach we use is similar to that described in ref. [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF]; i.e., we choose a histidine residue as a nuclophile for the reaction and we model the tetrahedral transition state for the reaction as PNPA-histidine structure constructed as a generalized rotamer of the histidine. An initial exploration of the available positions to place the anchor nucleophile, and their different energies lead us to fix our search to position 39. This position has additional advantages: good accessibility properties and proximity to the native thioredoxin active site [START_REF] Vohnik | Conformation, stability, and active-site cysteine titrations of Escherichia Coli D26A thioredoxin probed by Raman spectroscopy[END_REF], a high A C C E P T E D M A N U S C R I P T plasticity region. Also, position 39 is close to wild-type charged residues (36K, 57K) which could in principle be electrostatically complementary to p-nitrophenyl acetate.

However, this does not represent a bias in our computational analysis, as all mutations are allowed at both positions and lysine is not favored in the rotamer construction process. All subsequent computations described in this work have been carried out on the A39H background.

The new binding site is in principle defined by the positions neighboring to 39, excluding those with wild-type proline or glycine : 1, 4, 26, 28, 30, 36, 38, 42, 43, 46, 55, 57, 59 and 99, which are allowed to mutate to any amino acid except proline, cysteine and methionine (see figure 1A). All other positions were forced to conserve their wild-type amino acid but were free to change rotamer to avoid border effects. No previous stabilizing mutation was introduced. This gave a total of 15 designed positions and 93 fixed amino acid positions, yielding a sequence space of 10 18 .

Side-chain rotamer library enhancement for tetrahedral transition-state modeling. We do not introduce in the force field biasing extra terms in order to obtain interactions known to favor the reaction. On the other hand, we have enhanced the rotamer library to improve the modeling of protein-ligand (transition state model) interactions. We have used the backbone-dependent rotamer library by Dunbrack and Karplus [START_REF] Dunbrack | Backbone-dependent rotamer library for proteins: application to side-chain prediction[END_REF] including all rotamers with probability P>90%. The rotamer of the combined histidine (Nδ as nucleophile) and tetrahedral intermediate was modeled with the same initial dihedral angles as in [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF]: χ 1 and χ 2 from the histidine in the library, χ 3 (±30,±90,±150), χ 4 (±-60,±180), and χ 5 (±90). To provide a better docking of the ligand to thioredoxin, we have added a refined version of the standard library following the protocol described in [START_REF] Tortosa | Active sites by computational protein design[END_REF]. This procedure refines the library specifically in the ligand-protein interface by searching only for conformations that provide a better

A C C E P T E D M A N U S C R I P T

interaction and improve ligand recognition. For each pair of rotamers involved in binding we first minimize the structure (100 iterations of steepest descent) of the interacting rotamers. A low-dielectric environment (ε=1) in CHARMM was found to significatively speed-up calculations by biasing the search towards H-bond formation.

During minimization, backbone charges are set to zero to restrict the optimization to rotamer-rotamer interactions. In a second step, we minimize the resulting structure (100 steps of ABNR, ε=8) to obtain the refined conformation of both residues, which are stored as new rotamers for further use. To simplify the size of the calculations, we required a maximum interacting energy of -5.0 kcal/mol and performed a clustering in dihedral angle with a granularity factor of 10 degrees. This resulted in a total number of 2886 rotamers and a size of the combinatorial space of 10 95 possible structures.

Computational analyses. We performed the initial exploration of the available positions to place the anchor nucleophile which included the enhancement of the library described above and resulted in a CPU time consumption of 8000h in processors Power4 and Power4+ (1.3GHz and 1.7GHz) using the facilities provided by IDRIS supercomputer center. Once we fixed the anchor nucleophile position (H39) and performed the clustering in dihedral angle of the resulting library, we proceeded to recalculate the interaction energies of the remaining 2886 rotamers, (2000 CPU hours provided by the BSC supercomputer center: Power PC 970MP processors at 2.3 GHz).

The optimization procedure took 1000 CPU hours at BSC.

Site-directed mutagenesis, protein expression, purification and sample preparation

Oligonucleotides used for mutagenesis were purchased from Invitrogen Corp.

Mutations were introduced in the sequence of wild-type thioredoxin from E. coli, The protein was purified by size exclusion chromatography and anion exchange chromatography as previously described [START_REF] Perez-Jimenez | The efficiency of different salts to screen charge interactions in proteins: a Hofmeister effect?[END_REF]. The purity of the protein was checked by SDS-PAGE electrophoresis. Samples were prepared by exhaustive dialysis against 10 mM sodium phosphate, pH 6.95 buffer. Protein concentration was determined by UV absorbance (ε 280 = 8310 M -1 cm -1 ).

previously

Activity measurements

Two kinds of activity measurements were carried out in the present study. The nucleophilic hydrolysis of p-nitrophenyl acetate (PNPA) into p-nitrophenol (PNP) was performed as previously described in [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF] with minor modifications. All experiments were carried out at 22 ºC. Stock solutions of 75 mM PNPA were freshly prepared in acetonitrile by weight. The reaction was initiated by the addition of PNPA to a solution of thioredoxin previously tempered. PNP formation was followed by absorbance at 400.5 nm (ε 400.5 = 12400 M -1 cm -1 ) . In all cases, kinetics were measured against PNPA solutions in buffer, used as a blank to take into account the spontaneous hydrolysis of the reagent. Initial rates were calculated from the slope of the time-dependence of the absorbance and data were analyzed by non-linear least-squares fitting of the Michaelis-Menten equation (see Figure 3) using SigmaPlot (SPSS Inc.). Additionally, the natural oxidoreductase activity was determined using a turbidimetric assay consisting of the thioredoxin-catalyzed rate of reduction of insulin disulfides by dithiothreitol [START_REF] Holmgren | Thioredoxin catalyzes the reduction of insulin disulfides by dithiothreitol and dihydrolipoamide[END_REF].

Reaction mixtures contained thioredoxin (up to 6-8 µM) in potasium phosphate 0.1 M,

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT 10 EDTA 2 mM pH 6.5 and insulin 0.5 mg/mL. Reaction was initiated by addition of DTT, 1 mM final concentration, previously dissolved in water. The increase in the A 650 due to insulin aggregation upon thioredoxin mediated reduction was monitored for 1 h at 25 ºC (see Figure 4). Activities were calculated from the slopes of the maximum dA 650 /dt versus thioredoxin concentration. In all cases, absorbance measurements were corrected for the value given by a sample lacking thioredoxin and DTT.

Stability measurements

The thermal unfolding transitions were followed by differential scanning calorimetry of V3 against dialysis buffer at a scan rate of 90 K•h -1 . A VP-DSC microcalorimeter (Microcal Inc.) was used. Rescans of V3 were done to check the reversibility of the transition. Experimental data were analyzed on the basis of the twostate model as previously described [START_REF] Georgescu | Heat capacity analysis of oxidized Escherichia coli thioredoxin fragments (1--73, 74--108) and their noncovalent complex. Evidence for the burial of apolar surface in protein unfolded states[END_REF].

RESULTS AND DISCUSSION

Multi-objective optimization and construction of the Pareto Set

We approach enzyme design as a two-objective problem. We consider as one objective the effect on stability of the designed region, which is coupled to the primary function of thioredoxin due to spatial proximity. It is estimated from an approximation to the folding free energy as defined previously [START_REF] Jaramillo | Folding free energy function selects native-like protein sequences in the core but not on the surface[END_REF] and does not involve the presence are not specifically addressed by this procedure, since they are given by complex expressions including the equilibrium and rate constants for the several steps of the catalytic mechanism: see figure 3 in ref [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF].

To rank the sequences by simultaneously using two competing scores (akin to the folding free energy and the binding free energy of the transition-state complex) we use a method previously presented [START_REF] Suarez | Pareto optimization in computational protein design with multiples objectives[END_REF], based on standard techniques in computer science for multi-objective optimization. In a single objective protein design problem, the goal is to find the solution associated with the lowest cost function (the global minimum).

However, in a multi-objective optimization problem there is no single optimal solution, but a set of optimal solutions, each of them representing a different trade-off between the objectives. To compare two different solutions, a relation of dominance is defined as follows: A dominates B if, for all objectives, A has a better performance than B; as a result A has to be kept while B has to be discarded (even though B may lead to an active foldable protein). The set of all non-dominated solutions is called Pareto globally optimal set (or, simply, the Pareto Set) and finding this set is the goal of the multiobjective computational optimization problem [START_REF] Papalambros | Principles of Optimal Design[END_REF]. Note that, in the Pareto Set, only those solutions showing an optimal trade-off between the different objectives are kept.

As a result, the complexity of the problem is dramatically reduced without loss of relevant information (see illustration in Figure 2A).

In order to construct the Pareto Set, we have followed the procedure we have recently described in detail [START_REF] Suarez | Pareto optimization in computational protein design with multiples objectives[END_REF]. Briefly, we have used a protocol of Monte Carlo-Simulated Annealing (MCSA) cycles [START_REF] Gardiner | Digital computer studies of cell multiplication by Monte Carlo methods[END_REF][START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], combined with an exponential cooling schedule (initial and final temperature T 0 =1 R kcal/mol and T f = 0.01 R kcal/mol) to minimize a weighted sum of the scores corresponding to the two objectives:

∆G scalar (λ) = λ ⋅ ∆G stability + (1-λ) ⋅ ∆G newfunction (1)

A C C E P T E D M A N U S C R I P T

where ∆G stability and ∆G newfunction are the computational measures akin to the freeenergies of folding and binding of the transition-state model, respectively. The weight, λ, was allowed to take many different values in the range [0.05, 0.95] to ensure the exploration of different design regions. The rational behind this procedure is that each of the different λ values leads to a given point belonging to the Pareto Set of solutions (although neighboring λ values may lead to the same point). However, the use of λ values different from 0.5 produces solutions in a sub-optimal conformation for that sequence, so an additional minimization by an unbiased MCSA (λ =0.5) was performed for each of the obtained sequences to ensure they reach their minimum energy conformation. The Pareto Set of non-dominated solutions thus obtained is the result of the two-objective optimization problem and is shown in Figure 2B. Note that the tradeoff between the two competing objectives is clearly apparent in the Pareto Set, the most stable sequences showing poor binding and the best binding sequences being highly unstable.

Design of a thioredoxin variant with a promiscuous activity

The Pareto Set provides an ordered and non-redundant list of solutions, each of them with a different trade-off between the chosen objectives. We have examined this list in the intermediate stability/binding range (see Figure 2B) to identify relevant mutations for a three-mutation candidate for experimental validation. A consensus analysis immediately suggests mutations S1E, W28K and K57E (see alignment shown as Inset in Figure 2B). However, the Pareto Set shown in Figure 2 has been constructed on the basis of the average of the four possible conformers associated to the A39H mutation and the subsequent side-chain modeling (i.e., two different chiralities of the acetate and two possible ways of histidine attack). A more detailed examination reveals that the performance of the mutation S1E depends strongly of the conformer associated

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

13 to the A39H mutation, while mutations W28K and K57E are more robust in this regard.

These two mutations, together with the nucleophile insertion (A39H) were modeled with DESIGNER and seen to be reasonably close to the Pareto Set (see Figure 2B). We therefore chose this triple mutant variant (subsequently referred to as V3) for experimental validation. The structure predicted by DESIGNER for our mutant V3 is shown in Figure 1C, with our three designed mutations (W28K, A39H and K57E)

together with surrounding wild-type residues. It illustrates the balance between binding interactions with the ligand (36K, 28K) and stabilizing interactions with neighboring residues (1S-43D, 28K-57E, 30E-36K, 39H-57E).

Native and promiscuous activities of the designed thioredoxin variant

We have determined the promiscuous esterase activity by following the hydrolysis of PNPA into PNP catalyzed by our designed thioredoxin variant, V3, at 22 ºC in sodium phosphate buffer 10 mM, pH 6.9. We carried out the experiments using several concentrations of V3 as well as PNPA. The formation rate of PNP was found to be linear over a wide time interval. The reaction initial rates measured at different enzyme concentration values have been represented in Figure 3A. No deviations from the expected linear dependences with the enzyme concentration were found. The dependence of the reaction rate on PNPA concentration is consistent with a simple Michaelis-Menten mechanism (see Figure 3B) with K m =1.198•10 -3 M and k cat =3.69•10 -4 s -1 (for the uncatalyzed reaction k uncat =2.5•10 -6 s -1 , ref. [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF]). It must be noted that the values of K m and k cat are given by complex expressions including the equilibrium and rate constants for the several steps of the catalytic mechanism: see figure 3 in ref [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF].

As expected, wild type thioredoxin does not show any significant esterase activity.

In addition, the native oxidoreductase activity of thioredoxin and V3 was experimentally determined (Figure 4) by following the enzymatic-catalyzed reduction

A C C E P T E D M A N U S C R I P T

of insulin by DTT as described by Holmgren [START_REF] Holmgren | Thioredoxin catalyzes the reduction of insulin disulfides by dithiothreitol and dihydrolipoamide[END_REF]. Reduced insulin aggregates with a concomitant increase in solution turbidity that is reflected in the increase in absorbance at 650 nm observed after a certain lag-time (Figure 4A). Activities were calculated as the maximum of the dA 650 /dt slopes. These values change linearly with protein concentration (in the 0.2-3 µM concentration range). V3 is not only found to have native reductase activity but, remarkably, it shows a 1.5-fold enhancement in native activity with respect to wild type thioredoxin (Figure 4B). This enhancement is also evident in the lag-time values, which are significantly reduced in the case of V3 (Figure 4C). Recent single-molecule force-clamp spectroscopy studies [START_REF] Wiita | Probing the chemistry of thioredoxin catalysis with force[END_REF] have revealed two alternative forms of the catalytic mechanism of E. coli thioredoxin catalysis, the first requiring a re-orientation of the substrate disulfide bond and the second elongating the substrate disulfide bond. Given this complexity in the catalytic process, it does not seem possible at this stage to propose a reliable molecular interpretation of mutation effects on bulk-phase activity (such as the 1.5-fold enhancement observed with V3). The important point in the context of this work is, however, that the introduction of the promiscuous esterase activity in thioredoxin has not compromised the native oxidoreductase activity. This result does support that our multi-objective stability/binding design is a minimal-perturbation approach, less likely to lead to mutations affecting the original function, even when the high-plasticity positions close to the original active-site are used to create the second active-site associated to the new function (as we have done in this work).

A C C E P T E D M A N U S C R I P T

The balance between stability and promiscuous function in the designed thioredoxin variant

The thermal stability of V3 was determined by differential scanning calorimetry (DSC) at pH 7.0 in 10 mM phosphate buffer (see Figure 3C). The calorimetric transition was found to be reversible and well-described by the two-state reversible model.

Denaturation temperature (71.8 ºC) is significantly lower than that for wild type thioredoxin under the same solvent conditions: 86.6 ºC [START_REF] Georgescu | Heat capacity analysis of oxidized Escherichia coli thioredoxin fragments (1--73, 74--108) and their noncovalent complex. Evidence for the burial of apolar surface in protein unfolded states[END_REF]. Unfolding enthalpy is also found to be lower: 71 kcal/mol in V3 versus 114 kcal/mol in wild type thioredoxin.

However, these unfolding enthalpy values belong to the corresponding denaturation temperatures and, therefore, a significant part of the difference is to be attributed to the temperature dependence of the unfolding enthalpy.

Creating a new active-site in a protein involves introducing potential interactions that will only be satisfied upon substrate (or transition-state) binding. These unsatisfied interactions will have a destabilizing effect in the non-ligated protein. However, for a desing based on the Pareto Set of optimal stability/binding solutions, the stability decrease observed should correspond to the achieved binding affinity (reflecting an optimal trade-off between the competing objectives). We show below that this correspondence does qualitatively hold for the designed esterase activity in the V3 variant of thioredoxin.

Consider the binding of a ligand to a protein. The binding free energy is related to dissociation constant (K d ) through:

      - = ∆ d BINDING K RT G 1 •ln 0 (2) A C C E P T E D M A N U S C R I P T
where the superscript º indicates a standard state free-energy change. Standard state is defined by the pressure and concentration units associated to K d [START_REF] Murphy | Entropy in biological binding processes: estimation of translational entropy loss[END_REF]. As is custormary, we use K d in molar units for free-energy calculations and, therefore, our standard concentration state is 1 molar.

The binding free-energy can be split into two contributions, corresponding to the newly formed interactions in the protein-ligand complex (∆G INT ) and to the loss of translational entropy associated to reduction in the number of particles in solution upon binding (∆Gº TR ):

0 0 TR INT BINDING G G G ∆ + ∆ = ∆ (3)
where only the translational contribution depends on the standard-state concentration unit. Freire, Amzel and coworkers [START_REF] Murphy | Entropy in biological binding processes: estimation of translational entropy loss[END_REF][START_REF] Siebert | Loss of translational entropy in molecular associations[END_REF] have provided evidence that translational free-energies can be calculated to an acceptable approximation on the basis of Kauzmann's cratic entropies [START_REF] Kauzmann | Some factors in the interpretation of protein denaturation[END_REF]:

CRAT TR S T G ∆ - = ∆ • 0 (4)
where, for a P+L→PL process and 1 molar standard state, ∆S CRAT is -

8 cal•K -1 •mol -1 .
Equations 2-4 allow us to estimate the free-energy associated to protein-ligand interactions (∆G INT ) from the experimental value of the dissociation constant (K d ).

Thus, using as K d the Michaelis constant value determined from the fitting shown in figure 3B (K m =1.2•10 -3 M) we obtain -6.4 kcal/mol as an estimate of the free energy associated to the PNPA-V3 interactions. We assume now that those interactions are not satisfied in the non-ligated V3 variant and, therefore, that they contribute to its destabilization with respect to wild type thioredoxin. Schellman and coworkers [START_REF] Becktel | Protein stability curves[END_REF] showed many years ago that, to a first approximation, the effect of a free energy perturbation (g) on equilibrium denaturation temperature can be obtained from:

A C C E P T E D M A N U S C R I P T m m m H T g T ∆ = ∆ • (5)
where T m is the "non-perturbed" denaturation temperature value and ∆H m is the nonperturbed denaturation enthalpy at the temperature T m . For thioredoxin under the solvent conditions of the stability determinations reported here, T m =86.6 ºC=359.8 K and ∆H m =114 kcal/mol. Then, using -6.4 kcal/mol as the value of the free-energy perturbation, we obtain from equation 5 that ∆T m =-20 degrees. This value is in excellent qualitative agreement with the fact that the denaturation temperature for V3 is 15 degrees lower than that for wild type thioredoxin (see Figure 3B).

The above calculation is meant mainly for illustration. Thus for instance, it neglects unfolded-state effects and it has been carried using the PNPA binding free energy derived from the Michaelis constant value, which may in fact be given by a complex expression including equilibrium and rate constants for the serveral steps of the catalytic mechanism (see figure 3 in ref. [START_REF] Bolon | Enzyme-like proteins by computational design[END_REF]). Furthermore, it would have been desirable to perform the calculation on the basis of the activation free-energy value (i.e.

the free energy value for transition-state "binding"). This, however, would have required the assumption of a value for the front-factor in transtion-state-theory rate equation, which would have added considerable uncertainty to the calculation. Caveats notwithstanding, the calculation suggests in a clear manner an important general point, which is ellaborated below:

We may expect a clear difference between re-designing a native active site for a new function [START_REF] Allert | Computational design of receptors for an organophosphate surrogate of the nerve agent soman[END_REF] and introducing a second active-site for the new function. In the first case, destabilizing interactions consistent with the new function are substituted for those corresponding to the original function. In the second case, the destabilizing interactions associated to the new-function are introduced without eliminating those involved in the

A C C E P T E D M A N U S C R I P T
original "native" active site (which is to remain functional) and, therefore, creating a second active site will very likely cause a significant decrease in protein stability which may compromise folding. Note that, according to Schellman equation (eq. 5) a freeenergy perturbation of the order of several kcal/mol (required for order-of-magnitude enhancement in binding affinity or catalytic rate) translates into a decrease in denaturizing temperature of several tens of degrees. Certainly, design of protein promiscuity should benefit from previous stability enhancement but, as illustrated by this work, it also requires the minimization of the associated stability penalty by using an efficient computational multi-objective methodology that selects the optimal tradeoff solutions.

Concluding remarks

During evolution, many different protein features may simultaneously be under selective pressure; these include properties related to stability, catalysis, interactions, solubility, prevention of aggregation, etcetera. Another example of the interplay between multiple features in evolution is the development of a new function from a preexisting weak promiscuous function; in this case, to avoid gene duplication at early stages of the development of the new activity, the protein has to maintain both activities during a certain evolutionary stage.

In contrast, computational protein design usually employs a single-objective scoring function to rank the different sequences in the optimization procedure. However, the need to satisfy multiples objectives in protein engineering has always been realized, but rarely exploited, perhaps because of the lack of general methodological guidelines to the multi-feature optimization problem. This work makes proposals regarding such guidelines and provides the relevant experimental validation. In particular, this work

A C C E P T E D M A N U S C R I P T

shows with a specific experimental scenario how the computational analysis of several competing objectives can be posed to lead to a set of optimal solutions (not a single global minimum), with each possible criterion to choose a specific solution representing a necessary trade-off between the competing objectives.

From a different viewpoint, multi-objective analyses, such as that described here, can offer valuable insight on the diversity of the solution space, the choices that natural evolution confronts when selecting enzymes under different pressures and the role of enzyme promiscuity in shaping evolution.

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T

  cloned into pET-30a, using the QuikChange (TM) Site-Directed Mutagenesis method from Stratagene ©. The V3 triple mutant (W28K/A39H/K57E) A C C E P T E D M A N U S C R I P T was over expressed in BL21 (DE3) pLys during 4 hours. Cells were harvested and resuspended in 30mM Tris, 1mM EDTA, pH 8.3 (TE buffer), lysed by sonication and centrifuged at 10,000×g for 10 min. Nucleic acids in the soluble fraction were precipitated overnight by addition of a streptomycin sulphate 10% solution (Sigma-Aldrich) prepared in TE buffer and removed by centrifugation at 100,000×g for 30 min.

  of the PNPA tetrahedral intermediate. The other objective is a quantity akin to the activation free energy of the reaction which is estimated from calculations on a protein structure which includes a model of the transition state of the chemical reaction (the tetrahedral PNPA intermediate in the case of interest here) with reference to the structure without the transition-state model bound. Note that the values of K m and k cat A C C E P T E D M A N U S C R I P T

FIGURE 1 .

 1 FIGURE 1.-Structural analysis of the design of a promiscuous activity in thioredoxin.

FIGURE 2 .

 2 FIGURE 2.-Multi-objective computational optimization and construction of the Pareto

FIGURE 3 .

 3 FIGURE 3.-Promiscuous activity and stability of the designed V3 variant of

FIGURE 4 .

 4 FIGURE 4.-Native oxidoreductase activity of the wild-type form and the designed V3

  FIGURE 2A

Acknowledgments P.T. acknowledges the financial support from EMBO long-term fellowship program and FP6-NEST Pathfinder Synthetic Biology project num. 043340. A.J. acknowledges the financial support from HPC-Europa, from the "Programa de Incentivo a la Investigación (Universidad Politécnica de Valencia)" and the use of the BSC and IDRIS supercomputer facilities to perform the calculations reported here. This research was supported by Feder Funds and Grant BIO2006-07332 (Spanish Ministry of Education and Science) to J.M.S.-R. and Grant CVI-1668 (Junta de Andalucía) to B.I.-M.

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T