Xiao Hu

Dr Paul Thevenon

Dr Christophe Macabiau

Attitude Determination and RTK Performances Amelioration Using Multiple Low-Cost Receivers with Known Geometry

Since 2017, he is a Ph.D. student within the SIGnal processing and NAVigation (SIGNAV) research group of the TELECOM laboratory inside ENAC. His current research interests focus on Precise GNSS Positioning and hybridization of GNSS with other sensors.

INTRODUCTION

The concept of autonomous driving has become more and more a central topic for the automobile industry, where a precise position and attitude information is essential for this kind of application. Over the past decade, the universal GNSS has been dramatically utilized in various domains, such as aviation, marine, precise agriculture, geodesy and surveying, automotive, etc. However, the accuracy or integrity that a low-cost GNSS receiver can provide in a constrained urban or indoor environment is far from satisfactory for applications where decimeter or centimeter accuracy and error bounds are mostly envisioned. To reach this level of accuracy, techniques using raw carrier phase measurements have been developed. Carrier phase measurements are more precise than code measurements by a factor of a hundred [START_REF] Misra | Global positioning system: signals, measurements, and performance[END_REF]. Nevertheless, they are also less robust and include a so-called integer ambiguity that requires to implement a specific integer ambiguity resolution (IAR) process in order to use them for positioning. In some harsh environments, severe multipath and losses of lock of the receiver tracking loops create carrier cycle slips, which result in sudden changes of these ambiguities. If not detected, the cycle slips will create a bias on the carrier pseudorange measurement resulting in a reduction of position accuracy. Even if detected, the IAR process has to be, at least partially, re-initialized, leading also to a loss of positioning accuracy. To increase confidence and accelerate the IAR process by limiting the search space, restrictions can be established from the use of an array of two or more receivers with prior known and fixed geometry which includes the length of the baseline vectors between the antennas of the receiver array and the orientation of the vectors.

In recent years, several studies have focused on the use of an array of receivers for attitude determination [START_REF] Iafrancesco | GPS/INS Tightly coupled position and attitude determination with low-cost sensors Master Thesis[END_REF], Gabrielle et al. [START_REF] Gnss | GNSS-Based Attitude Determination[END_REF] studies the viability of attitude determination by using a new ambiguity-attitude estimator which improves the probability of successful integer ambiguity resolution. Nonetheless, the improvement of RTK positioning performance was not studied. Daniel et al. [START_REF] Medina | On the Kalman filtering formulation for RTK joint positioning and attitude quaternion determination[END_REF] developed a method for the recursive estimation of the positioning and attitude problems using GNSS carrier phase observations from an array of receivers, but they calculated the position of each receiver separately thus they did not take advantage of the known geometry. Farhad et al. [START_REF] Aghili | Attitude determination and localization of mobile robots using two RTK GPSs and IMU[END_REF] used an adaptive KF for 3-dimensional attitude determination and position estimation of a mobile robot by fusing the information from a system of two RTK GPSs and an IMU, however, they also did not consider the known geometry of the receiver as a constraint to help improve the positioning performance. Zheng et al. [START_REF] Zheng | Enhancing Ambiguity Resolution Performance Using Attitude Determination Constraints[END_REF] presented a methodology for integrating carrier phase attitude determination and positioning systems by considering one of the receiver pairs in the attitude determination system also used as the rover for the relative positioning system. Nevertheless, their attitude determination and positioning systems remained independent which did not much ameliorate the success rate of IAR for the RTK. Khodabandeh et al. [START_REF] Khodabandeh | Single-Epoch GNSS Array Integrity: An Analytical Study[END_REF] introduced a concept of array-based between-satellite single difference satellite phase biases determination to accelerate the single-receiver IAR, but they did not take into consideration the attitude information of the vehicle which cannot analyze the influence with attitude consideration. Peirong Fan et al. [START_REF] Fan | Precise and Robust RTK-GNSS Positioning in Urban Environments with Dual-Antenna Configuration[END_REF] proposed a dual-antenna constraint RTK algorithm to improve the system AR success rate, which combines GNSS measurements of both antennas by making use of the baseline vector constraint between them. However, the attitude information of the vehicle is still not taken into account.

To the author's best knowledge, very few publications can be found that address the use of an array of receivers to improve the accuracy of the array position or for some internal steps of precise position computation for RTK processing with vehicle attitude determination, such as cycle slip detection or integer ambiguity resolution. The objective of this contribution is then to explore the possibility of achieving precise positioning with a low-cost architecture: using multiple low-cost receivers with known geometry to enable the vehicle attitude determination and RTK performances amelioration. The concept of RTK positioning using an array of GNSS receivers has already been introduced in our previous work [START_REF] Hu | Improvement of RTK performances using an array of receivers with known geometry[END_REF]. In [START_REF] Hu | Cycle-slip Detection and Repair Using an Array of Receivers with Known Geometry for RTK Positioning[END_REF], we discussed the improvement of cycle-slip detection and repair for RTK processing by using an array of two receivers with known geometry.

In this paper, while we refer to our earlier work, the focus is different. In both our previous contributions, the proposed precise position and attitude determination algorithm was verified with only simulated measurements to enable the conduction of several specific comparison scenarios between our multi-receiver system and the traditional single receiver system. The present paper is the extension and improvement of our previous work, which aims at validating its application to actual situations by using the real measurement collected from different environments with different length of array antenna baselines, in terms of the positioning accuracy, the fixing rate, and the attitude determination accuracy.

SYSTEM GEOMETRY AND CONFIGURATION

With the intention of performing precise attitude estimation, a dual antenna set-up has been considered, where two GNSS antennas with a known baseline length are mounted on the vehicle's rooftop to get an attitude estimation. Furthermore, the absolute position accuracy is augmented using the real-time kinematics (RTK) approach, in which the vehicle is positioned relative to a third receiver as the virtual reference station (VRS), whose position is static and known. By knowing the position of the VRS, the vehicle can be positioned absolutely, too. In the positioning algorithm, we strongly rely on carrier phase positioning which, thanks to its low noise characteristics, may enable decimeter-level positioning. Figure 1 shows the typical geometry of our measurement set-up.

Figure 1: Geometry of the Model including the definition of the attitude of the vehicle

According to figure 1, the two GNSS antennas on the car's rooftop span the array antenna baseline named 𝐛 12 , which one can resolve for the Euler attitude angles to get the vehicle's orientation (heading 𝜓 and pitch 𝜃).

For every time epoch, we estimate both the RTK position and receiver array attitude. The 𝐛 13 , baseline spanned between one vehicle antenna and the antenna of the VRS, instead, enables locate the position of the car relative to the VRS antenna.

𝐛 13 = 𝑥 ⃗ 1 -𝑥 ⃗ 3 (4)

MATHEMATICAL MODELS

The realization of GNSS-PVT navigation is mainly based on a Kalman Filter which is the most popular choice for its optimality and simplicity to be implemented. In this study, an Extended Kalman Filter (EKF) based position and attitude determination algorithm is proposed.

State transition model

The state transition or state-space model describes how the states or parameters of the system vary over time based on a specific linear model.

In our EKF modeling, the state parameter transition between subsequent epochs is given by: 𝐱 𝑛 = 𝛟 𝑛-1 𝐱 𝑛-1 + 𝐰 𝑛 (5) Where:

𝐱 𝑛 refers to the state vector at epoch 𝑛 𝛟 𝑛-1 refers to the transition matrix of epoch 𝑛 -1 𝐰 𝑛 refers to the so-called system noise vector at epoch n Together with the process noise vector, one can define the process noise covariance matrix as:

𝐐 𝑛 = 𝐸[𝐰 𝑛 𝐰 𝑛 T ] (6)
This matrix has then the variances of the state parameter's estimates based on the system model.

The estimated parameters are collected inside the state vector. The state vector collects 5 vehicle state parameters and 2*(𝑁 𝑠𝑎𝑡 -1) satellite state parameters which are: the 3D position of GNSS receiver 1 relative to GNSS receiver 3 (𝐛 13 𝑇 ), the pitch angle of the vehicle (𝜃), the heading angle of the vehicle (𝜓), the double-difference integer ambiguities of the visible satellites seen by GNSS receiver pair 1-3 (𝐍 13 𝑇 ) and the double-difference integer ambiguities of the visible satellites seen by receiver pair 2-3 (𝐍 23 𝑇 ). Note that in a given epoch, we consider only the satellites which are mutually visible by the 3 receivers in this section, to simplify the notations of the models. In our implementation, a different set of satellites can be visible for receiver 2 and 3. The state vector is then given by:

𝐱 = [𝐛

Transition model for satellite-related state parameters

In the case of the satellite-related parameters, they are assumed all as constant over subsequent epochs with a very small noise compared to the position and attitude related state parameters.

𝐍 13,𝑛 = 𝐍 13,𝑛-1 + 𝐰 𝑎𝑚𝑏,𝑛 (9) 𝐍 23,𝑛 = 𝐍 23,𝑛-1 + 𝐰 𝑎𝑚𝑏,𝑛 where:

𝐰 𝑎𝑚𝑏,𝑛 is a centered Gaussian vector with a covariance matrix 𝐐 𝑎𝑚𝑏 = 𝜎 𝑎𝑚𝑏 2 𝐈 𝑁 𝑠𝑎𝑡 -1 𝐈 𝑁 is the identity matrix of size 𝑁

The resulting state transition matrix 𝛟 is then given by an identity matrix and different values of process noise variance are added to complete the model.

𝚽 = [ 1 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 1 ] (10)
And the corresponding process noise matrix 𝐐 is given as follows:

𝐐 = [ 𝐐 13 𝟎 3×1 𝟎 3×1 𝟎 3×(𝑁 𝑠𝑎𝑡 -1) 𝟎 3×(𝑁 𝑠𝑎𝑡 -1) 𝟎 1×3 𝜎 𝜃 2 0 𝟎 1×(𝑁 𝑠𝑎𝑡 -1) 𝟎 1×(𝑁 𝑠𝑎𝑡 -1) 𝟎 1×3 0 𝜎 𝜓 2 𝟎 1×(𝑁 𝑠𝑎𝑡 -1) 𝟎 1×(𝑁 𝑠𝑎𝑡 -1) 𝟎 (𝑁 𝑠𝑎𝑡 -1)×3 𝟎 (𝑁 𝑠𝑎𝑡 -1)×1 𝟎 (𝑁 𝑠𝑎𝑡 -1)×1 𝐐 𝑎𝑚𝑏 𝟎 𝑁 𝑠𝑎𝑡 -1 𝟎 (𝑁 𝑠𝑎𝑡 -1)×3 𝟎 (𝑁 𝑠𝑎𝑡 -1)×1 𝟎 (𝑁 𝑠𝑎𝑡 -1)×1 𝟎 𝑁 𝑠𝑎𝑡 -1 𝐐 𝑎𝑚𝑏 ] (11)

Measurement model

The measurement model describes how the individual sensor measurements are related to system states. In general, for every epoch n, the measurement vector 𝐳 𝒏 , which contains all measured values, can be described as a function of the state vector 𝐱 𝒏 as:

𝐳 𝒏 = 𝐡 𝑛 (𝐱 𝒏 ) + 𝐯 𝑛 ( 12
)
with 𝐡 𝑛 the function that relates one or more states with each measured value and 𝐯 𝑛 the measurement noise vector, which describes the expected Gaussian noise of every measured value. As for the process noise covariance matrix, the definition of the measurement noise covariance matrix follows as:

𝐑 𝑛 = 𝐸[𝐯 𝑛 𝐯 𝑛 T ] (13)
In our model, the measurement vector comprises the following measured values: Double-difference (DD) code phase measurement vector of receiver 1, DD code phase measurement vector of receiver 2, DD carrier phase measurement vector of receiver 1, and DD carrier phase measurement vector of receiver 2.

𝐳 𝑛 = [(𝐏 13 ) 𝑇 (𝐏 23 ) 𝑇 𝜆(𝛟 13 ) 𝑇 𝜆(𝛟 23 ) 𝑇 ] 𝑇 [START_REF] Teunissen | On the Foundation of the Popular Ratio Test for GNSS Ambiguity Resolution[END_REF] In this measurement model, the position of receiver 2 is expressed in terms of the position of receiver 1 and the baseline vector between the 2 receivers of the array, such that it contains the known array baseline length information and the attitude information that we want to estimate. The individual double-difference corrected pseudo-range and phase GPS measurement for our short baseline case (less than 3 km) can be modeled as: where:

𝑃 13 𝑘𝑙 = (
• 𝑃 𝑛𝑝 𝑘𝑙 refers to the double-difference code phase measurement vector of receiver pair 𝑛𝑝 and satellite pair 𝑘𝑙

• 𝜆𝜙 𝑛𝑝 𝑘𝑙 refers to the carrier phase measurement vector of receiver pair 𝑛𝑝 and satellite pair 𝑘𝑙

• 𝐞 𝑘𝑙 refers to the difference between the Line-of-Sight vector of satellite 𝑘 and l towards the receiver

• 𝑁 𝑛𝑝 𝑘𝑙 refers to the double-difference integer ambiguity of receiver pair 𝑛𝑝 and satellite pair 𝑘𝑙

• 𝑛 𝑃,𝑛𝑝 𝑘𝑙 , 𝑛 𝜙,𝑛𝑝 𝑘𝑙 refers to the noise measurement of the double-difference code and phase measurement respectively

To reflect the difference in accuracy between the code measurement and the carrier phase measurement, a fixed scale factor of 𝑎 = 1/100 is applied to account for the much better accuracy of carrier phase measurements.

An elevation-dependent measurement noise variance between all satellites is defined to complete the measurement model, defining the measurement covariance matrix 𝐑. Firstly, the measurement covariance matrix 𝐑 SD for the single difference measurements will have the following shape: ) is a diagonal matrix with the terms 1/ sin(𝑒𝑙 𝑘 ) on its diagonal 𝑒𝑙 𝑘 is the elevation of satellite 𝑘 in radians

𝐑 code,SD = 2 *
The measurement covariance matrix 𝐑 for the double-difference measurements can then be deduced using the following formula:

𝐑 = 𝐃 [ 𝐑 code,SD 𝟎 2 * 𝑁 𝑠𝑎𝑡 𝟎 2 * 𝑁 𝑠𝑎𝑡 𝐑 carrier,SD ] 𝐃 𝑻 (17)
Where

𝑫 = [ -1 0 ⋯ 0 1 0 -1 ⋱ ⋮ 1 ⋮ ⋱ ⋱ 0 ⋮ 0 ⋯ 0 -1 1 ] ( 18
)
is the single-differencing matrix used for computing the double difference. The location of the column full of 1's correspond to the reference satellite for this particular epoch.

The relationship between the state and measurement vector (function 𝐡 𝑛 (𝐱 𝒏 )) is obviously not linear thus we need to linearize this non-linear measurement function and obtain the measurement matrix H (Jacobian matrix of 𝐡 𝑛 (𝐱 𝒏 )) to apply an Extended Kalman Filter. The non-linear measurement function 𝐡 𝑛 (𝐱 𝑛 ) is linearized to the measurement matrix 𝐇 𝑛 as:

𝐡 𝑛 (𝐱 𝑛 ) ≈ 𝐡 𝑛 (𝐱 ̂𝑛 -) + 𝐇 𝑛 (𝐱 𝑛 -𝐱 ̂𝑛 -) (19
) which is then used to calculate the Kalman gain matrix 𝐊 𝑛 and the a-posteriori covariance matrix 𝐏 𝑛 + , as will be seen later in this section. The measurement matrix is computed around the predicted position 𝐱 ̂𝑛as:

𝐇 𝑛 | 𝐱 𝑛 =𝐱 ̂𝑛 -= 𝜕 𝜕𝐱 𝑛 𝐡 𝑛 (𝐱 𝑛 )| 𝐱 𝑛 =𝐱 ̂𝑛 - (20)
So, the corresponding Matrix 𝐇 can be defined as follows: Two alternating steps which are the state prediction step and the state update step are then conducted to complete the proposed EKF algorithm.

𝐇 = [ 𝑒 𝑒 1 𝑒 𝑛 1 𝑒 𝑢 1 0 0 0 … 0 0 … 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ 𝑒 𝑒 𝑁𝑠𝑎𝑡-1 𝑒 𝑛 𝑁𝑠𝑎𝑡-1 𝑒 𝑢 𝑁𝑠𝑎𝑡-1 0 0 0 … 0 0 … 0 𝑒 𝑒 1 𝑒 𝑛 1 𝑒 𝑢 1 ℎ 𝜃 ℎ 𝜓 0 … 0 0 … 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ 𝑒 𝑒 𝑁𝑠𝑎𝑡-1 𝑒 𝑛 𝑁𝑠𝑎𝑡-1 𝑒 𝑢 𝑁𝑠𝑎𝑡-1 ℎ 𝜃 ℎ 𝜓 0 … 0 0 … 0 𝑒 𝑒 1 𝑒 𝑛 1 𝑒 𝑢 1 0 0 𝜆 𝐺𝑃𝑆 … 0 0 … 0 ⋮ ⋮ ⋮ ⋮ ⋮ 0 ⋱ 0 0 ⋱ 0 𝑒 𝑒 𝑁𝑠𝑎𝑡-1 𝑒 𝑛 𝑁𝑠𝑎𝑡-1 𝑒 𝑢 𝑁𝑠𝑎𝑡-1 0 0 0 … 𝜆 𝐺𝑃𝑆 0 … 0 𝑒 𝑒 1 𝑒 𝑛 1 𝑒 𝑢 1 ℎ 𝜃 ℎ 𝜓 0 … 0 𝜆 𝐺𝑃𝑆 0 0 ⋮ ⋮ ⋮ ⋮ ⋮ 0 ⋱ 0 ⋮ ⋱ ⋮ 𝑒 𝑒 𝑁𝑠𝑎𝑡-1 𝑒 𝑛 𝑁𝑠𝑎𝑡-1 𝑒 𝑢 𝑁𝑠𝑎𝑡-1 ℎ 𝜃 ℎ 𝜓 0 … 0 0 0 𝜆 𝐺𝑃𝑆 ] (21

RTK PROCESSING

The cycle slip detection and repair schema based on multi-epoch measurements is first conducted, as introduced in [START_REF] Hu | Cycle-slip Detection and Repair Using an Array of Receivers with Known Geometry for RTK Positioning[END_REF], Cycle slip detection and repair are crucial to maintain continuity of carrier phase observations and to benefit the precise GPS carrier phase observations for high-precision GPS positioning. There are currently many methods to detect and repair the cycle slips. However, it is still a challenging issue, especially in the case of the single-frequency measurements.

In this paper, the Differential Phases of Time Cycle Slip Resolution Method [START_REF] Xu | GPS: theory, algorithms and applications[END_REF] is considered. It is based on the observation of the differential phases between two adjacent epochs, which should exhibit the actual ambiguity, plus some clock errors and remaining noise term as shown in (38):

𝜆𝛥 𝑡 𝛷 = 𝛥 𝑡 𝜌 -𝛥 𝑡 (𝛿𝑡 r -𝛿𝑡 𝑘 )𝑐 + 𝜆𝛥 𝑡 𝑁 + 𝜀 p (22)
Where

• 𝛥 𝑡 𝑥 is the time differential operator applied to parameter 𝑥 • 𝜌 refers to the geometric distance between the satellite and the receiver antennas • 𝛷 refers to the carrier phase measurement • 𝛿𝑡 r , 𝛿𝑡 𝑘 refer to the clock errors of the receiver and satellite respectively • 𝑐 is the speed of light • 𝜀 p refers to remaining errors of the time-differenced phase measurement Except for the ambiguity term, all other terms on the right side are of low variation. Any cycle slips will lead to a sudden jump in the time difference of the phases. Based on the past observation of differential phases measurement, a prediction of the current differenced data can be obtained by polynomial extrapolation or interpolation. The residual between the prediction and the observation can then be used as a detector metric, to be compared to the detection threshold to decide if there are any cycle slips.

After the cycle slips detection process, a cycle-slip validation and size determination process is conducted to verify the determined sizes of cycle-slip. Cycle slips can be repaired using integer vector estimation similarly to ambiguity resolution in the position domain.

After repairing the existing cycle slips, the RTK processing began to be implemented. From the previously described EKF process, we firstly obtain a float estimation of the double-difference integer ambiguity. The accuracy of the position state estimate is further improved by fixing the DD ambiguities to integer numbers by using the well-known LAMBDA [START_REF] Teunissen | The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation[END_REF] [13] algorithm.

One selects the integer candidates based on the sum of squared errors to get a fixed solution. The candidate with the lowest error norm is chosen once the ratio of the Maximum A Posteriori error norm between the second-best candidate and the best candidate is bigger than a threshold. It is a pre-defined threshold or the critical value that the squared norm of ambiguity residuals of the best and second-best candidates should overpass to validate the integer estimation. In our paper, an empirical fixed value of 3 is taken as in [START_REF] Teunissen | On the Foundation of the Popular Ratio Test for GNSS Ambiguity Resolution[END_REF].

Once the IAR process is declared successful, a new position is computed using the DD carrier phase measurements corrected by the validated DD integer ambiguities. This final position is a fixed solution. If the IAR process is not declared successful, the final position is kept as the float solution.

EXPERIMENTAL DATA COLLECT SET-UP AND SCENARIOS

Data collect set-up

In this section, the proposed precise position and attitude determination algorithm is verified with real measurements from two low-cost GNSS receivers and one normal GNSS Receiver. First, the hardware set-up is described. Subsequently, the measurement and process noise will be provided, and the obtained measurement results will be analyzed.

To investigate the feasibility of our proposed precise positioning and attitude determination algorithm, we set up a measurement campaign using 4 low-cost GNSS antennas. Measurements have been taken by recording GPS pseudo-range and carrier phase measurements simultaneously. We took measurements by putting the 4 patch antennas at various distances between each other, from 60 cm up to about 2.0 meters. Besides, we performed the recordings in several different environments including an urban environment, a suburban environment, and an open sky environment.

0] m Table 1 -Summary of data collects

As shown in Table 1, several data collects are processed in this study.

Description of data collection 1 -Open sky environment

To get an open sky and stable environment, the first measurement session took place at the ENAC's football field, as shown in Figure 3, the two receivers are static and their positions are fixed on the football field with favorable satellite visibility, this scenario is mainly used for validation of the algorithm implementation and give a reference of our multi-receiver system performance.

Description of data collection 2 -Suburban environment

The second measurement session took place at the ENAC premises, in Toulouse, France. The true trajectory provided by the Novatel SPAN equipment and the corresponding satellite visibility during the data collection are shown in Figure 4.

Description of data collection 3 -Urban environment

The third measurement session was performed when the vehicle was from ENAC to the city center, in Toulouse, France. The whole trajectory in Google Earth and the corresponding satellite visibility during the data collection are shown in Figure 5.

EXPERIMENTAL RESULTS AND DISCUSSION

Investigation plan

The data collections allow us to investigate a certain number of results that were obtained through simulations in previous papers [START_REF] Hu | Improvement of RTK performances using an array of receivers with known geometry[END_REF], [START_REF] Hu | Cycle-slip Detection and Repair Using an Array of Receivers with Known Geometry for RTK Positioning[END_REF], and some more. Our data collection set-up allows to vary the distance between the 2 rover antennas, named array baseline length in the rest of the paper, as well as the type of environment. In this section, the impact of these 2 parameters (array baseline and type of environment) on the following points will be addressed for the following 3 points:

1. correlation of measurement error collected by an array of receivers 2. improvement of cycle slip detection and repair 3. improvement of positioning accuracy, in terms of horizontal/vertical/attitude accuracy and fixed solution availability

Point 1 can provide some information about the possibility to change the SD observation covariance matrix described in Equation ( 16), by replacing the diagonal matrix with one with non-diagonal terms. By having measurements taken from close antennas, it is assumed that some of the measurement errors are correlated, notably the multipath. The degree of correlation shall depend on the array baseline length. Taking into account the correlation between different measurements will improve the tuning of the EKF and hopefully, the estimation performances.

Point 2 investigates the detection of cycle slips. In order to have a reference about the true occurrence of the cycle slip on our data collects, a post-processing method using dual-frequency measurements [15] will be used. Depending on the observed performance of the cycle slip detection, the process noise of the ambiguity evolution model, involved in the definition of 𝐐 𝑎𝑚𝑏 in equation ( 11) can be optimized.

Finally, using the optimal tuning with regards to points 1 and 2, point 3 is conducted so as to determine the optimal positioning performances using real data.

In conclusion, this investigation plan provides some EKF tuning parameters, and also which array baseline length provides the best positioning performances in all the types of environment tested.

Experimental results on the correlation of measurement errors

As the measurements are coming from signals received by the closely placed antennas, it is safe to consider a certain level of correlation between these noise measurements. For example, the multipath error may be similar to the measurement performed by the 2 closely mounted receivers from the same satellite.

By removing the corresponding geometric distance term from DD code and phase observations, the correlation coefficient between the same DD code and phase observations measured by the different antenna can be computed. Fig. 6 gives an example of the comparison of the code and phase measurement errors for one satellite.

Figure 6: comparison of the code and phase measurement errors for one pair of satellites

Table 2 reports the result of the correlation coefficient of the code and phase measurement errors for different array baseline lengths and various environments. The numbers presented here are the average values of the correlation coefficient of all the common satellites tracked by both the 2 rover receivers. One can notice from the table that the DD phase measurement errors are very correlated and there is also a medium correlation between the code measurement. We speculate that this might be due to the characteristic of the GNSS double-difference measurements. It is important to highlight the fact that there is a number of common parts inside the DD observations of the two closely placed rover receivers, includes the information via the reference station and the reference satellite, especially for the carrier phase measurement. The parts that are not correlated are quite small when comparing to the common parts, which can explain the high correlation relation between the noise measurement of receiver 1 and 2 observed in Table 2. Another finding was that one can observe a law depending on the baseline length: the shorter the distance, the more significant the correlation.

To make our model closer to the real situation thus ameliorate the EKF performance, we then updated the observation covariance matrix in the EKF model by replacing the diagonal matrix with one with non-diagonal terms by adding the correlation coefficient found in Table 2.

Experimental results on the cycle slip detection performance

The cycle slip is mainly due to an event external to the receiver, such as an obstruction or severe multipath environment. In this case, there could be a strong correlation between the occurrence of cycle slips in the 2 receivers who are very close to each other.

That is why at this point, we want to verify if the occurrence of CS is simultaneous on all receivers.

We firstly run a dual-frequency algorithm [15] to obtain true cycle slips. Then compare the number of correct/incorrect CS detection/repair for all single receivers. Figure 7 shows a typical result of the CS detection for one dataset by using these two methods.

Figure 7: CS detection for one dataset by using these two methods.

Table 3 gives the detailed quantitative performance of the cycle slip detection process. The detection rate here describes the ratio of the number of cycle slips successfully detected to the total CS number, while the correct detection rate represents the proportion of correct CS alarms of that method. From the result one can conclude that our proposed single frequency method still has many limitations, only part of the cycle slip can be detected, but the correct correction rate is quite favorable. We also notice the occurrence of several severe CSs is simultaneous on all receivers by both methods, which appear consistent with our previous model thus confirms the hypothesis that we anticipated in [START_REF] Hu | Cycle-slip Detection and Repair Using an Array of Receivers with Known Geometry for RTK Positioning[END_REF]. Figure 9 shows the comparison between the reference trajectory and the dual receiver solution. In some extremely harsh situation, the number of the satellite is not enough to enable a PNT solution for our algorithm, that is why there are some data breakpoints in the graph for the dual receiver solution, while this situation is still improved when compared to single receiver situation thanks to our doubled observations. The dual receiver system performs better than the single receiver situation and provides better solution availability thanks to the doubled observations redundancy.

Experimental results of the RTK performance comparison for different array baseline and various data collections

In our algorithm, for every estimation epoch, if the IAR process is declared successful, a new position is computed using the DD carrier phase measurements corrected by the validated DD integer ambiguities. This position then is used to update the fixed solution. If the IAR process is not declared successful, we keep the float solution.

In the following scenarios, the results are compared between the single receiver system and our proposed multi-receiver system, in terms of the Success fixed-rate, Horizontal positioning (east and north directions) error statistics (mean, standard deviation, and 95% bound), and array attitude error statistics for all 3 data collections. Figure 12 gives an illustration of the Horizontal positioning statistics calculation.

Figure 12: Illustration of Horizontal positioning error statistics calculation

Table 4 gives the RTK performance for our open sky data collect 1, as we can see from the table, the dual receiver array system provides better performance than the single receiver RTK solution, thus demonstrating the interest of such an approach.

Table 4 -Performance comparison for data collection 1 -Open Sky

In the second step, we want to analyze the robustness against the harsh environment of multi-receiver RTK systems. The data used for this study was collected in Toulouse suburban area (dataset 2) and urban area (dataset 3).

The sub-urban data was collected when the vehicle was driven at the ENAC campus. The reference trajectory was provided with cm-level accuracy. The maximum standard deviation values, up to 10 cm, occur around 500 epochs, which correspond also to the zone having minimum visible satellites. The environment is quite favorable with at least 8 satellites in view for most of the time.

Table 5 gives the performance comparison between the proposed dual-receiver system and the single receiver system in the suburban environments and for different array baseline lengths. Again, a better accuracy result is obtained in the dual receiver situation. Moreover, there is no huge difference in the accuracy of the positioning results for all the dual receivers' situations. however, as expected, the attitude accuracy does improve sightly as the array baseline increases. The urban area data were collected around the city center. A dm-level trajectory accuracy was obtained by the SPAN receiver.

The number of visible satellites decreased a lot. A clear uncertainty increase in the position solution was observed in the reference result during the section where the number of satellites was less than 6. Table 6 presents the performance comparison between the proposed dual-receiver system and the single receiver system in the urban environments and for different array baseline length.

Table 6 -Performance comparison for different array baseline for data collection 3 -Urban

As may be seen from the table, the planned comparisons revealed that the use of the receiver array improves the fixed rate, positioning accuracy, and also the solution availability for all considered array baseline values, this allows the conclusion that the use of an array of receivers with known geometry to improve the RTK performance is feasible and effective.

CONCLUSION

This paper is the extension and improvement of our previous work [START_REF] Hu | Improvement of RTK performances using an array of receivers with known geometry[END_REF], [START_REF] Hu | Cycle-slip Detection and Repair Using an Array of Receivers with Known Geometry for RTK Positioning[END_REF]. In this contribution, we present a method that includes an array of receivers with known geometry to enable the vehicle attitude determination and enhance the RTK performance in different environments. Taking advantage of the attitude information and the known geometry of the array of receivers, we are able to improve some internal steps of precise position computation.

In order to optimize our processing, the correlation of the measurement errors affecting observations taken by our array of receivers has been determined. Then, the performance of our real-time single frequency cycle-slip detection and repair algorithm has been assessed, by comparing it to the cycle slips detected by a dual-frequency algorithm performed in post-processing. These two investigations yielded important information so as to tune our Kalman Filter. Additionally, our method achieves a relatively accurate estimation of the attitude of the vehicle which provides additional information beyond the positioning.

Finally, we demonstrate through real data processing results that our multi-receiver RTK system is more robust to degraded satellite geometry, in terms of ambiguity fixing rate, and get a better position accuracy under the same conditions when comparing with the single receiver system. Our experiments correlate favorably with our previous simulation results and further support the idea of using an array of receivers with known geometry to improve the RTK performance.

  ) with ℎ 𝜃 = (𝑒 𝑘𝑙 ) 𝑇 * (-𝑙) * [ -sin𝜃sin𝜑 -sin𝜃cos𝜑 cos𝜃 ] and ℎ 𝜓 = (𝑒 𝑘𝑙 ) 𝑇 * (-𝑙) * [ 𝑐𝑜𝑠𝜃cos𝜑 -𝑐𝑜𝑠𝜃sin𝜑 0] respectively.

Figure 2 :Figure 2

 22 Figure 2: Real data collection set-up: 4 GNSS u-blox antennas and one Novatel SPAN Receiver on the car rooftop Figure 2 shows a typical set-up of the GNSS antennas for one of the measurement sessions. According to Figure 1, the two GNSS antennas on the car's rooftop span the 𝐛 12 vehicle antenna baseline, which one can resolve for the pitch and heading attitude angles to get the vehicle's orientation. The 𝐛 13 baseline spanned between one vehicle antenna and the VRS antenna, instead, is able to position the car relative to the VRS antenna. If the VRS antenna location is known, the absolute position of the car can be determined. The measurement test was performed with a vehicle on which the following hardware was mounted: • 4 u-blox F9P GPS Receivers with 1 Hz data rate • 4 L1 patch antennas mounted on the roof of a vehicle along its longitudinal axis • 1 Septentrio AsteRx-U Receiver on the roof of a building as the reference station for RTK processing • 1 Novatel Propak6 SPAN Receiver tightly-coupled with tactical grade IMU, on multi-baseline post-processing RTK to give the reference position and attitude of the vehicle

Figure 3 :

 3 Figure 3: Receiver fixed position for dataset 1

Figure 4 :

 4 Figure 4: Trajectory and Corresponding Satellite visibility for dataset 2

Figure 5 :

 5 Figure 5: Trajectory and Corresponding Satellite visibility for dataset 3

  Due to the imperfect detection of the cycle slip, a bigger value 𝜎 𝑎𝑚𝑏 = 0.1 was chosen to account for possible undetected cycle slips. Slight improvement in positioning results is achieved when comparing with a smaller value 𝜎 𝑎𝑚𝑏 = 0.01.Typical Experimental results of the positioning and attitude determination performanceTypical results of the state estimations (positioning and attitude of the vehicle) for the different data collection cases are given in the following figures. In figure8, the 2 rows show the estimation of the horizontal (East and North) positioning coordinate respectively with respect to the true position of the receiver.

Figure 8 :

 8 Figure 8: Illustration of positioning estimation for dataset 3

Figure 9 :

 9 Figure 9: Comparison between the reference trajectory and Dual receiver solutionFigure10gives the illustration of the 3D RTK Positioning estimation error result for our multi-receiver method, as one can notice, the algorithm succeeded in outputting the right positioning result, with an accuracy of about 0.5 m for all three directions, which is acceptable for a harsh environment.

Figure 10 :

 10 Figure 10: 3D RTK Positioning Estimation error illustrationFigure11shows the estimation of the pitch and heading angles of the vehicle for dataset 1. One can see from the figure that the error between the estimated result and the true value is extremely small (less than 2 deg), which can provide us with a relatively accurate vehicle posture by using our proposed method.

Figure 11 :

 11 Figure 11: Illustration of vehicle attitude estimation for dataset 1

  The array antenna baseline vector 𝐛 12 has the following coordinates in the local navigation (East, North, Up) ENU frame:

				cos(𝜃) sin(𝜓)
	𝐛 12 (𝜓, 𝜃) = ‖𝐛 12 ‖. [	cos(𝜃) cos(𝜓)	]	(1)
				sin(𝜃)
	Solving this equation for the attitude angles heading 𝜓 and the pitch 𝜃 yields:
	𝜓(𝑡 𝑖 ) = atan (	(𝐛 12 ) 𝑒 (𝑡 𝑖 ) (𝐛 12 ) 𝑛 (𝑡 𝑖 )	)	(2)
	𝜃(𝑡 𝑖 ) = atan (	(𝐛 12 ) 𝑢 (𝑡 𝑖 ) √(𝐛 12 ) 𝑒 2 (𝑡 𝑖 ) + (𝐛 12 ) 𝑛 2 (𝑡 𝑖 )	)

Transition model for position and attitude related state parameters In

  our EKF modeling, for the position and attitude related state parameters, we suppose that they follow a random walk model, meaning that the speed and the angular rate are a zero-mean Gaussian process.

	𝐛 13,𝑛 = 𝐛 13,𝑛-1 + 𝐰 13,𝑛
	𝜃 𝑛 = 𝜃 𝑛-1 + 𝑤 𝜃,𝑛		(8)
	𝜓 𝑛 = 𝜓 𝑛-1 + 𝑤 𝜓,𝑛
	where:
	𝐰 13 is a centered Gaussian vector with a covariance matrix 𝐐 13
	𝑤 𝜃 is a centered Gaussian variable with a standard deviation 𝜎 𝜃
	𝑤 𝜓 is a centered Gaussian variable with a standard deviation 𝜎 𝜓
	13 𝑇	𝜃 𝜓 𝐍 13 𝑇	𝐍 23 𝑇 ] 𝑇	(7)

  𝐞 𝑘𝑙 ) 𝑇 𝐛 13 + 𝑛 𝑃,13

	𝑘𝑙
	cos 𝜃 sin 𝜓
	𝑃 23 𝑘𝑙 = (𝐞 𝑘𝑙 ) 𝑇 (𝐛 13 -|𝐛 12 | [ cos 𝜃 cos 𝜓	]) + 𝑛 𝑃,23 𝑘𝑙	(15)
	sin 𝜃
	𝜆𝜙 13 𝑘𝑙 = (𝐞 𝑘𝑙 ) 𝑇 𝐛 13 + 𝜆𝑁 13 𝑘𝑙 + 𝑛 𝜙,13 𝑘𝑙
	𝜆𝜙 23 𝑘𝑙 = (𝐞 𝑘𝑙 ) 𝑇 (𝐛 13 -|𝐛 12 | [ cos 𝜃 cos 𝜓 cos 𝜃 sin 𝜓	]) + 𝜆𝑁 23 𝑘𝑙 + 𝑛 𝜙,23 𝑘𝑙
	sin 𝜃

Table 2 -Correlation coefficient of the code and phase measurement errors

 2 Array baseline length data collect #1 Open Sky data collect #2 Sub-urban data collect #3 Urban 0.6 m

		0.576	0.999	0.538	0.999	0.475	0.999
	1.3 m	0.568	0.986	0.499	0.997	0.392	0.998
	2.0 m	0.545	0.963	0.486	0.996	0.359	0.997

Table 3 -Cycle slip detection performances

 3

		data collect #1 Open Sky	data collect #2 Sub-urban	data collect #3 Urban
		Detection	Correct	Detection	Correct	Detection	Correct
		Rate	Detection Rate	Rate	Detection Rate	Rate	Detection Rate
	Reference algorithm	1	1	1	1	1	1
	Receiver 1	0.38	0.84	0.36	0.86	0.31	0.86
	Receiver 2	0.39	0.89	0.38	0.85	0.26	0.85
	Receiver 3	0.36	0.83	0.27	0.84	0.28	0.84
	Receiver 4	0.37	0.86	0.32	0.79	0.29	0.75

Table 5 -Performance comparison for different array baseline for data collect 2 -Suburban

 5

	Variation of scenario 1	Performance indicator	Unit	Mean	Standard deviation	95% bound
	Single receiver	Horizontal positioning error Success fixed rate	m %	0.2880 69.35	0.1438 N/A	0.5301 N/A
	Dual receiver	Horizontal positioning error Absolute heading error	m deg	0.2185 0.4666	0.1240 0.8054	0.4518 0.8971
	𝑙 = 2 m	Absolute pitch error Success fixed rate	deg %	0.3524 84.69	0.7755 N/A	0.9355 N/A

ACKNOWLEDGEMENTS

This work has been supported by the China Scholarship Council (CSC). It is the Chinese Ministry of Education's non-profit organization that provides support for international academic exchange with China.

Code Phase Code Phase Code Phase


Véronique Roig-Zamboni

Beatrice Cobucci-Ponzano

Roberta Iacono

Maria Carmina Ferrara

Stanley Germany

Yves Bourne

Giancarlo Parenti

Marco Moracci

Gerlind Sulzenbacher
email: gerlind.sulzenbacher@afmb.univ-mrs.fr

Structure of human lysosomal acid α-glucosidase-a guide for the treatment of Pompe disease

est destinée au dépôt

P

ompe disease (also known as glycogen storage disease type 2 or acid maltase deficiency) is an autosomal recessive disorder caused by mutations in the gene that encodes the hydrolase acid α-glucosidase (GAA), member of glycoside hydrolase family GH31 [START_REF] Lombard | The carbohydrate-active enzymes database (CAZy) in 2013[END_REF] , and involved in the lysosomal breakdown of glycogen. Functional deficiency of GAA results in lysosomal accumulation of glycogen and cellular damage in all tissues, particularly cardiac and skeletal muscle [START_REF] Van Der Ploeg | Pompe's disease[END_REF][START_REF] Raben | Acid alpha-glucosidase deficiency (glycogenosis type II, Pompe disease)[END_REF] . Pompe disease is characterized by a broad phenotypic spectrum that ranges from a slowly progressing late-onset phenotype to a devastating classical infantile-onset, but in all cases, progressive muscle hypotonia and loss of motor, respiratory and cardiac functions leads to respiratory failure [START_REF] Van Der Ploeg | Pompe's disease[END_REF][START_REF] Parenti | Pompe disease: from new views on pathophysiology to innovative therapeutic strategies[END_REF] . Recombinant human GAA (rhGAA) produced in CHO cells has been approved in 2006 for enzyme replacement therapy (ERT) to treat Pompe disease and has proven to be beneficial for patient's survival and to stabilize the disease course [START_REF] Van Den Hout | Long-term intravenous treatment of Pompe disease with recombinant human alpha-glucosidase from milk[END_REF][START_REF] Kishnani | Chinese hamster ovary cell-derived recombinant human acid alpha-glucosidase in infantile-onset Pompe disease[END_REF][START_REF] Chen | Reversal of cardiac dysfunction after enzyme replacement in patients with infantile-onset Pompe disease[END_REF][START_REF] Rossi | Long-term enzyme replacement therapy for pompe disease with recombinant human alpha-glucosidase derived from chinese hamster ovary cells[END_REF] . Since then, ERT is the only approved treatment for Pompe disease, but despite the clinical benefits, patient's response is very variable and the efficiency of the treatment is limited by insufficient targeting and uptake in muscle tissues, immunogenic reactions and build-up of autophagic compartments in myocytes [START_REF] Van Der Ploeg | A randomized study of alglucosidase alfa in late-onset Pompe's disease[END_REF][START_REF] Cardone | Abnormal mannose-6-phosphate receptor trafficking impairs recombinant alpha-glucosidase uptake in Pompe disease fibroblasts[END_REF][START_REF] Kishnani | Cross-reactive immunologic material status affects treatment outcomes in Pompe disease infants[END_REF][START_REF] Fukuda | Autophagy and mistargeting of therapeutic enzyme in skeletal muscle in Pompe disease[END_REF] . Hence, the quest for alternative therapeutic strategies based on different approaches and rationale has become compulsory and pharmacological chaperone therapy (PCT) has been proposed as an alternative or complementary approach to ERT [START_REF] Parenti | Pompe disease: from new views on pathophysiology to innovative therapeutic strategies[END_REF][START_REF] Parenti | Treating lysosomal storage diseases with pharmacological chaperones: from concept to clinics[END_REF][START_REF] Parenti | Pharmacological chaperone therapy: preclinical development, clinical translation, and prospects for the treatment of lysosomal storage disorders[END_REF] . PCT is based on the concept that smallmolecule ligands may act by blocking conformational fluctuations of a partially misfolded protein, rescuing its functional state and allowing escape from the endoplasmic reticulum-associated protein degradation (ERAD) machinery (Fig. 1) [START_REF] Fan | A counterintuitive approach to treat enzyme deficiencies: use of enzyme inhibitors for restoring mutant enzyme activity[END_REF] . Hundreds of disease-causing mutations in the gaa gene have been identified, including insertions, deletions, splice site, nonsense, and missense mutations (Supplementary Table 1). Missense mutations result in production of full-length GAA, likely to not fold as efficiently as the stable wild-type enzyme, and patients affected by these mutations are thus potential candidates for PCT [START_REF] Flanagan | The pharmacological chaperone 1-deoxynojirimycin increases the activity and lysosomal trafficking of multiple mutant forms of acid alpha-glucosidase[END_REF] . While not more than ~10-15% patients are estimated to be amenable to PCT [START_REF] Flanagan | The pharmacological chaperone 1-deoxynojirimycin increases the activity and lysosomal trafficking of multiple mutant forms of acid alpha-glucosidase[END_REF] , it has been shown that pharmacological chaperones act as enzyme enhancers when co-administered with rhGAA, by favoring enzyme delivery, stability and maturation [START_REF] Porto | The pharmacological chaperone N-butyldeoxynojirimycin enhances enzyme replacement therapy in Pompe disease fibroblasts[END_REF] , making this PCT independent from the type of mutation carried by patients.

Here we report the high-resolution structures of mature rhGAA and its complexes with the GAA inhibitors acarbose, 1-deoxynojirimycin (DNJ), N-hydroxyethyl-DNJ (NHE-DNJ) and with the allosteric pharmacological chaperone N-acetylcysteine [START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF] . These structures give insight into substrate recognition, support at the molecular level the action of currently known pharmacological chaperones, and provide a molecular framework for the rationalization of mutations in clinical isolates of individuals affected by Pompe disease.

Results

Structural overview. GAA is synthesized as a 110 kDa glycoprotein, which is targeted to the lysosome via the mannose-6-phosphate receptor and undergoes in the late endosomal/ lysosomal compartment a series of proteolytic and N-glycan processing events to yield a mature active form composed of four tightly associated peptides [START_REF] Van Hove | High-level production of recombinant human lysosomal acid alpha-glucosidase in Chinese hamster ovary cells which targets to heart muscle and corrects glycogen accumulation in fibroblasts from patients with Pompe disease[END_REF][START_REF] Moreland | Lysosomal acid alpha-glucosidase consists of four different peptides processed from a single chain precursor[END_REF] . Because crystallization of the commercial Myozyme ® precursor form of rhGAA (Q57-C952) did not afford crystals diffracting beyond ~7 Å and protein disorder predictors [START_REF] Linding | Protein disorder prediction: implications for structural proteomics[END_REF] indicated the presence of disordered peptide regions, we performed in situ proteolysis with α-chymotrypsin to remove putative flexible surface loops hampering formation of productive crystal contacts. The proteolytic treatment yielded a polypeptide of ~5 kDa lower mass than the rhGAA precursor (Fig. 2a), which crystallized readily. This allowed us to collect the diffraction data extending to 1.9 Å and solve the structure of rhGAA by molecular replacement (Table 1). The proteolytically digested form had activity comparable to that of the precursor (2.34 ± 0.06 and 2.26 ± 0.16 U mg -1 for the precursor and mature The crystal structure of rhGAA corresponds to the mature form of human placental GAA and rhGAA [START_REF] Moreland | Lysosomal acid alpha-glucosidase consists of four different peptides processed from a single chain precursor[END_REF] (Fig. 2b), suggesting that treatment with α-chymotrypsin afforded removal of protease-labile regions (Q57-T80, G116-M122, R199-A204 and A782-R794) in a manner comparable to the proteolytic maturation occurring in vivo. The overall architecture of rhGAA is similar to that of the human glycoside hydrolase family GH31 1 homologues, the intestinal brush-border enzymes maltaseglucoamylase (MGAM) [START_REF] Sim | Human intestinal maltase-glucoamylase: crystal structure of the N-terminal catalytic subunit and basis of inhibition and substrate specificity[END_REF][START_REF] Ren | Structural insight into substrate specificity of human intestinal maltase-glucoamylase[END_REF] and sucrase-isomaltase (SI) [START_REF] Sim | Structural basis for substrate selectivity in human maltaseglucoamylase and sucrase-isomaltase N-terminal domains[END_REF] (Fig. 2c and Supplementary Fig. 1). A N-terminal trefoil Type-P domain is followed by a β-sheet domain, the catalytic (β/α) 8 barrel bearing two inserts after β-strands β3 (insert I) and β4 (insert II), and proximal and distal β-sheet domains at the C-terminus. The precursor form of rhGAA contains ~14 kDa of carbohydrate and seven glycosylation sites have been predicted and characterized [START_REF] Moreland | Lysosomal acid alpha-glucosidase consists of four different peptides processed from a single chain precursor[END_REF] . We observed glycan structures of various lengths for five of them in the crystal structure, notably at N140, N233, N390, N470 and N652 (Fig. 2c and Supplementary Fig. 2). We detected residual difference electron density at N882, insufficient to model a glycan structure, and did not observe difference electron density near N925.

Active site and inhibitor binding. The narrow substrate-binding pocket of rhGAA is located near the C-terminal ends of β-strands of the catalytic (β/α) 8 domain and shaped by a loop from the N-terminal β-sheet domain and both inserts I and II (Fig. 2c). Family GH31 enzymes perform catalysis via a classical Koshland double displacement reaction mechanism with retention of the anomeric carbon configuration in the product. From early characterization of the catalytic site of GAA [START_REF] Hermans | Human lysosomal alpha-glucosidase. characterization of the catalytic site[END_REF] and by homology with mechanistically dissected members of family GH31 [START_REF] Lovering | Mechanistic and structural analysis of a family 31 alpha-glycosidase and its glycosyl-enzyme intermediate[END_REF][START_REF] Kim | Expanding the thioglycoligase strategy to the synthesis of alpha-linked thioglycosides allows structural investigation of the parent enzyme/substrate complex[END_REF] , the catalytic nucleophile and acid/base can be assigned to D518 and D616, respectively. To get insight into the interaction of rhGAA with the documented pharmacological chaperone 1-deoxynojirimycin [START_REF] Flanagan | The pharmacological chaperone 1-deoxynojirimycin increases the activity and lysosomal trafficking of multiple mutant forms of acid alpha-glucosidase[END_REF][START_REF] Kishnani | Duvoglustat HCl increases systemic and tissue exposure of active acid α-glucosidase in pompe patients co-administered with alglucosidase α[END_REF] and its derivative N-hydroxyethyldeoxynojirimycin [START_REF] Bruckmann | Systematic structure-activity study on potential chaperone lead compounds for acid α-glucosidase[END_REF] (DNJ and NHE-DNJ, respectively), we soaked rhGAA crystals with these active site directed iminosugar inhibitors and obtained for both complexes diffraction data extending to 2.0 Å resolution (Table 1). DNJ binds in a 4 C 1 conformation in the substrate-binding subsite -1 [START_REF] Davies | Nomenclature for sugar-binding subsites in glycosyl hydrolases[END_REF] and is stabilized by hydrogen bonds to the side-chains of D404, D518, R600, D616, and H674. Further stabilizing interactions are provided by hydrophobic contacts with W376, I441, W516, M519, W613, and F649 (Fig. 3a,b). Apart from a 20°tilt in the side-chain of W376, no major conformational changes occur upon DNJ binding in the rhGAA active site. This is not the case for binding of NHE-DNJ, which adopts the same pose as DNJ, but where the hydroxyethyl substituent induces substantial conformational changes to the side chains of M519 and W481 (Fig. 3c,d). In this ligand-induced conformation, rhGAA residue W481, located on the tip of insert I, establishes a stabilizing interaction with the hydroxyethyl substituent of NHE-DNJ. Similar conformational changes can be observed upon binding of NHE-DNJ to the N-terminal subunit of MGAM (NtMGAM) [START_REF] Sim | New glucosidase inhibitors from an ayurvedic herbal treatment for type 2 diabetes: structures and inhibition of human intestinal maltaseglucoamylase with compounds from Salacia reticulata[END_REF] (Supplementary Fig. 3a).

Substrate recognition and specificity. To gather a perception of substrate recognition by rhGAA, we acquired the diffraction data extending to 2.45 Å resolution for a rhGAA crystal soaked with acarbose, a non-cleavable α-1,4-tetrasaccharide substrate analogue (Table 1). In this complex, the valienamine ring lodged in subsite -1 adopts a 2 H 3 half-chair conformation, but despite this conformational difference with respect to DNJ bound in the 4 C 1 half-chair conformation, the hydrogen-bonding pattern in subsite -1 is essentially identical for the two compounds (Fig. 3e, f and Supplementary Fig. 3b). The non-hydrolysable "interglycosidic" nitrogen occupies the catalytic centre and is hydrogen bonded to the catalytic acid/base D616. The 6-deoxyglucosyl moiety in subsite + 1 establishes via its 2-and 3-hydroxyl groups hydrogen bonds with the conserved R600 and with D282 originating from a loop in the N-terminal β-sheet domain. Though not conserved in sequence, residues of this loop interact invariably with carbohydrate ligands in all the available crystal structures of family GH31 members. The maltose unit of acarbose in subsites + 2 and + 3 does not make any direct interactions with rhGAA and is rather stabilized by crystal lattice packing interactions and a water mediated contact with the side chain of W618, located at the rim of the substrate-binding pocket. This paucity of interactions suggests that rhGAA possesses only two productive substrate-binding sites, subsites -1 and + 1, but interestingly, the acarbose maltose moiety bound to rhGAA adopts the same overall pose as when bound to NtMGAM [START_REF] Sim | Human intestinal maltase-glucoamylase: crystal structure of the N-terminal catalytic subunit and basis of inhibition and substrate specificity[END_REF] suggesting a functional role also for subsites + 2 and + 3 in rhGAA (Supplementary Fig. 3c).

Glycogen is a large polymer build-up by linear chains of α-1,4linked glucose residues carrying α-1,6-linked branches. In the cytosol the energy storage particle is degraded by the concomitant action of glycogen phosphorylase and a debranching enzyme. Within the lysosome, no debranching enzyme is present and GAA must assure the hydrolysis of both α-1,4and α-1,6glycosidic linkages. However, rhGAA shows clear preference for the former linkage, as the specificity constant on maltose is 32-fold higher when compared to the specificity constant on isomaltose (Supplementary Table 2). To uncover the structural basis for the dual substrate specificity, we attempted soaking of a α-1,6-linkage between subsite -1 and + 1 and that productive hydrogen bonds between the glycopyranose moiety in subsite + 1 and D282 are maintained (Fig. 4a,b). However, due to the increased length of the α-1,6-linkage with respect to an α-1,4linkage, the stabilizing hydrogen bond with R600 is weakened, which might explain lower activity of the enzyme on α-1,6-linked branches.

The secondary substrate-binding domain. In the rhGAAacarbose complex, we could observe the acarvosine moiety of a second acarbose molecule lodged in a pocket within the N-terminal trefoil Type-P domain, situated on the same face of rhGAA where the active site is located and 25 Å away from acarbose bound therein (Fig. 4c). The valienamine moiety hydrogen-bonds with its 4-and 6-hydroxyl groups to the mainchain nitrogen and oxygen of C127, the main-chain oxygen of A93 and the side-chain of D91. Further stabilizing interactions are provided by stacking interactions with P125 and W126. P125 and D91 are conserved or replaced by conservative substitutions in the family GH31 enzymes comprising a trefoil Type-P domain (Supplementary Fig. 4a). Adjacent to the secondary carbohydratebinding site the surface loop G116-M122 had been trimmed by proteolytic treatment. Conceivably, the removal of the segment G116-M122, unique to the lysosomal representatives of family GH31 (Supplementary Fig. 4a,b), might contribute to the formation of the secondary substrate-binding pocket. One could envisage that the secondary acarbose-binding site on rhGAA represents a genuine substrate-binding site possibly enhancing the processivity of the enzyme. This site could also represent a pharmacophore for in silico screening for novel pharmacological chaperones. Most notably, a trisaccharide derived from acarbose is bound to the N-terminal domain of family GH31 α-1,4-glucan lyase from Gracilariopsis lemaneiformis 33 in a position close to the secondary substrate-binding site described here (Supplementary Fig. 5).

N-acetylcysteine is an allosteric pharmacological chaperone.

The known pharmaceutical drug N-acetylcysteine (NAC) has been shown to act as allosteric pharmacological chaperone enhancing the stability of mutated endogenous GAA and rhGAA used for ERT, without affecting the catalytic activity [START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF] . To harness the therapeutic potential of NAC for Pompe disease via structural studies, we pursued the crystal structure of rhGAA in complex with NAC. The 1.83 Å resolution structure of the complex obtained by crystal soaking experiments (Table 1), reveals two NAC-binding sites remote from the active site (Fig. 4d,e), providing structural evidence that NAC is indeed an allosteric chaperone, as anticipated by functional studies [START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF] . We ascertained via activity assays and differential scanning fluorimetry that NAC stabilizes equally well rhGAA and the proteolytically matured enzyme produced in this study, and that the chaperone is specific towards GAA, as NAC has no stabilizing effect on a GH31 homologue from rice (Supplementary Fig. 6a,b and Supplementary Tables 3 and4). In the rhGAA-NAC complex structure a NAC molecule, designated as NAC1, is located about 30 Å away from the active site, at the interface between the (β/α) 8 barrel and the distal β-sheet domain (Fig. 4d). NAC1 binds to the main chain carboxyl of H432 via its amide nitrogen and the carboxyl function makes a water mediated contact with the side chain of D513. The Cβ-Sγ and the acetyl groups establish stacking interactions with the guanidine group of R437 and the loop G434-G435, respectively. A second, partially (75%) occupied NAC molecule (NAC2), is lodged at the interface between the (β/α) 8 barrel and the proximal β-sheet domain at a distance of 25 Å from the active site (Fig. 4e). The carboxyl function of NAC2 makes a water-mediated contact with the main-chain oxygen of L756, the thiol stacks against the side chain of Q757 and the acetyl group stacks against the side chain of H742. The two NAC molecules establish fewer and weaker interactions with rhGAA than the chaperones DNJ and NHE-DNJ. This reflects the higher concentrations of NAC needed for chaperone activity (mM vs, µM) as exemplified by the K D of 11.57±0.74 mM that we obtained by differential scanning fluorimetry and by the K i of 3.4 and 3.0 µM for DNJ and NHE-DNJ, respectively [START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF][START_REF] Bruckmann | Systematic structure-activity study on potential chaperone lead compounds for acid α-glucosidase[END_REF] . A sigmoidal saturation curve best fitting the experimental points supports the NAC allosteric fully (NAC1) and partially (NAC2) occupied binding sites (Supplementary Fig. 6c). The stacking interactions established between the acetyl groups of NAC1 and NAC2 and rhGAA must contribute significantly to the binding energy, since amino acids related to NAC, such as N-acetylserine and N-acetylglycine have equally a stabilizing effect on rhGAA, whilst the non-acetylated counterparts have no effect [START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF] . The structure of the rhGAA-NAC complex shows a decrease in atomic thermal displacement parameters of residues surrounding the NACbinding sites with respect to that of unbound rhGAA, pointing toward an overall inter-domain stabilizing function driven by NAC (Supplementary Fig. 7a-d). In the crystal NAC appears also to exercise an anti-oxidative effect as residue C938 is oxidized to the sulfenic acid form in all the crystal structures described here, except for the structure of the rhGAA-NAC complex (Supplementary Fig. 8a,b).

NAC and iminosugars show different chaperoning profiles. Some mutations leading to Pompe disease are responsive to both NAC and the iminosugars DNJ and NB-DNJ, while others are targeted specifically by one or the other pharmacological chaperone [START_REF] Flanagan | The pharmacological chaperone 1-deoxynojirimycin increases the activity and lysosomal trafficking of multiple mutant forms of acid alpha-glucosidase[END_REF][START_REF] Porto | The pharmacological chaperone N-butyldeoxynojirimycin enhances enzyme replacement therapy in Pompe disease fibroblasts[END_REF][START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF] . The vast majority of GAA mutants responsive to the chaperone activity of iminosugars are located either in reach of the active site or on structural elements contributing to its global architecture. In Pompe patient, fibroblasts and mouse models the GAA mutants most responsive to NAC (A445P, Y455F, and L552P) [START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF] are located on inserts I and II, where the corresponding wild-type residues contribute to domain stabilizations. A445 together with F487 defines the boundaries of insert I, and the A445P mutant most likely destabilizes the entire scaffold of insert I (Supplementary Fig. 9a). The side-chain hydroxyl of Y455 on insert I interacts with a long loop of the catalytic domain whilst the side chain of L552 on insert II makes hydrophobic interactions with residues from the N-terminal β-sheet domain (Supplementary Fig. 9b,c). The respective Y455F and L552P mutants have most certainly an overall destabilizing effect, which is clearly rescued by NAC. Hence, the chaperone action of iminosugars can be considered as a conformational rescue of structural perturbed active sites, whereas the effect of the allosteric chaperone NAC appears to be due to the stabilization of overall or local conformational fluctuations.

Discussion

Although lysosomal storage diseases have low incidence in overall population, they represent a heavy burden to patients' health and have an important socio-economic impact. ERT with recombinant human rhGAA has become available as drug for the standard care for Pompe disease a decade ago and has proven to be beneficial for patient's survival and to stabilize the disease course [START_REF] Van Den Hout | Long-term intravenous treatment of Pompe disease with recombinant human alpha-glucosidase from milk[END_REF][START_REF] Kishnani | Chinese hamster ovary cell-derived recombinant human acid alpha-glucosidase in infantile-onset Pompe disease[END_REF][START_REF] Chen | Reversal of cardiac dysfunction after enzyme replacement in patients with infantile-onset Pompe disease[END_REF][START_REF] Rossi | Long-term enzyme replacement therapy for pompe disease with recombinant human alpha-glucosidase derived from chinese hamster ovary cells[END_REF] . Undoubtedly, structural information was desired to harness the therapeutic value of rhGAA, but remarkably, this live-saving drug could not be crystallized so far and its 3D-structure remained elusive. Yet, shortly after submission of this manuscript the Garman group released two isomorphous structures of recombinantly expressed human GAA in complex with glucose at the Protein Data Bank, with accession numbers 5KZW and 5KZX.

In this study we employed in situ proteolysis to afford well-diffracting crystals of rhGAA and interestingly the crystallographic model we obtained corresponds to the mature forms of recombinant and placental GAA found in the endosome/lysosomal compartment [START_REF] Moreland | Lysosomal acid alpha-glucosidase consists of four different peptides processed from a single chain precursor[END_REF] . One of the peptide segments removed by proteolytic treatment (G116-M122) is specific to the lysosmal representatives of the GH31 enzyme family and its removal might be required for the generation of the secondary substrate-binding | DOI: 10.1038/s41467-017-01263-3 | www.nature.com/naturecommunications site in the trefoil Type-P domain. De facto, none of the over a dozen native or complex structures of the human intestinal homologues MGAM and SI, where the peptide range corresponding to rhGAA G116-M122 is present, contains a carbohydrate or carbohydrate-like ligand such as glycerol in the area corresponding to the rhGAA secondary substrate-binding site. The secondary substrate-binding site might favor the adherence of the enzyme onto the glycogen particle during hydrolytic processing and enhance its ability to detangle the α-1,6-linked branch points, which in cytosolic glycogen metabolism are taken care of by a debranching enzyme, not present in the lysosome. Except for the human intestinal homologues, the other structurally characterized GH31 members do not contain a trefoil Type-P domain, but interestingly a secondary substrate-binding site is equally present in the N-terminal region of G. lemaneiformis α-1,4-glucan lyase, which degrades starch following an elimination mechanism instead of hydrolysis. With an increasing amount of forthcoming structural and functional data it should be interesting to see whether within glycoside hydrolase family GH31 relationships can be established between particular modes of substrate binding and degradation. Despite the clinical benefits of ERT with rhGAA, patient's response is very variable and mostly limited by insufficient enzyme targeting to the lysosomal compartment, high cost and heavy burden of treatment. In recent years PCT, presenting advantages such as low cost and burden, and high distribution volume, has emerged as an attractive alternative, either alone or in combination with ERT, for the treatment of LSDs [START_REF] Parenti | Pharmacological chaperone therapy: preclinical development, clinical translation, and prospects for the treatment of lysosomal storage disorders[END_REF][START_REF] Parenti | Pharmacological chaperone therapy for lysosomal storage diseases[END_REF] . For Pompe disease, for which clinical trials combining ERT and treatment with NB-DNJ and DNJ have recently been concluded [START_REF] Kishnani | Duvoglustat HCl increases systemic and tissue exposure of active acid α-glucosidase in pompe patients co-administered with alglucosidase α[END_REF][START_REF] Parenti | A chaperone enhances blood α-glucosidase activity in Pompe disease patients treated with enzyme replacement therapy[END_REF] , PCT is particularly promising as pharmacological chaperones are small molecules that are expected to penetrate membranes and reach therapeutic levels in target tissues, including skeletal muscle (38-55% of human body weight), more easily than rhGAA. A potential drawback of most pharmacological chaperones under study for the rescue of lysosomal glycosidases is that they are active site-directed iminosugars and inhibitors of the target enzymes. Furthermore, glucose configured iminosugars targeting GAA have the additional drawback of inhibiting also the intestinal starch-digesting enzymes and the endoplasmic reticulum glycoprotein processing enzyme GAA II. The identification of second-generation chaperones able to stabilize target enzymes without inhibitory action and undesired off-target effects is thus highly desirable. In this study, we were able to harness by structural and biochemical analyses the therapeutic value of the previously identified pharmacological chaperone NAC. On the basis of docking studies using a structural model of GAA derived from the X-ray structure of NtMGAM, the authors of the original report on the allosteric chaperoning function of NAC [START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF] proposed a NAC-binding site at the interface between the catalytic (β/α) 8 barrel and the distal β-sheet domain of GAA. Although located in proximity of the observed NAC1 binding site, this "virtual" pocket derived from a homology-based model is not present in the experimental 3D structure of rhGAA, emphasizing that homology models can provide guidance to virtual screening up to a certain extent, but ultimately experimental data are needed for precise mapping of binding sites. Interestingly, an astonishingly high number of small solvent molecules originating from the crystallization medium (glycerol, ethylene glycol, sulphate, chloride) are present in the different rhGAA structures presented here. Based on the principle of the multiple solvent crystal structures approach [START_REF] Mattos | Locating and characterizing binding sites on proteins[END_REF] , the different ligand binding sites might provide guidance to pharmacophore mapping in high-throughput screening experiments devoted to the identification of novel allosteric chaperones. We envisage that the structural framework provided here will allow mapping of disease mutations and guide clinicians in therapeutic choices and furthermore supply an accurate molecular scaffold for in silico screening to map novel druggable sites and identify novel pharmacological chaperones with complementary and synergistic effect on adverse GAA mutations.

Methods

Sample preparation. rhGAA (algucosidase alfa, Myozyme ® ) produced in CHO cells in the precursor form was from Genzyme, Cambridge, MA. As source of enzyme, the residual amounts of the infusions of Myozyme ® administrated for the treatment of Pompe patients at the Department of Translational Medical Sciences, Federico II University, Naples, Italy, were used. For the infusions, Myozyme ® was conditioned in a buffer composed of 110 mM mannitol, 36 μM polysorbate 80 (Tween80), 3 mM Na 2 HPO 4 •7H 2 O and 21 mM NaH 2 PO 4 •H 2 O, pH 6.1. One aliquot of the infusion residuals was dialyzed against a buffer composed of 10 mM HEPES pH 7.0 and 150 mM NaCl, whereas a second aliquot was dialyzed against 10 mM Na-citrate buffer pH 5.5, supplemented with 150 mM NaCl. Extensive crystallization trials with rhGAA in the original buffer and in the HEPES and citrate buffer conditions provided crystals diffracting at the best to about 7 Å resolution. Protein disorder prediction with the computational tool DisEMBL [START_REF] Linding | Protein disorder prediction: implications for structural proteomics[END_REF] (http://dis.embl.de/) indicated the existence of disordered regions within the rhGAA sequence and limited proteolysis was thus envisaged to remove potential disordered surface loops deleterious for the formation of well-ordered crystals. To this aim, rhGAA was dialyzed against a buffer composed of 10 mM HEPES pH 7.6, 200 mM NaCl and 2 mM CaCl 2 and α-chymotrypsin and subtilysin were tested at different protease to protein mass ratios ranging from 1:2 to 1:100, at different incubation times and incubation temperatures. At selected time points, aliquots were removed from the reaction mixture and analysed by SDS-PAGE. Best and reproducible results were obtained by α-chymotrypsin treatment overnight at 20 °C and with a protease/rhGAA mass ratio of 1:2, affording a band on SDS-PAGE (8%) with a mass difference of ~5 kDa with respect to the rhGAA precursor (Fig. 2a). The proteolysed enzyme was separated from the proteolytic reaction mixture by size-exclusion chromatography in 10 mM HEPES pH 7.0 and 100 mM NaCl using a HiLoad 16/60 Superdex 200 column, and concentrated to approximately 3.7 mg ml -1 using a centrifugal filter unit with a 30-kDa cutoff. N-terminal sequencing of proteolysed rhGAA was performed using a PROCISE ® 494 sequencing system.

Activity tests on the rhGAA forms. The activity assays were performed at 37 °C by adding 5 µg of precursor or mature forms of rhGAA to 0.2 ml of reaction mixture containing 20 mM 4NP-Glc and 100 mM Na-acetate pH 4.0. After 2 min the reactions were blocked by adding 0.8 ml of 1 M Na-carbonate pH 10.2. Absorbance was measured at 420 nm at room temperature and the extinction coefficient to calculate enzymatic units was 17.2 mM -1 cm -1 . One enzymatic unit is defined as the amount of enzyme catalyzing the conversion of 1 μmol substrate into product in 1 min, under the indicated conditions. Spontaneous hydrolysis of the substrates was subtracted by using appropriate blank mixtures without enzyme. All the kinetic data were calculated as the average of at least two experiments and were plotted and refined with the program Prism 5.0 (GraphPad Software, San Diego, CA USA).

Crystallization and structure solution. Crystallization trials were performed by the sitting-drop vapour-diffusion method using a Mosquito ® (TTP labtech) crystallization robot and commercial screens. First crystal hits emerged from a condition of the Structure Screen 1 (Molecular Dimensions) containing 2.0 M ammonium sulphate, 0.1 M HEPES pH 7.5 and 2% v/v PEG400, and after several optimization rounds well diffracting crystals were obtained by the hanging-drop vapour-diffusion method by mixing 1 µl protein solution at 3.7 mg ml -1 with 0.5 µl reservoir solution composed of 1.9 M ammonium sulfate, 0.1 M HEPES pH 7.0 and 2% v/v PEG400. Crystals belong to space group P2 1 2 1 2 1 with average unit cell axes of 97 × 103 × 128 Å and one molecule per asymmetric unit. Ligand soaking was performed by the addition of small amounts of powder of NAC, DNJ, NHE-DNJ or acarbose to the crystallization droplets (1-2 snowflake-like crystals of approximate dimensions of 0.5 mm), followed by incubation for ~1 h. Native rhGAA crystals were flash-cooled directly in liquid nitrogen, whereas crystals soaked with the various ligands were cryo-protected with reservoir solution supplemented with 15% (v/v) glycerol prior cooling. X-ray diffraction data were acquired at beam line Proxima1 at the Synchrotron Soleil, Gif-sur-Yvette, and at beam line ID23-2 at the European Synchrotron Radiation Facility, Grenoble. All diffraction data were processed with XDS 37 and scaled and merged using the CCP4 [START_REF] Winn | Overview of the CCP4 suite and current developments[END_REF] suite of programs Pointless, Aimless [START_REF] Evans | How good are my data and what is the resolution?[END_REF] and Truncate.

The structure of rhGAA was solved by molecular replacement with the program Phaser 40 using the coordinates of Homo sapiens NtMGAM (PDB ID 2QLY, 44% sequence identity with GAA) as search model (rotation function Z-score: 7.0; translation function Z-score: 21.4; number of clashes from packing analysis: 0; log likelihood gain: 1383). Subsequent auto-building with the program Buccaneer 41 yielded a model that was 93% complete with R free of 31.2%. Maximum-likelihood refinement, including TLS, and model adjustment were carried out with the programs Refmac [START_REF] Murshudov | Refinement of macromolecular structures by the maximum-likelihood method[END_REF] and Coot [START_REF] Emsley | Features and development of Coot[END_REF] , respectively. Random sets of ~5% of reflections, depending on the resolution limit, were set aside for cross-validation purposes. The composition of cross-validation data sets was systematically taken over from the parent data set. Model quality was assessed with internal modules of Coot [START_REF] Emsley | Features and development of Coot[END_REF] and with the Molprobity server [START_REF] Chen | MolProbity: all-atom structure validation for macromolecular crystallography[END_REF] . Ligands were fitted into unbiased F o -F c difference electron density maps calculated after 10 cycles of rigid body refinement with Refmac [START_REF] Murshudov | Refinement of macromolecular structures by the maximum-likelihood method[END_REF] . The data collection and refinement statistics are listed in Table 1 with representative electron density in Fig. 3, 4 and Supplementary Figs. 2 and8. Figures representing structural renderings were generated with the PyMOL Molecular Graphics System (DeLano, W.L. The PyMOL Molecular Graphics on http://www. pymol.org/). Chemical structures were drawn with ChemDraw 9.0 Pro (Cambrigesoft Corp, Cambridge, MA). Sequence alignments were made using MUSCLE [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF] and graphical rendering of the alignments, taking into account structural information, were made with ESPript [START_REF] Robert | Deciphering key features in protein structures with the new ENDscript server[END_REF] . Differential scanning fluorimetry. Thermal stability scans of rhGAA were performed as described in Porto et al., 2012 [START_REF] Porto | Pharmacological enhancement of α-glucosidase by the allosteric chaperone N-acetylcysteine[END_REF] . Briefly, the enzyme was diluted 48-fold (0.1 mg ml -1 ) in 25 mM Na-phosphate buffer, pH 7.4, 150 mM NaCl and 1:1000 SYPRO Orange dye (Life Technologies). Thermal scans were performed in triplicate in absence or in presence of 10 mM NAC, with steps of 1 °C min -1 in the range 25-95 °C in a Real Time LightCycler (Bio-Rad). Fluorescence was normalized to the maximum value within each scan to obtain relative fluorescence. Melting temperatures were calculated according to Niesen et al., 2007 [START_REF] Niesen | The use of differential scanning fluorimetry to detect ligand interactions that promote protein stability[END_REF] . The standard deviations for each melting temperature were calculated from three replicates.

The dissociation constant (K D ) of NAC was measured by thermal stability scans of rhGAA according to Vivoli et al., 2014 [START_REF] Vivoli | Determination of protein-ligand interactions using differential scanning fluorimetry[END_REF] . DSF scans were performed as described above, in the range 0-28 mM NAC. The melting temperature values were plotted as function of ligand concentration. The experimental data were best fitted according to a simple cooperative model equation reported in Vivoli et al., 2014 [START_REF] Vivoli | Determination of protein-ligand interactions using differential scanning fluorimetry[END_REF] by using the software GraphPAD Prism (GraphPad Software, San Diego, CA, USA). To measure the thermal stability of GAA from rice type V (resuspended in 2.8 M ammonium sulphate) (Sigma-Aldrich), the enzyme was diluted 34.5-fold (0.1 mg ml -1 ) and the thermal scans were performed as described above for rhGAA.

Kinetic constants of rhGAA on different substrates. All substrates were purchased from Sigma-Aldrich. Kinetic constants on 4NP-Glc (1-35 mM) were measured as described above. Kinetic constants on maltose (1-150 mM) and isomaltose (15-200 mM) were measured in 100 mM Na-acetate buffer pH 4.0 at 37 °C, by using 5 and 15 μg of enzyme, in a final volume of 0.2 and 0.1 ml, respectively. After 2 min, the amount of glucose produced in the reactions was determined using a glucose oxidase-peroxidase system (D-Glucose kit, Megazyme). One unit of activity was defined as the amount of enzyme releasing 1 μmol of glucose per minute at the conditions described. Spontaneous hydrolysis of the substrates was subtracted by using appropriate blank mixtures without enzyme. All the kinetic data were calculated as the average of at least two experiments and were plotted and refined with the program Prism 5.0 (GraphPad Software, San Diego, CA, USA).

Kinetic constants on bovine and rabbit glycogens were measured by using the substrates ranging from 1 to 123 mg ml -1 in 50 mM citrate-phosphate buffer pH 4.0 at 37 °C by using 5 μg of enzyme in the final volume of 0.1 ml. After 2 min, 0.1 ml of Na-carbonate 0.25 M was added to the reaction mixtures. To determine the amount of reducing sugars released 0.1 ml were analysed by the Somogyi-Nelson assay [START_REF] Smogyii | Notes on sugar determination[END_REF] . One unit of activity was defined as the amount of enzyme releasing 1 μmol of reducing equivalents per minute at the conditions described.

Fig. 1

 1 Fig.1The effect of pharmacological chaperones on misfolded lysosomal enzymes and on recombinant enzymes used in ERT. Lysosomal enzymes are assisted by molecular chaperones during synthesis. Mutated enzymes fail to fold correctly and are intercepted by the quality control (QC) system of the endoplasmic reticulum (ER). a Pharmacological chaperones favor proper folding of mutated enzymes, prevent their recognition by the quality control system and stabilize the enzyme during transport to their destination. b Pharmacological chaperones can enhance the effect of recombinant enzymes administered in ERT by favoring trafficking to lysosomes and increasing enzyme stability

Fig. 2

 2 Fig. 2 Structure of mature rhGAA. a Proteolytic treatment of rhGAA. b Schematic representation of the sequence of GAA. Myozyme ® rhGAA used in ERT starts at residue Q57. Domains corresponding to the rhGAA structure are colored as in c, with the regions removed by treatment with α-chymotrypsin colored in white. c Cartoon representation of the structure of rhGAA consisting of the trefoil type-P domain (salmon), the N-terminal β-sheet domain (slate), the catalytic GH31 (β/α) 8 barrel domain (green) with insert I (gold) and insert II (pink), and the proximal (orange) and distal (teal) β-sheet domains. Catalytic residues (magenta) and glycan chains (grey) are depicted as sticks. Dashed lined circles highlight peptide portions removed by proteolysis with α-chymotrypsin before crystallization. (SP, signal peptide; PP, propeptide)

Fig. 3

 3 Fig. 3 Ligand binding to rhGAA. DNJ (a), NHE-DNJ (c) and acarbose (e) colored in orange bound to rhGAA, color-coded as in Fig. 2c, overlapped onto unbound rhGAA (grey) in (a) and (c). Hydrogen bonding interactions are represented as dashed lines. In e, substrate-binding subsites are numbered and residues of a symmetry-related molecule are depicted in white sticks. Unbiased F o -F c difference electron density maps, calculated before incorporation of the ligands into the models and contored at 3.0 σ, are shown in b, d and f

Fig. 4

 4 Fig. 4 rhGAA substrate recognition and specificity and allosteric chaperone binding sites. a Model of isomaltose (steelblue), derived from an overlap with B. obeum α-glucosidase in complex with isomaltose (PDB ID 3MKK), superposed onto the rhGAA-acarbose complex (rmsd of 1.50 Å for 318 aligned Cα positions). b Model of isomaltose bound to rhGAA in surface representation. c The secondary substrate-binding site of rhGAA with unbiased F o -F c difference electron density map calculated before incorporation of acarvosine (orange) into the model and contoured at 2.0 (lightblue) and 3.0 (blue) σ. The surface loop removed by proteolytic treatment is represented by dashed lines. d NAC1 (orange) in the fully occupied binding site. e NAC2 (orange) in the partially occupied binding site. In d and e, unbiased F o -F c difference electron density maps, calculated before incorporation of NAC into the model, are shown in lightblue (2.0 σ) and blue (3.0 σ)

  NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-01263-3 ARTICLE NATURE COMMUNICATIONS | 8: 1111

Table 1

 1 Data collection and refinement statistics (molecular replacement) a, b

		rhGAA	rhGAA-NAC	rhGAA-DNJ	rhGAA-NHE-DNJ	rhGAA-acarbose
	Data collection
	Space group	P2 1 2 1 2 1	P2 1 2 1 2 1	P2 1 2 1 2 1	P2 1 2 1 2 1	P2 1 2 1 2 1
	Cell dimensions
	a, b, c (Å)	98,102,129	97,102,129	97,103,129	97,102,128	97,103,128
	Resolution (Å)	48.83-1.90	47.62-1.83	53.80-2.00	47.47-2.00	51.49-2.45
		(1.93-1.90)	(1.86-1.83)	(2.04-2.00)	(2.04-2.00)	(2.54-2.45)
	R merge	0.103 (1.711)	0.072 (0.569)	0.159 (0.549)	0.134 (0.602)	0.188 (0.885)
	R pim	0.050 (0.823)	0.040 (0.311)	0.125 (0.426)	0.091 (0.421)	0.116 (0.554)
	CC (1/2)	0.992 (0.661)	0.999 (0.859)	0.965 (0.675)	0.992 (0.722)	0.986 (0.644)
	I/σI	13.7 (1.2)	17.8 (3.3)	6.2 (2.0)	8.4 (2.3)	8.0 (1.9)
	Completeness (%)	100 (100)	99.6 (92.0)	95.3 (95.5)	99.5 (100)	94.7 (97.0)
	Redundancy	9.0 (9.2)	7.5 (7.1)	3.1 (3.0)	4.5 (4.2)	5.3 (5.2)
	Wilson B (Å 2 )	28.51	21.25	14.22	16.73	25.97
	Refinement
	Resolution (Å)	48.83-1.90	47.62-1.83	53.80-2.00	47.47-2.00	51.49-2.45
	No. reflections	97,044	107,574	79,052	81,565	42,513
	R work	16.19 (33.0)	15.37 (23.60)	18.34 (22.70)	15.26 (21.30)	17.42 (26.0)
	R free	18.87 (33.80)	17.83 (26.00)	21.81 (27.40)	18.38 (24.60)	21.23 (30.20)
	No. atoms
	Protein	6652	6669	6657	6653	6650
	Ligand/ion	272	329	256	303	339
	Water	668	776	773	633	290
	B-factors (Å 2 )
	Protein	31.46	24.01	17.81	20.74	28.87
	Ligand/ion	59.56	56.20	49.34	43.67	54.76
	Water	38.83	34.38	27.34	30.58	28.29
	R.m.s. deviations
	Bond lengths (Å)	0.009	0.008	0.009	0.007	0.007
	Bond angles (°)	1.407	1.327	1.366	1.297	1.346
	Ramachandran
	Favored (%)	97.96	97.74	97.72	97.24	97.74
	Allowed (%)	2.04	2.14	2.28	2.64	2.14
	Disallowed (%)	0	0.12	0	0.12	0.12

a One single crystal was used for each data collection b Values in parentheses are for the highest-resolution shell ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-01263-3 4 NATURE COMMUNICATIONS | 8: 1111

| DOI: 10.1038/s41467-017-01263-3 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | 8: 1111| DOI: 10.1038/s41467-017-01263-3 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | 8: 1111| DOI: 10.1038/s41467-017-01263-3 | www.nature.com/naturecommunications

Acknowledgements

We are grateful to S. Cottaz for providing glucopyranosyl-α-(1,6)-thio-maltotriose. We thank the European Synchrotron Radiation Facility (ESRF) and Synchrotron Soleil for beam time allocation and the beam line staff for assistance with data collection. This work was supported in part by the CNRS and the French Infrastructure for Integrated Structural Biology (FRISBI) ANR-10-INSB-05-01. Telethon support to G.P. is gratefully acknowledged.

Data availability. The atomic coordinates and structure factors of rhGAA and complexes with NAC, DNJ, NHE-DNJ and acarbose have been deposited in the Protein Data Bank with accession numbers 5NN3, 5NN4, 5NN5, 5NN6 and 5NN8, respectively. All the other data are available from the corresponding author on request.

Author contributions

Additional information

Supplementary Information accompanies this paper at doi:10.1038/s41467-017-01263-3.

Competing interests:

The authors declare no competing financial interests.

Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/ Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


published or not. The documents may come

Demographic inference using skyline plots on approximate bayesian computation

Miguel Navascués, Concetta Burgarella

We have implemented the same approach under an Approximate Bayesian Computation analysis. The purpose of this is to facilitate analysis of dataset were the MCMC approach is inappropriate, such as analysis of DNA sequences with intragenic recombination. The performance of the method has been evaluated with simulated datasets (see table ). Simulations for the ABC procedure were performed with FASTSIMCOAL (Excoffier & Foll 2011 Bioinformatics), using a model with 10 periods of constant size. Sample were of 50 diploid individuals typed at 100 microsatellite loci. Summary statistics were calculated with ARLSUMSTATS (from Arlequin). Relative root of mean squared error, relative bias and proportion of times the true value fall outside the 95%HPD interval were estimated from 50 simulations of target data.

  Yew (Taxus baccata) Burgarella et al. (2012) Mol Ecol 21:3006-21 Distribution in the Iberian Peninsula ABC Skyline Plots from seven microsatellite markers (model with 3 periods) probability density Bayesian Skyline Plots are representations of the posterior probability density of the effective population size in function of time. A model with several periods of constant population size is used to explore the demographies that best explain the data. Instead of estimating the posterior probability density for parameters of the model (i.e. size at each interval) the probability density of population size in function of time (i.e. size at each generation) is estimated. This method is currently implemented in BEAST (Drummond & Rambaut 2007 BMC Evol Biol), which uses an MCMC estimation of likelihood under the coalescent.

  UMR CBGP, INRA Montpellier, navascues@supagro.inra.fr UMR AGAP, INRA Montpellier, burgarella@supagro.inra.fr

		PRIOR	POSTERIOR (1 random simulation)	POSTERIOR (50 simulations)	ERROR
	CONSTANT SIZE
		EXPANSION
		BOTTLENECK
	EXPANSION	+ BOTTLENECK
				MEDIAN
				95%HPD


own prior to analysis. Further, in vitro for the control and knowledge of s, particularly the severity of population ng insight into the expected role of genetic utation rate and short generation time of on to occur on experimentally tractable t-up allows for an additional benefithole-genome data. This added temporal important component in the puzzle of nd demography -as it becomes possible ts describing the change in frequency [18] beneficial mutations. Thus, time-sampled y of any individual allele to be used to of natural selection, rather than simply variation as utilized by standard single spectrum based statistics [20]. ental approach and the above reasoning, loped novel statistical tests of selection for data. We infer effective population size nd develop novel analytical-, maximum ximate Bayesian -based approaches to ions of genetic drift and selection in this lly, based on this population genetic trate that IAV development of drug pectations of Fisher's Geometric Model, h to predicting viral evolution in response ments and novel selective pressures. S1). The frequency spectra indicated that most single nucleotide polymorphisms (SNPs) were segregating at low frequency and the

Figure 1 .TemporalF

 1 Figure 1. Experimental set-up. IAV was adapted from chicken egg to MDCK cells (passages 1-3) and then serially passaged in MDCK cells in the presence (red) or absence (white) of escalating concentrations of oseltamivir (passages 4-12) in replicate experiments. doi:10.1371/journal.pgen.1004185.g001

Acknowledgements

Comité d'organisation ETEE 2015 -Empirisme et Théorie en Écologie et Évolution Collaborators: Arnaud Becheler (currently at UMR EGCE, Gif-sur-Yvette) & Renaud Vitalis: Evaluation of temporal F ST outlier method Laurène Gay & Joëlle Ronfort (UMR AGAP): Medicago truncatula Project SelfAdapt funded by Métaprogramme ACCAF (INRA) Miguel Navascués (INRA, UMR CBGP) Detecting selection with temporal samples 30/06/2015

Selection on standing variation W = 0.61,FPR = 2.92 × 10 -4

Abstract

In a single isolated populations, allele frequencies will change through time subject to the processes of selection (acting on specific loci) and genetic drift (acting on the whole genome). Genetic data collected at different times can be used to make inferences on the effective population size (i.e. strength of drift) and to detect outlier loci, whose changes in allele frequencies are unlikely to be only the product of the inferred demography. However, the presence of self-fertilization may pose a problem for the detection of loci under selection. Selfing reduces the effective size of populations and the effective recombination among loci (promoting hitch-hiking). We investigated the effect of the presence of partial selfing reproduction in the power and false discovery rate for the detection of selected loci. In addition, we characterized the footprint of selection along the chromosome containing the selected site.