Mehdi Eshagh
email: eshagh@kth.se

Least-squares modification of extended Stokes' formula and its second-order radial derivative for validation of satellite gravity gradiometry data

Keywords: Biased, unbiased, optimum estimator, global root mean square error, truncation error, error spectra, signal spectra

The gravity anomalies at sea level can be used to validate the satellite gravity gradiometry data. Validation of such a data is important prior to downward continuation because of amplification of the data errors through this process. In this paper the second-order radial derivative of the extended Stokes' formula is employed and the emphasis is on least-squares modification of this formula to generate the second-order radial gradient at satellite level. Two methods in this respect are proposed: a) modifying the second-order radial derivative of extended Stokes' formula directly, b) modifying extended Stokes' formula prior to taking the second-order radial derivative. Numerical studies show that the former method works well but the latter is very sensitive to the proper choice of the cap size of integration and degree of modification.

Introduction

The satellite gravity gradiometry (SGG) is a technique to measure second-order derivatives of the Earth's gravity field from space. It is expected to determine the geopotential coefficients to higher degrees and orders than those are obtained from other satellite techniques. The SGG data can be used to study the geophysical/geodynamical phenomena as well. Quality of the data is important, as occurrence of any error in the data will lead to a wrong interpretation and unrealistic conclusions for the phenomena. Therefore, the quality of SGG data should be controlled prior to use, or in other words, the data should be validated.

Different methods of validating SGG data have been proposed. A simple way could be the direct comparison of the real SGG data with the synthesized gravitational gradients using an existing Earth's gravitational model (EGM). Another idea is to use regional gravity data to generate the gradients at satellite level. [START_REF] Haagmans | An alternative concept for validation of GOCE gradiometry results based on regional gravity[END_REF] and [START_REF] Kern | Determination of gravity gradients from terrestrial gravity data for calibration and validation of gradiometric GOCE data[END_REF] used the extended Stokes formula (ESF) and extended Hotine formula to generate the gravitational gradients using terrestrial gravity data. [START_REF] Denker | Computation of gravity gradients for Europe for calibration/validation of GOCE data[END_REF] used the least-squares spectral combination technique to generate and validate the gravitational gradients. Bouman et al. (2003) has set up a calibration model based on instrument (gradiometer) characteristics to validate the measurements. [START_REF] Mueller | Computation of calibration gradients and methods for in-orbit validation of gradiometric GOCE data[END_REF] used the terrestrial gravity anomalies to generate the gravitational gradients, and after that [START_REF] Wolf | Kombination globaler potentialmodelle mit terresrischen schweredaten fur die berechnung der zweiten ableitungen des[END_REF] investigated the deterministic approaches to modify the integrals and validate the SGG data. In fact, the spectral weighting scheme [START_REF] Sjöberg | Least-squares combination of satellite harmonics and integral formulas in physical geodesy[END_REF]and 1981[START_REF] Wenzel | Zur Geoidbestimmung durch kombination von schwereanomalien und einem kugelfuncationsmodell mit hilfe von integralformeln[END_REF] was used by Wolf A c c e p t e d M a n u s c r i p t (2007). Stochastic methods of modifying Stokes' formula, or in other words leastsquares modification (LSM) can be used for the extended Stokes formula as well; see Sjöberg (1984a[START_REF] Sjöberg | Least-squares modification of Stokes' and Vening Meinesz' formulas by accounting for errors of truncation, potential coefficients and gravity data[END_REF][START_REF] Sjöberg | Refined least-squares modification of Stokes' formula[END_REF][START_REF] Sjöberg | A general model for modifying Stokes' formula and its leastsquares solution[END_REF]. Least-squares collocation can be used for validation purposes. [START_REF] Tscherning | Calibration of GOCE gravity gradient data using smooth ground gravity[END_REF] considered this method and concluded that the gradients can be predicted with an error of 2-3 mE in the case of an optimal size of the collection area and optimal resolution of data. Zielinski and Petrovskaya (2003) proposed a balloon-borne gradiometer to fly at 20-40 km altitude simultaneously with satellite mission and proposed downward continuation of satellite data and comparing them with balloon-borne data. Bouman and Koop (2003) presented an along-track interpolation method to detect the outliers. Their idea is to compare the along-tack interpolated gradients with measured gradients. If the interpolation error is small enough the differences should be predicted reasonably by an error model. [START_REF] Pail | Local gravity field continuation for the purpose of in-orbit calibration of GOCE SGG observations[END_REF] proposed a combined adjustment method supporting high quality gravity field information within the well-surveyed test area for continuation of local gravity field upward and validating the SGG data. [START_REF] Bouman | Calibration of GOCE SGG data using high-low STT, terrestrial gravity data and global gravity field models[END_REF] stated that there are some limitations in generating the gravitational gradients using terrestrial gravimetry data and EGMs. When an EGM model is used, high degrees and orders should be taken into account and the recent EGMs seem to be able to remove the greater part of the systematic errors. In their regional approach they concluded that the bias of the gradients can accurately be recovered using leastsquares collocation. Also, they concluded that the method of validation using highlow satellite-to-satellite tracking data fails unless a higher resolution EGM is available. [START_REF] Kern | Determination of gravity gradients from terrestrial gravity data for calibration and validation of gradiometric GOCE data[END_REF] and [START_REF] Kern | Outlier detection algorithms and their performance in GOCE gravity field processing[END_REF] presented an algorithm for detecting the outliers in the SGG data in the time domain.

The second-order radial derivative (SORD) of extended Stokes' kernel (ESK) is isotropic and azimuth-independent. The isotropy is an important property in modifying ESF otherwise it will not be an easy task. Two methods of generating the SGG data are investigated in this paper, in the first method (Method 1), the SORD of ESF is modified (derivative prior to modification) and in the second method (Method 2) ESF is modified and after that the SORD is taken (modification prior to derivative). Modification of ESF and its SORD based on the biased LSM (BLSM), unbiased LSM (ULSM) and optimum LSM (OLSM) are the main subject of this study which is a new issue in the scope of SGG. Obviously, Methods 1 and 2 will not deliver the same results, but we are going to test in which cases these methods are comparable. We select the SORD of ESF as its kernel function is isotropic, in such a case, we can use both methods to generate the second-order radial gradient and compare the results. The importance of this study is mostly related to Method 2 although Method 1 (based on the LSM) is new as well. If we can find the cases, where Method 2 performs well, the method can be used to some how modify the horizontal derivatives of ESF having non-isotropic kernels, to generate the other gradients. A similar study was done by [START_REF] Wolf | Kombination globaler potentialmodelle mit terresrischen schweredaten fur die berechnung der zweiten ableitungen des[END_REF] but just based on deterministic approaches. However we concentrate on generating the second-order radial gradient based on the LSM approaches.

The disturbing potential can be expressed by an integral which is well-known as ESF. This integral formula is [START_REF] Heiskanen | Physical Geodesy. W.H Freeman and company[END_REF]:

A c c e p t e d M a n u s c r i p t ( ) ( ) ( ) , 4 R T P S r g Q d σ ψ σ π = ∆ ( ) ( ) ( ) 2 2 1 , c o s 2 n n n n S r r P ψ ψ ∞ = + = Ω ∑ , ( 1c
)
is the spectral form of ESK with the spectrum:

( )

1 2 1 n n R r n r +   Ω =   -  . ( 1d
)
Equation ( 1a) shows that the integration should be performed globally, which means that ( )

g Q ∆
with a global coverage is required. Therefore we should look for an approach to modify the integral in such a way that the contribution of the far zone data is minimized. Different methods for modifying Stokes' formula have been presented, but here the concentration is on the stochastic approaches of Sjöberg (1984a) and[START_REF] Sjöberg | Least-squares modification of Stokes' and Vening Meinesz' formulas by accounting for errors of truncation, potential coefficients and gravity data[END_REF]. In fact, the theory behind this part of the study was presented by him, but just on Stokes' integral for geoid determination. However, we are going to test the capability of these stochastic approaches in modifying ESF and its SORD and generating the SGG data for validation purposes. In the following we investigate the LSM of ESF.

LSM of ESF

The general estimator of the disturbing potential based on ESF is very similar to the general geoid estimator of [START_REF] Sjöberg | A general model for modifying Stokes' formula and its leastsquares solution[END_REF]; and the only difference is related to the kernel function and its spectrum. Let us start the discussion by this general disturbing potential estimator:

( ) ( ) ( ) ( ) ( ) 0 T E G M 2 , 4 2 L L n n n R R T P S r g Q d b r g P σ ψ σ π = = ∆ + ∆ ∑ ∫∫ , (2a)
A c c e p t e d M a n u s c r i p t

where L is the maximum degree of modification, ( ) n b r is a parameter which differs with the type of the LSM, and

( ) ( ) ( ) ( ) 2 2 1 , , c o s 2 L L n n n n S r S r s r P ψ ψ ψ = + = -∑ , ( 2b
)
is the modified ESF and ( ) n s r are the modification parameters, which are estimated. The closed form formula of this function is (Heiskanen and Moritz 1967, p. 93, Eq. 2-162):

( ) 2 2 2 2 c o s , 3 c o s 5 3 l n 2 R R Rl R r R l S r l r r r r ψ ψ ψ - +   = + - - +     , ( 2c
)
where

2 2 2 cos l r R Rr ψ = + - , (2d)
is the spatial distance between the points P and Q.

( )

EGM n g P ∆
is the Laplace harmonic expansion of ( ) Heiskanen and Moritz 1967, p. 97). In order to show from which sources the gravity anomaly is derived we separate them into The LSM parameters ( ) r s r are derived based on solving the following system of equations [START_REF] Sjöberg | A general model for modifying Stokes' formula and its leastsquares solution[END_REF]:

g P ∆ (
( ) ( ) 2 M kr r k r a s r h r = = ∑ , k = 2, 3,…, M, (2e)
where mathematical forms of kr a and

( ) k h r depend on type of the LSM which is used.

Equation (2e) differs with the system of equations in which the modification parameters of the Stokes formula is used for geoid determination. As Eq. (2e) shows both sets of the modification parameters and the truncation coefficients are altitudedependent and variable with the elevation of the computation point P. We will investigate the changes in these parameters and coefficients in Section 4. The mathematical formula of the elements of coefficient matrix and the right hand side vector of Eq. (2e) depend on the method of the LSM. In the following we summarized them in three propositions.

Proposition 1

The BLSM parameters for the disturbing potential estimator at satellite level are derived by setting )

2 0 0 2 nr nk n n n E E c ψ ψ σ ∞ = + + ∑ ( ) ( ) ( ) 2 0 , k k k k h r r Q r ψ σ   = Ω - +   ( )( ) ( ) ( ) 2 2 0 0 2 , n n n k n n k n Q r c r E ψ σ σ ψ ∞ =   + + -Ω   ∑ . k , r = 2, 3,…, M where ( ) ( ) 0 0 2 1 2 rk rk k E e ψ ψ + = , ( 3a
) ( ) ( ) ( ) 0 0 cos cos sin rk r k e P P d π ψ ψ ψ ψ ψ ψ = ∫ , ( 3b
) ( ) ( ) ( ) 0 0 , , c o s s i n k k Q r S r P d π ψ ψ ψ ψ ψ ψ = ∫ , ( 3c
) ( ) ( ) ( ) 0 0 2 2 1 , 2 k n n n n n Q r r e ψ ψ ∞ ′ = + = Ω ∑ . ( 3d
)
where 2 k σ is the error spectrum of the terrestrial gravimetric data and k dc is the error spectrum of the gravity anomaly obtained from an existing EGM.

For more details about the computation of 2 k σ the reader is referred e.g. to Sjöberg (1991), Ågren (2004) and [START_REF] Ellmann | Computation of three stochastic modification of Stokes' formula for regional geoid determination[END_REF]. The spectrum of the gravity anomaly is evaluated by using an EGM for those degrees below the maximum degree of modification, and analytical models like [START_REF] Kaula | The investigation of the gravitational fields of the moon and planet with artificial satellites[END_REF], [START_REF] Tscherning | Closed covariance expressions for gravity anomalies, geoid undulations and deflections of vertical implied by anomaly degree variance models[END_REF] are used for the degrees above that maximum degree. According to [START_REF] Ellmann | The geoid for the Baltic countries determined by the least squares modification of Stokes' formula[END_REF][START_REF] Ellmann | Computation of three stochastic modification of Stokes' formula for regional geoid determination[END_REF] and Ågren (2004) the Tscherning-Rapp model is superior with respect to the others. Relying on their conclusion we use this model in our numerical studies through this paper. (2e) with the following elements [START_REF] Sjöberg | A general model for modifying Stokes' formula and its leastsquares solution[END_REF]: (2e) with the following elements [START_REF] Sjöberg | A general model for modifying Stokes' formula and its leastsquares solution[END_REF]:

( ) ( ) ( ) ( ) 0 0 0 0 2 kr rk r kr rk r kr k nr nk n n a a d E d E d E E d δ ψ ψ ψ ψ ∞ = = = - - + ∑ , ( ) ( ) ( ) ( ) ( ) ( ) 2 2 0 0 0 2 , , k k k k k n n n n n k n h r r Q r d Q r d r E σ ψ ψ σ ψ ∞ =   = Ω - + -Ω   ∑ , k , r = 2, 3,…, M
( ) ( ) ( ) ( ) 0 0 0 0 2 kr rk r kr rk r kr k nr nk n n a a C E C E C E E C δ ψ ψ ψ ψ ∞ = = = - - + ∑ , k , r = 2, 3,…, M ( ) ( ) ( ) ( ) ( ) ( ) 2 2 0 0 0 2 , , k k k k k n n n n n k n h r r Q r C Q r C r E σ ψ ψ σ ψ ∞ =   = Ω - + -Ω   ∑ .
where

( ) 2 / 2 k k k k k k k c dc c dc k M C c k M σ  + ≤ ≤  = +  >   , ( 4a
) ( ) ( ) ( ) ( ) 0 0 0 2 , , L L n n n n k k Q r Q r s r E ψ ψ ψ = = -∑ . ( 4b
)

LSM of the SORD of ESF

The gravitational tensor is a symmetric tensor with 5 independent elements. This tensor is also traceless (trace of this tensor is zero) since the disturbing potential is harmonic outside the Earth's surface. Therefore, if we present a method for validating the simplest components, such as ( ) yy T P contains the blunder. In this section similar to the previous one we consider the LSM method of modifying the SORD of ESF to generate ( ) rr T P at point P at satellite level. We define the following general estimator ( ) where ( ) n s r ′′ are the modification parameters to be estimated, and (Reed 1973, Eq. 5.35)

rr T P : ( ) ( ) ( ) ( ) ( ) 0 T E G M 2 , 4 2 = ′′ = ∆ + ∆ ∑ ∫∫ L L rr rr n n n R R T P S r g Q d b r g P σ ψ σ π , ( 5a
( ) ( ) ( ) 2 3 2 2 5 3 3 3 1 4 1 10 , 1 c o s 1 8 rr t t t S r t D R D D D D ψ ψ    - +    = - - - -- +        1 cos 2 3 cos 15 6ln 2 t D t ψ ψ  - +   + - +       , ( 5c
)
where t R r = and 2 1 2 cos

D t t ψ = - + .
The LSM of SORD of ESF is very similar to that of ESF itself. If we estimate the modification parameter ( ) n s r ′′ , the left hand side of the system of equations (Eq. 2e) will not change and remain the same with those presented for the LSM of ESF. In such a case only the right hand side of Eq. (2e) will change. In the following, we present how the modification parameters ( )

( ) ( ) ( ) { } 2 0 , k k k k zz h r r Q r ψ σ ′′   = Ω - +   ( ) ( ) ( ) ( ) 2 2 0 0 2 , k n n k n n k zz n Q r c r E ψ σ σ ψ ∞ =  ′′    + + -Ω     ∑ , ( 6a
)
where

( ) ( ) ( ) 0 0 , , c o s s i n k r r k rr Q r S r P d π ψ ψ ψ ψ ψ ψ   =   ∫ , ( 6b
)
or in spectral form

( ) ( ) ( ) 0 0 2 2 1 , 2 k n n k rr n n Q r r e ψ ψ ∞ = + ′′   = Ω   ∑ , ( 6c
)
where , respectively, we can write the following estimator for the error of rr T : 7) [START_REF] Sjöberg | Refined least-squares modification of Stokes' formula[END_REF] proved that the integral part of this equation has the following spectral form

( ) ( )( ) ( ) 2 1 2 n n n n r r r + + ′′ Ω = Ω , ( 6d
( ) ( ) ( ) ( ) ( ) 0 T E G M 2 , 4 2 L L rr rr n n n R R T P S r Q d s r P σ δ ψε σ ε π = ′′ = + ∑ ∫∫ . (
( ) ( ) ( ) ( ) ( ) { } 0 T 0 2 , , 4
2

∞ * = ′′ ′′   = Ω - -  ∑ ∫∫ L T L rr n n n n rr n R R S r Q d r s r Q r σ ψ ε σ ψ ε π , (8a) where ( ) ( ) ( ) ( ) 0 0 0 2 , , L L n n k n k rr rr k Q r Q r s r E ψ ψ ψ = ′′     = -     ∑ . ( 8b
)
Substituting Eq. (8a) into Eq. ( 7) we obtain

( ) ( ) ( ) ( ) { } ( ) ( ) ( ) T E G M 0 2 2 , 2 2 ∞ * = = ′′ ′′ ′′   = Ω - - +   ∑ ∑ L L rr n n n n n n rr n n R R T P r s r Q r P s r P δ ψ ε ε . ( 8c
) By considering that the errors are random with zero stochastic expectation:

( ) ( ) T E G M 0 n n E E ε ε = = , ( 8d
)
we have

{ } 2 2 1 4 rr rr E T E T d σ δ δ σ π     = =       ∫∫ ( ) ( ) ( ) { } ( ) ( ) { } 2 2 2 2 2 0 0 2 , , 4 = ′′ ′′ ′′     = Ω - - + + +     ∑ L L L n n n n n n n n rr rr n R r s r Q r Q r c s r dc ψ σ ψ ( ) ( ) { } ( ) { } 2 2 2 2 0 0 1 , , 4 ∞ = + ′ ′     + Ω - +     ∑ L L n n n n n rr rr n M R r Q r Q r c ψ σ ψ . ( 8e
( ) ( ) ( ) ( ) ( ) { } ( ) 2 2 0 0 0 2 , , k k k k k k n n n n k rr rr n h r r Q r d Q r d r E σ ψ ψ σ ψ ∞ = ′′ ′′     = Ω - + -Ω     ∑ .
Proof. Similar to the proof presented for Corollary 1 we can write the error estimator as:

( ) ( ) ( ) ( ) ( ) { } ( ) 0 T E G M 0 2 , , 4 2 M L L rr rr n n n rr n R R T P S r Q d Q r s r P σ δ ψε σ ψ ε π = ′′   = + + -   ∑ ∫∫ ( ) ( ) { } ( ) ( ) ( ) 0 0 1 1 , , 2 2 L L T L T n n n n n rr rr n M n L R R Q r s r g P Q r g P ψ ψ ∞ = + = + ′′     - + ∆ - ∆     ∑ ∑ , (9a) where ( ) 0 , L n rr Q r ψ   
 was already presented in Eq. ( 8b). The mathematical expectation of Eq. (9a) shows that the estimator is unbiased through M, which is in fact the maximum degree of the EGM for generating the long wavelength structure of the estimator. By squaring the above equation and considering that the EGM and terrestrial data are not correlated and after taking the statistical expectation and the global average, as we did in Eq. ( 8e), we have the following spectral form for the global mean squares error:

{ } ( ) ( ) ( ) { } { 2 2 2 2 0 2 , 4 L rr n n n n rr n R E T r s r Q r δ ψ σ ∞ * = ′′ ′′   = Ω - - +   ∑ ( ) ( ) { } } 2 0 , L n n n rr Q r s r dc ψ * * ′′   + +   , (9b)
where

( ) ( ) 2 0 n n s r n L s r n L * ′′  ≤ ≤  ′′ =  >   and 2 n n n dc n L dc c n L * ≤ ≤  =  >  . ( 9c
( ) ( ) ( ) ( ) ( ) { } ( ) 2 2 0 0 0 2 , , k k k k k k n n n n k rr rr n h r r Q r d Q r d r E σ ψ ψ σ ψ ∞ = ′′ ′′     = Ω - + -Ω     ∑ .
A c c e p t e d M a n u s c r i p t

Proof. If the general error estimator Eq. ( 7) is written in the following spectral form:

( ) ( ) ( ) ( ) ( ) 0 T E G M 2 , 4 2 L L rr rr n n n R R T P S r Q d b r P σ δ ψε σ ε π = ′′ = + ∑ ∫∫ , ( 10a
)
then the spectral form of the above error estimator will be:

( ) ( ) ( ) ( ) { } ( ) ( ) ( ) T E G M 0 2 2 , 2 2 M L rr n n n n n n rr n n R R T P r s r Q r P b r P δ ψ ε ε ∞ * = = ′′ ′′ ′′   = Ω - - +   ∑ ∑ . ( 10b
)
The global mean square error of the estimator is:

{ } ( ) ( ) ( ) { 2 2 0 2 , 4 L rr n n n n rr n R E T b r s r Q r c δ ψ ∞ * * =   ′′ ′′   = - - +     ∑ ( ) ( ) ( ) } ( ) 2 2 2 2 0 2 , 4 M L n n n n n n rr n R r s r Q r b r dc ψ σ * =   ′′ ′′ ′′   + Ω - - +     ∑ , ( 10c
)
Taking the derivative of Eq. ( 10c) with respect to ( ) k s r ′′ and after further simplification the corollary is proved.

Numerical investigations

The EGM96 [START_REF] Lemoine | Geopotential model EGM96[END_REF]) is used to generate gravity anomaly and its error spectra. A correlation length of 0.1 and a standard error of 5 mGal are considered to generate the error spectra of terrestrial gravity anomaly. The EGM96 to degree and order M = 150 is used equal to the maximum degree of modification, L = 150. The BLSM, ULSM and OLSM methods are used to modify ESF and its SORD. The instability of Eq. (2e) to estimate the modification parameters in the ULSM and OLSM has been investigated by Ågren (2004), [START_REF] Ellmann | The geoid for the Baltic countries determined by the least squares modification of Stokes' formula[END_REF][START_REF] Ellmann | Computation of three stochastic modification of Stokes' formula for regional geoid determination[END_REF]. Here we take advantage of the truncated singular value decomposition method for regularizing Eq. (2e). According to [START_REF] Ellmann | The geoid for the Baltic countries determined by the least squares modification of Stokes' formula[END_REF] one can use any other regulator in this respect. Ågren ( 2004) also showed that the instability of this system of equations is harmless. The main goal of these studies was to modify Stokes' formula to determine a geoid. However, as was explained earlier, in order to construct a system of equations like Eq. ( 2e), to estimate the modification parameters, the derivative of global mean squares error of the estimator is taken with respect to ( ) k s r or ( ) k s r ′′ (which are altitude-dependent) directly, which lead to the same coefficients matrix as that was obtained by Ågren (2004) and [START_REF] Ellmann | Computation of three stochastic modification of Stokes' formula for regional geoid determination[END_REF]. Consequently, this matrix will have the same properties as that is used in modifying Stokes' formula and the conclusions made in these studies are also valid in our study. The important difference is due to the computation of the truncation coefficients and the spectra of ESF and/or its SORD.

A c c e p t e d M a n u s c r i p t hand side of Eq. (2e); and since the coefficient matrix of this system of equations is unstable then the harmlessness condition of the instability of the system is not preserved. In this case, the estimated modification parameters are defected by the integration error of the truncation coefficients. After some tests, we concluded to use the spectral form of the truncation coefficients (Eqs. 3d and/or 6c), involving Paul's coefficients ( )

0 nn e ψ ′
. The Paul's coefficients are not difficult to compute using the recursive formulas presented e.g. by [START_REF] Paul | Recurrence relations for integrals of associated Legendre functions[END_REF]. The only problem is that summation cannot be performed to infinity. We can select a large number (like 7000) as a truncation number. We tested and compared different truncation numbers for this series and we found that as long as the truncation number is above 5000 the truncation coefficients are more or less the same.

The truncation coefficients are altitude-dependent therefore the modification parameters will be altitude-dependent as well. If we look at Eqs. ( 2a) and (5a) we will find out that the second terms of both estimators which contribute the long wavelength structure of the estimators involve ( ) The modification parameters are the main factors for the changes in ESK. As was already explained these parameters are estimated in such way that ESK gets close to zero at the end of the cap size of integration, or in other words, it becomes sensitive to the existing data around the computation point. Effectiveness of modification can be presented by plotting the original and modified ESK in different geocentric angles in one plot. It should be reminded that Method 1 (taking derivative before modification) and Method 2 (modification before taking derivative) are considered in this subsection. Similar to the numerical studies performed in previous subsections the maximum degree of modification is L = 150 and the cap size of integration is 3 . Figures (3a), ( 3b) and (3c) show the modified ESK by Method 1 with the BLSM, ULSM and OLSM, respectively. In these figures different cases for the LSM are considered. First we assume that there is no error neither in terrestrial data ( 2 0

n σ = )
nor in the EGM based data ( 0 n dc = ). It means that we 100% rely on the quality of the data. Therefore the only error to minimize will be the truncation error of the SORD of ESF. In the second case, we consider these errors in the LSM. The third case considers no error in terrestrial data and the last case no error in the EGM data. The same situation holds for Method 2 in Figs. (3b), (3d) and (3f). Figures (3a), (3c) and (3c) express that the LSM is feasible and meaningfull in all cases. As Figs. (3b), (3d) and (3f) represent in some cases Method 2 works and in some does not. This method does not work when only the EGM error is considered. It is natural because we assumed that the terrestrial data are errorless then the estimator considers very high weight for them. By taking the SORD, the error of the EGM data is amplified; on the other hand, the estimator is subjected to errorless terrestrial data and concentrates on the data which contains a small portion of the gravitational signal (high frequency). Consequently, the ESK has to change largely to fulfill this constraint. 2e) with and without error of data. We present these figures to show that the global RMSE of the estimator is decreasing with increasing the cap size of modification (as it is expected). Figures (5c) and (5e) represent the global RMSE of the estimator Eq. (5a), which is the SORD of ESF, or in other words, Method 1. Similar pattern is observed for both cases in which the error of data is included or excluded. The errors are also decreasing by the increasing the cap size of modification. It is also reasonable as by adding more data in least-squares estimation the errors decay. Method 2 is the case that ESF is modified using the LSM methods and after taking the SORD of modified ESF the errors are estimated. Figure (5d) illustrates that when the SORD of the modification parameters in inserted to the mathematical formula of the global RMSE of the SORD of ESF the errors are still decreasing. Since we observed large variations in the errors we had to plot the vertical axis of the errors in logarithmic A c c e p t e d M a n u s c r i p t scale for better visualization. As can be seen in the figure, the errors are considerably large before the geocentric angle 2.5 for the modified estimators by the ULSM and OLSM with Method 2, while the error of the BLSM remains smaller than 0.01 E to the cap size 10 . The consequence of considering no error for the data will be similar. The main important conclusion that can be made here is related to the proper choice of cap size of modification in Method 2. One should be very careful to select the cap size prior to using this method; otherwise the method will yield wrong results. Figures (5d) and (5f) show that Method 2 works as long as the cap size is larger than 2.5 for the modified estimator by the ULSM and OLSM. Method 2 will work without problem for the BLSM case. A c c e p t e d M a n u s c r i p t

Generation of rr T over Fennoscandia by Methods 1 and 2

As was explained in the introduction of this section, the EGM96 is assumed as the true gravitational field. A regular grid of gravity anomalies were generated with the 30 30

′ ′ ×
resolution and the goal is to estimate rr T at 250 km level. Methods 1 and 2 are tested in different cap sizes of integration, resolutions of the terrestrial data, error of the terrestrial and the EGM data and the degree of modification. Fennoscandia limited between latitudes 50 and 75 and longitudes 0 and 35 is selected as the test region. rr T is generated in a grid with the same resolution of the gravity anomalies at sea level using Methods 1 and 2 and the LSM techniques. The statistics of the gravity anomalies and rr T are presented of the test area in Table 1. Gravity anomalies have more or less a smooth pattern in Fennoscandia. The maximum and minimum of the anomalies are situated in the Norwegian mountains and in the middle of Baltic Sea, respectively. The maximum and minimum of rr T are also in these two regions. Figure 6 shows the maps of the anomalies and rr T in Fennoscandia. Having a general view of the regime of the gravity anomaly and the gradient rr T , we continue the studies on modification and estimation of the gradients over this region.

Test of cap size of modification

The gradient rr T is generated directly from the EGM96 and is supposed to be the true gradient. This gradient is also generated using Methods 1 and 2. The numbers " 1 " and " 2 " in the tables of this section deal with these two methods. The generated gradients using these methods are compared with that were obtained directly from the EGM96 and the statistics of their differences (the errors of rr T ) are presented. Two cases are considered in all parts of the numerical studies, where the errorless and erroneous data are used in the LSM of ESF and its SORD. Table 2 shows the statistics of the errors of rr T for the caps sizes 3 and 2 . As the table shows when Method 1 is used and the errors of the data are set to zero, the BLSM yields the gradient with an error of 2 mE, while the ULSM and the OLSM do with a 13 mE error which is 6 times larger. If the cap size of modification reduces to 2 the BLSM's error reaches to 10 mE but error of the other LSMs is larger (35 times). Now let us consider Method 2. As we see in the table, this method can somehow estimate the gradient. If the BLSM is used, the error is 4 mE and it is roughly about 12 mE (three times larger) for the other LSM methods. By selecting a cap size of 2 , the ULSM and OLSM do not perform well and the errors are unrealistically large, while the BLSM still works with 6 mE error. If we return to Method 1 and include the data errors in the LSM, the error will be 4.6 mE for the BLSM and 7 mE for the ULSM and OLSM. The interesting issue is the improvement achieved by including the data error. When we consider zero error for the data, the estimator consider the same and very high weights for the data in the cap size and the EGM based data. However, when we add some errors to the terrestrial data, the errors are balanced in the LSM and the terrestrial data contributes according to their accuracies. One may ask why the BLSM based estimator was not improved. Actually existence of ( )

0 , L n Q r ψ in ( ) n b r
′′ increases the sensitivity of the estimator to the cap size in the ULSM and OLSM. In other words, some properties of the integral part of the estimator contribute to the second part, while this is not the case for the BLSM. In fact by reducing the weight of the data in the cap size of the integration the solution goes closer to that gradient obtained directly from the EGM and this is why we see better results in these two cases in the presence of the data errors. As we observe Method 2 does not work when the cap size is 2 and the sign " ! " we putted in the table means very big and unrealistic number which are due to weakness of the method in this cap size. In the case of using a cap size of 3 Method 2 works with an error of about 0.05 E. Now let us select a cap size of 2.5 . Table 2 shows the results of integration based on this cap size. Method 1 with errorless data yield 10 mE error based on BLSM and 20 mE based on the others. Method 2 yields 7 mE error with the BLSM while for the ULSM and OLSM is about 150 mE. In the case of including the data error again Method 1 with the BLSM seems to work with an error of 11 mE. Method 2 yields large errors when the ULSM and OLSM are used. We putted the sign " ! " instead of the very large numbers in the table. Methods 1 and 2, including the data error, yield the gradients with errors of 11 mE and 50 mE, respectively. One consequence of this section is that a cap size of 3 is suitable for integrating the terrestrial data. Selection of larger cap size will not be interesting. As Fig. 5 shows the SORD of ESK is very close to zero for larger geocentric angles and modification with such a cap size is not appropriate. [START_REF] Eshagh | On satellite gravity gradiometry[END_REF] studied the meaningfulness of the modification with a 3 cap size and found 0.1 E difference between the cases where the formula is modified and is not modified. We can learn from Tables 2 and3 that Method 2 will work if we rely 100% on the quality of the existing data. In this case the LSM methods perform like the deterministic methods in which just the truncation error of the formula is minimized. [START_REF] Eshagh | On satellite gravity gradiometry[END_REF] has also compared the Molodensky's modification [START_REF] Molodensky | Methods for study of the external gravity field and figure of the Earth[END_REF]) and compared it with the BLSM and found very similar results. Since errorless data is considered in the LSM of ESF then taking SORD, which acts as an amplifier operator, will not change the estimate significantly as there is no data error to be amplified. The BLSM seems to be the most suitable A c c e p t e d M a n u s c r i p t method for generating the gradient at satellite level. It seems smoother LSM than the others. Although the differences between these methods may not be significant theoretically but practically they are as we can see most of the friends of the LSM methods use the ULSM (Ågren 2004) or OLSM [START_REF] Ellmann | Computation of three stochastic modification of Stokes' formula for regional geoid determination[END_REF][START_REF] Kiamehr | Precise gravimetric geoid model for Iran based on GRACE and SRTM data and the least-squares modification of Stokes' formula with some geodynamic interpretations[END_REF][START_REF] Daras | Determination of a gravimetric geoid model of Greece using the method of KTH[END_REF][START_REF] Abdallah | Determination of a gravimetric geoid model of Sudan using the KTH method[END_REF] in modifying Stokes' formula for geoid determination. The reason could be the sensitivity of these LSM methods to the terrestrial data which are responsible to recover high frequencies of the geoid. As we showed in Table 3 Method 2 can generate the gradient with a 6 mE error which is smaller than that is estimated by Method 1 (10 mE). It should be mentioned that proper regularization of the system of equations Eq. ( 2e) in ULSM and OLSM is not meaningful as the instability of the system is harmless and the differences are not due to regularization. It is also similar to the results of [START_REF] Wolf | Kombination globaler potentialmodelle mit terresrischen schweredaten fur die berechnung der zweiten ableitungen des[END_REF] who concluded that the caps size 2.5 is the best. As it is clear from her thesis, she used the deterministic approaches and the spectral weighting scheme in her computations which is in agreement with Method 2 with the BLSM and with errorless data.

Test of resolution and error of terrestrial data

In order to investigate the dependence of Methods 1 and 2 to the terrestrial data we consider one test on the resolution and another on the error of terrestrial and the EGM data. Table 4 shows the statistics of the errors of the estimated rr T by Methods 1 and 2 with the resolution 15 15 ′ ′ × and with errorless and erroneous data. Insignificant changes are observed in the table, as there is no information higher than the nyquist frequency of the EGM96. The resolution of this EGM is 30 30 ′ ′ × corresponding to degree and order 360. Consequently, the enhancement of the results (if any) will be very small as the attenuation factor reduces higher frequencies at satellite level. Clearly, if an EGM with higher resolution is used the data with higher resolution may not be used. In the second part of our investigation, we work with the 30 30 ′ ′ × data and test the influence of the error of terrestrial data. At first, we consider no error in the EGM based data namely n dc =0 and modify the estimators assuming 5 mGal error in terrestrial data. That part of the table which is related to this study is separated by 2 n σ meaning that there is terrestrial data error in the LSM. We see very considerable reduction in the error of the estimated gradients. As we see Methods 1 and 2 yield the same error of 2 mE and the error of Method 1 based on ULSM and OLSM is about 3 mE. Such a good result could be expected as by including the error of the terrestrial data in the LSM process the terrestrial data are down-weighted and the methods yield closer results to that gradient directly obtained from the EGM. Reversely, if we consider no error in the terrestrial data and up-weight the data, the contribution of the EGM is reduced. Since the long wavelength structure of the gravitational field has the major contribution in satellite data, therefore putting high weights for a small portion of data, which are restricted to a certain cap size, cannot estimate the gradient well at satellite level. This is why we observe large errors for Method 1 with erroneous data. The situation will be even worse for Method 2 as the errors are considerably amplified by taking SORD. This subsection shows the importance of the quality of the EGM which is used.

Test of degree of modification

Degree of modification shows the contribution of the long wavelength structure of the EGM used in the LSM. Here we consider four degrees of modification 100, 125, 175 and 200. We already investigated the degree 150 in the previous subsections. Table 5 shows the statistics of the errors of rr T estimated by Methods 1 and 2 with degrees of modification of 100 and 125. When L = 100 we see large error for Method 1 and in the same order as that is for all the LSM methods considering errorless data. Method 2 delivers the error in the same level of Method 1 for BLSM, but the errors are twice larger for the others. When the errors of the data are included, Method 1 performs better with an error of 5 mE for all the LSM methods. However, the errors of Method 2 are considerably larger than that of Method 1. If the degree of modification increases to 125 the results seem better in Method 1 with errorless data. Method 2 is good just with the BLSM. Method 1 with the BLSM and in the presence of erroneous data performs similarly with an error of 7 mE. As was explained before, the ULSM and OLSM are more sensitive to the terrestrial data than the BLSM and when the degree of modification, or in other words, the EGM contributes further and Method 1 using the BLSM yields smaller error, while when higher degrees of the EGM are included the their errors are cumulatively increased as well. T . Method 1 yields rr T with a 2.7 mE error based on the BLSM and 9 mE based on the ULSM and OLSM when the data are errorless and L = 175. If such errorless data are considered for Method 2 the errors of the method will be 1.8 mE for the BLSM and 8 mE for the others. If the degree of modification increases to 200 we will not see significant improvement in the methods with the BLSM, but the errors will be reduced by about 2 mE for the ULSM A c c e p t e d M a n u s c r i p t and OLSM. Similar results can be achieved if the erroneous data are considered in the modification, but the error of the method will increase by about 2 mE. The other LSMs will improve the estimation about 1 mE. The error of Method 2, in the presence of the erroneous data, is about 42.0 mE for the L = 175 and L = 200, which means that after a certain degree Method 1 based on the BLSM will not improve the estimation. However, this is not the case for the ULSM and OLSM as very significant differences are observed in the errors. The error of Method 2 is about 14 mE when L =175 while it is 178 mE when L =200. This error is 8.6 mE in the former case while it is 212 mE in the latter. One can conclude that it is possible to find an optimum degree of modification for the ULSM and OLSM. Therefore, we can expect to find such a degree for Method 2 based on the ULSM and BLSM if the modification degree is properly selected. We have further investigated Method 2 at different degrees of modification. We found that the minimum achievable error using Method 2 and the ULSM is 11.8 ~ 12 mE at L = 170 and 6.6 ~ 7 mE around L = 189 for the OLSM. Although Method 1 yields good results in most cases but the main important assumption of this method is the isotropy of ESK. In this study we used the SORD of ESF because in the both methods we have isotropic kernels and possibility of comparing the methods. Method 2 works in special cases and one should be careful in the use of this method. Selection of the proper cap size and degree of modification is very important in this method. If these parameters are selected correctly Method 2 can be used for generating the gradient with error within 0.01 E. This method can also be used in generating the other elements of the tensor of gravitation.

Conclusions

In this paper the LSM methods were used for ESF and its SORD to generate rr T at 250 km level. The mathematical formulas of elements of the system of equations from which the modification parameters are obtained based on the LSM methods were developed. We have shown that if we organize the system of equations in such a way that the altitude-dependent modification parameters and/or its SORD are directly solved the properties of the system of equations is the same as that is for the LSM of ordinary Stokes' formula for geoid computation. Two methods were presented to generate rr T which were named Methods 1 and 2. Method 1 seems working very well depending on the accuracy of the data proper choice of the cap size of the integration and degree of modification. This method produced the best results when the degree of modification is L = 150 and the resolution of data is 30 30 ′ ′ × in a cap size 3 and if we consider an error-free EGM. Definitely by increasing the degree of modification better results would be achieved. Method 2 is very sensitive to proper choice of these parameters otherwise the results will be wrong. The reason is the amplification of the errors by taking SORD from the modified ESF. In this study we found that if the cap size of integration is around 3 and the degree of modification is about 175 we can generate rr T with 7 mE error using the ULSM and OLSM. If we assume no error in the LSM except the truncation error of ESF, Method 2 with the BLSM will work with sufficient accuracy. This method is very similar to deterministic modification and therefore we can conclude Method 2 with a deterministic modifier should work.

  solving the system of equations Eq. (2e) with the following elements[START_REF] Sjöberg | A general model for modifying Stokes' formula and its leastsquares solution[END_REF]

  error coefficients of ESK and the Paul coefficients (Paul 1978), respectively. Proposition 2 The ULSM parameters of the disturbing potential estimator at satellite level are derived, if we select (

  The OLSM parameters of the disturbing potential estimator at satellite level are derived, if we set (

Corollary 1

 1 At the first step, we start with the following corollary which is related with the BLSM of the SORD of ESF. The BLSM parameters are derived by setting ( equations Eq. (2e) with the following right hand side

  was introduced in Eq. (1d).

  The only difference between the biased geoid estimator of[START_REF] Sjöberg | Refined least-squares modification of Stokes' formula[END_REF] and the biased rr T estimator is related to the truncation coefficients and spectrum of the SORD of ESF. If we consider the errors of terrestrial and satellite data by T ε and EGM n ε

  result to zero and rearranging the known and the unknown parameters ( ) k s r ′′ , the corollary is proved.

  system of equations Eq. (2e) by considering

  ) By differentiating Eq. (9b) with respect to the modification parameters to zero and after further simplifications the corollary is proved.

  as well (because of altitude-dependency of the truncation coefficients and the modification parameters).

Fig

  Fig. 1. a) and b)

  s c r i p t presented them in Figs. (1b) and (1d). Figures (1c) and (1d) are ( ) n b r for the BLSM. Because of having very similar plots for the ULSM and the OLSM we did not present their ( ) n b r . However, the consequence of their presentation will be the same as that of the BLSM. Figure (1b) which is a zoomed part of Fig. (1a), shows that the magnitude of truncation coefficients and modification parameters decreases by increasing altitude.

Fig

  Fig. 2.

  Fig. 3. a), c) and e) modified ESK by Method 1 using BLSM, ULSM and OLSM, respectively, b), d) and f) corresponding modified ESK by Method 2, L = 150 and 0 3 ψ =

Figures

  Figures (5a) and (5b) show the global RMSE of Eq. (2e) with and without error of data. We present these figures to show that the global RMSE of the estimator is decreasing with increasing the cap size of modification (as it is expected). Figures (5c) and (5e) represent the global RMSE of the estimator Eq. (5a), which is the SORD of ESF, or in other words, Method 1. Similar pattern is observed for both cases in which the error of data is included or excluded. The errors are also decreasing by the increasing the cap size of modification. It is also reasonable as by adding more data in least-squares estimation the errors decay. Method 2 is the case that ESF is modified using the LSM methods and after taking the SORD of modified ESF the errors are estimated. Figure(5d) illustrates that when the SORD of the modification parameters in inserted to the mathematical formula of the global RMSE of the SORD of ESF the errors are still decreasing. Since we observed large variations in the errors we had to plot the vertical axis of the errors in logarithmic

Fig. 5 .

 5 Fig. 5. Global RMSE of a) disturbing potential with erroneous data, b) with errorless data, of c) rr T with erroneous data Method 1, d) rr T with erroneous data Method 2, e) rr T with errorless data Method 1, f) rr T with errorless data Method 2.

Fig. 6 .

 6 Fig. 6. a) Gravity anomalies (mGal) and b) rr T (E ) in Fennoscandia

Table 1 .

 1 Statistics of gravity anomalies and rr T in Fennoscandia

		Max	Mean	Min	Std
	rr T ( E)	0.430	-0.002	-0.455	0.187
	g ∆ ( mGal)	73.844	3.582	-54.643	17.284

Table 2 .

 2 Statistics of errors of rr T for the cap sizes 0 3 =

									ψ	and 2 with 30 30 ′ ×	′	resolution for the
				gravity anomalies and gradient and with L=150. Unit : 1 mE
					ψ	0	=	3			ψ	0	=	2
				Max	Mean			Min	Std	Max	Mean	Min	Std
	Errorless data	1 2	BLSM ULSM OLSM BLSM ULSM	4.2 14.6 14.7 17.6 40.2	-0.4 -3.9 -3.9 1.4 3.2			-11.6 -60.7 -61.5 -7.2 -19.6	2.2 12.7 13.0 4.0 11.1	13.7 689.4 710.0 19.1 3309.5	-3.1 -96.1 -100.3 1.5 319.2	-46.9 -999.5 -1041.0 -13.9 -2433.0	10.3 349.2 364.2 6.1 1200.5
			OLSM	41.6	3.3			-21.5	11.6	3504.6	337.9	-2545.9	1268.2
			BLSM	13.2	-0.2			-12.4	4.6	25.0	-6.0	-88.8	20.0
		1	ULSM	22.9	1.0			-17.8	7.5	2469.0	241.8	-1788.7	892.5
			OLSM	23.3	1.0			-18.2	7.6	2406.4	235.8	-1721.7	868.2
			BLSM	125.8	9.2			-85.4	41.6	163.0	12.8	-154.3	64.1
			ULSM	176.9	12.5			-112.7	55.1	!

Table 3 .

 3 Statistics of error of rr T for the cap size 0 2.5 =

				ψ				with 30 30 ′ ×	′	resolution for the gravity
		anomalies and gradient with L=150. Unit: 1 mE
					ψ	0	=	2.5
				Max	Mean	Min	Std
	Errorless data	1 2	BLSM ULSM OLSM BLSM ULSM	13.1 33.1 33.3 22.9 401.3	-3.1 -6.2 -6.3 1.9 34.7		-50.1 -83.3 -85.3 -13.6 -287.4	10.3 21.3 21.9 6.72 144.3
			OLSM	424.5	36.9		-303.4	153.2
	Erroneous data	1 2	BLSM ULSM OLSM BLSM ULSM	15.9 14132.2 117.7 1004.6 502.3 -3.6 -56.2 10.7 1388.4 117.0 -985.6 497.8 148.0 11.0 -147.8 52.4 ! ! ! !
			OLSM	!	!			!	!

Table 4 .

 4 Statistics of errors of rr T for testing resolution and error of terrestrial data in Methods 1 and 2

					for 0 3 = ψ	and L =150. Unit: 1 mE
					15 15 ′ ′ ×				30 30 ′ ′ ×
				Max	Mean	Min	Std	Max Mean	Min	Std
	Errorless data	1 2	BLSM ULSM OLSM BLSM ULSM OLSM	4.5 17.0 16.7 18.0 39.7 41.0	-0.5 -4.0 -4.1 1.4 3.0 3.0	-13.5 -64.8 -65.7 -7.6 -21.0 -22.9	2.5 13.6 13.9 4.2 11.3 11.8	6.0 6.8 9.5 6.0 71.8 74.3	0.5 0.17 0.3 0.5 5.7 6.1	-6.3 -8.0 -9.4 -6.3 -39.9 -43.1	2.1 2.8 3.3 2.1 20.8 22.3	σ	2 n
	Erroneous data	1 2	BLSM ULSM OLSM BLSM 122.2 13.8 20.8 21.3 ULSM 167.4 OLSM 168.9	-0.5 0.5 0.6 8.2 10.7 10.9	-14.9 -18.4 -18.8 -88.4 -108.6 52.0 4.6 7.1 7.3 40.6 335.8 18.3 50.4 -11.3 55.2 -12.2 55.2 -12.2 ! ! -109.0 52.6 ! !	-165.6 -175.9 -175.9 -263.2 110.1 37.6 40.6 40.6 ! ! ! !

n dc

Table 6

 6 is similar to Table5but for the modification degrees 175 and 200. Some improvements are seen in the estimated rr

Table 5 .

 5 Statistics of errors of rr T estimated based on Methods 1 and 2 for 0 3 =

	ψ	and resolution

Table 6 .

 6 Statistics of errors of rr T estimated based on Methods 1 and 2 for 0 3 =

	ψ	and resolution

Acknowledgements

The author would like to appreciate Professor Lars E. Sjöberg for the scientific discussion about least-squares modification. The Swedish National Space Board (SNSB) is cordially acknowledged for the financial support, project no. 63/07:1. Mr. Mohammad Bagherbandi is appreciated for his criticisms on the draft version of the paper. Professor Randell Stephenson and his reviewing board are acknowledged.

A c c e p t e d M a n u s c r i p t

Here, we consider two different approaches to generate rr T at satellite level, either modifying the disturbing potential estimator, Eq. (2a) and after that taking the SORD of the modified estimator or modifying the rr T estimator, Eq. (5a), directly. Obviously, these methods will not yield the same results, as the former method minimizes the global mean square error of the disturbing potential and not its SORD. However, we will investigate in which cases they yield similar results in practice.

This section consists of four parts. As was mentioned, the modification parameters and the truncation coefficients are altitude-dependent. Subsection 4.1 will consider this matter for ESF and its SORD. The modified ESK and its SORD are dependent on altitude, error spectral of the EGM or terrestrial data as well as the truncation error of the integral. Subsection 4.2 assumes a constant altitude and investigates different cases of modification based on the error spectra of the data in Methods 1 and 2. Also it will illustrate 4.1 ( )

altitudes

So far most of the geodesists worked on modification of Stokes' formula. In such a case the truncation coefficients of Stokes' kernel are not altitude-dependent and it is enough to generate the modification parameters once during the integration of Stokes' integral. In airborne gravimetry the truncation coefficients are altitudedependent but the aircraft's height is restricted to few kilometers, say e.g. to 5 km, which is considerably lower than satellites' altitude. In this way of gravimetry the integration is performed after downward continuation of the airborne data. Few persons considered modification at airborne level; see e.g. [START_REF] Sjöberg | Considering data gaps in geoid modeling by modifying Stokes's formula[END_REF]. If regeneration of the airborne data is required, one can estimate the truncation coefficients just once as the aircraft height is approximately constant. The situation at satellite level is complicated; firstly the altitude of a satellite is very large. Secondly due to very large orbital perturbations the satellite altitude cannot be regarded constant. [START_REF] Shepperd | A recursive algorithm for evaluating Molodeski-type truncation error coefficients at altitude[END_REF] used a similar mathematical technique as used by Hagiwara (1978) to derive a recursive formula for generating the altitudedependent truncation coefficients. He obtained the formula but it was impractical because of instabilities with respect to elevation (although his mathematical derivations were correct). [START_REF] Shepperd | A recursive algorithm for evaluating Molodeski-type truncation error coefficients at altitude[END_REF] stated that this formula can work for the altitudes not higher than 20 km. Therefore one concludes that [START_REF] Shepperd | A recursive algorithm for evaluating Molodeski-type truncation error coefficients at altitude[END_REF] formula works for airborne gravimetry aspects. The other idea is to carry out a numerical integration process to solve the integral formula of the truncation coefficients, namely Eqs. (3c) or (6b). We numerically tested this matter and concluded the non-efficiency of this idea. Our studies showed that the integration error defects the result and increasing the integration step size would be very time consuming. One may use this method for the cases where the altitude is constant and the numerical integration should be performed once. However, this method should not be used in the ULSM and OLSM as the truncation coefficients are in the right

Methods 1 and 2 in a global RMSE point of view

Equations ( 2a) and (5a), or ESF and its SORD are modified in a least-squares sense, which means that the modification parameters are estimated in such a way that the global RMSE of the estimators is minimized. Substituting the estimated modification parameters into the mathematical expression of the global RMSE yields the global RMSE of each estimator. This error can be a criterion to judge about quality of the LSM and Methods 1 and 2.

List of abbreviations

Biased least-squares modification……………………………………………....BLSM Earth's gravitational model………………………………………………………EGM Ågren J. (2004) Regional geoid determination methods for the era of satellite gravimetry, Numerical investigations using synthetic Earth gravity models, Doctoral thesis in Geodesy, Royal Institute of Technology, Stockholm, Sweden.


Francisco Muñoz-Martínez

Cristina Torres

Santiago Castanys
email: castanys@ipb.csic.es

Francisco Gamarro
email: gamarro@ipb.csic.es.

Biomedicina López-Neyra

CDC50A plays a key role in the uptake of the anticancer drug perifosine in human carcinoma cells

Functional aminophospholipid translocases are composed of at least two proteins: an alpha subunit from the P4 subfamily of P-type ATPases and a beta subunit from the CDC50-Lem3p family. Over-expression and knockdown of the human beta subunit CDC50A in KB cells enhanced and decreased, respectively, the uptake of both fluorescent aminophospholipid analogues and the anticancer alkylphospholipid perifosine. Confocal microscopy showed that CDC50A-V5 was localized at the endoplasmic reticulum and the Golgi complex of both KB (perifosinesensitive) and KB PER-R (perifosine-resistant, alkyl-phospholipid uptake deficient) cells, but was only widely distributed in the early and late endosomes in KB cells.

Biotinylation of cell surface proteins allowed CDC50A-V5 to be detected in the plasma membrane of KB cells but not in KB PER-R cells, thereby suggesting a defect in CDC50A trafficking that could explain the inability of KB PER-R to uptake perifosine. Over-expression of CDC50A in HeLa and HEK293T cells did not increase uptake, since the protein was retained at the endoplasmic reticulum and Golgi. However, when CDC50A was co-expressed with the P4-ATPase Atp8b1, the two proteins co-localized at the plasma membrane and the uptake of aminophospholipids and perifosine increased strikingly in both cell lines. These findings suggest that CDC50A plays a key role in perifosine uptake in human cells, presumably by forming a functional plasma membrane translocator in combination with a P4-ATPase. KEYWORDS: perifosine uptake, CDC50A, Atp8b1, P4-ATPase, aminophopholipid translocase, resistance to perifosine, cancer chemotherapy.

INTRODUCTION

Perifosine is an oral alkyl-phospholipid (ALP) Akt inhibitor that is currently being tested in phase II clinical trials. Unlike most kinase inhibitors, which target the adenosine triphosphate-binding region, perifosine targets the pleckstrin homology domain of Akt, thereby preventing its translocation to the plasma membrane [START_REF] Kondapaka | Perifosine, a novel alkylphospholipid, inhibits protein kinase B activation[END_REF].

Single-agent activity with perifosine has been observed in sarcoma and Waldenström macroglobulinemia patients [START_REF] Gills | Perifosine: update on a novel Akt inhibitor[END_REF]. However, the disappointing response rates of common solid tumors to perifosine as a single agent [START_REF] Chee | The AKT inhibitor perifosine in biochemically recurrent prostate cancer: a phase II California/Pittsburgh cancer consortium trial[END_REF][START_REF] Leighl | A Phase 2 study of perifosine in advanced or metastatic breast cancer[END_REF][START_REF] Argiris | A phase II trial of perifosine, an oral alkylphospholipid, in recurrent or metastatic head and neck cancer[END_REF][START_REF] Knowling | A phase II study of perifosine (D-21226) in patients with previously untreated metastatic or locally advanced soft tissue sarcoma: A National Cancer Institute of Canada Clinical Trials Group trial[END_REF] have diminished expectations and prompted further investigation into its mechanism of action.

The ability to identify patients likely to respond to targeted therapeutics would decrease the size, duration, and cost of clinical trials, thus resulting in more efficient translation to improved patient outcomes. In this sense, Hennessy and co-workers [START_REF] Hennessy | Pharmacodynamic markers of perifosine efficacy[END_REF] reported a strong correlation between proportional modulation of PI3K pathway biomarkers and quantified perifosine efficacy, thereby suggesting that the early measurements of such biomarkers will optimize the selection of responsive patients and guide perifosine dosage during cancer treatment. Furthermore, the identification of other potential biomarkers for perifosine efficacy in cancer patients would allow those patients most likely to respond to perifosine, either alone or in combination with other therapies, to be selected. ALP uptake is a prerequisite for triggering the intracellular events that finally cause cell death in different cellular systems, including yeasts [START_REF] Hanson | Lem3p is essential for the uptake and potency of alkylphosphocholine drugs, edelfosine and miltefosine[END_REF][START_REF] Riekhof | The yeast plasma membrane P4-ATPases are major transporters for lysophospholipids[END_REF], the protozoan parasite Leishmania donovani [START_REF] Pérez-Victoria | Phospholipid translocation and miltefosine potency require both L. donovani miltefosine transporter and the new protein LdRos3 in Leishmania parasites[END_REF] and cancer cells [START_REF] Gajate | Intracellular triggering of Fas, independently of FasL, as a new mechanism of antitumor ether lipid-induced apoptosis[END_REF]. Indeed, a clear correlation between ALP uptake and apoptosis induction has been reported previously [START_REF] Mollinedo | Selective induction of apoptosis in cancer cells by the ether lipid ET-18-OCH3 (Edelfosine): molecular structure requirements, cellular uptake, and protection by Bcl-2 and Bcl-X(L)[END_REF].

Two pathways for the entry of ALP into cancer cells, namely raft-mediated endocytosis in lymphoma cells [START_REF] Van Der Luit | A new class of anticancer alkylphospholipids uses lipid rafts as membrane gateways to induce apoptosis in lymphoma cells[END_REF] and translocator-mediated uptake in carcinoma Data from our previous work [START_REF] Muñoz-Martínez | The anti-tumor alkylphospholipid perifosine is internalized by an ATP-dependent translocase activity across the plasma membrane of human KB carcinoma cells[END_REF] strongly supported the hypothesis that a putative aminophospholipid translocase (APLT; P-type ATPase P4 subfamily) could be responsible for the uptake of ALPs in human KB carcinoma cells. Indeed, the plasma membrane P4-ATPase homologues of yeast and Leishmania have been shown to play a key role in the uptake of miltefosine and edelfosine [START_REF] Riekhof | The yeast plasma membrane P4-ATPases are major transporters for lysophospholipids[END_REF][START_REF] Pérez-Victoria | Functional cloning of the miltefosine transporter. A novel P-type phospholipid translocase from Leishmania involved in drug resistance[END_REF] in combination with a second protein from the Cdc50p-Lem3p family that functions as a beta subunit [START_REF] Hanson | Lem3p is essential for the uptake and potency of alkylphosphocholine drugs, edelfosine and miltefosine[END_REF][START_REF] Pérez-Victoria | Phospholipid translocation and miltefosine potency require both L. donovani miltefosine transporter and the new protein LdRos3 in Leishmania parasites[END_REF][START_REF] Sanchez-Cañete | Low plasma membrane expression of the miltefosine transport complex renders Leishmania braziliensis refractory to the drug[END_REF][START_REF] Saito | Cdc50p, a protein required for polarized growth, associates with the Drs2p Ptype ATPase implicated in phospholipid translocation in Saccharomyces cerevisiae[END_REF]. As fourteen putative APLTs [START_REF] Folmer | P4 ATPases -lipid flippases and their role in disease[END_REF] and only three beta subunits (CDC50A, B and C) [START_REF] Katoh | Identification and characterization of CDC50A, CDC50B and CDC50C genes in silico[END_REF] have been identified in humans, we decided to try to elucidate the role of these human beta subunits as functional elements of the putative plasma membrane perifosine translocator (and other related ALPs) in cancer cells, and found that CDC50A plays a key role in the uptake of perifosine in human cancer cells, presumably by forming a functional plasma membrane translocator in combination with a P4-ATPase.

MATERIALS AND METHODS

2

MTT survival assays, [ 14 C]-perifosine uptake, accumulation of fluorescent phospholipid analogues and measurement of endocytosis rates.

Unless otherwise indicated, these four different assays were performed as described in our previous work [START_REF] Muñoz-Martínez | The anti-tumor alkylphospholipid perifosine is internalized by an ATP-dependent translocase activity across the plasma membrane of human KB carcinoma cells[END_REF], with the only modification that, in the MTT-based colorimetric assay to determine perifosine sensitivity, the cell lines were incubated with the anticancer drug for only 2 h instead of 24 h before adding MTT to the culture medium. This modification allowed us to better link the observed cytotoxicity with the amount of perifosine incorporated in a short time by each cell line, and therefore to establish whether cell death actually correlated with intracellular accumulation of perifosine. (GenBank accession no. BC117946), the murine orthologue of the human gene that encodes for the putative APLT ATP8B1 (95% identity), was purchased from Source BioScience Geneservice (Cambridge, UK) as a full open-reading frame (ORF) cloned into pCR-BluntII-TOPO vector. The construct was isolated and purified using the Perfectprep ® plasmid isolation kit (Eppendorf), and the ORF subsequently amplified using the Expand High Fidelity PCR System with the primers: 8b1-F (5'-ACCATGGGCACAGAAAGAGACTCG); 8b1-R (5'-TCACAGATCCTCTTCAGAGATGAGTTTCTGCTCGCTCTCCACTGTGCGCCGGTA ), (myc-tag coding sequence underlined). The resulting PCR product was cloned into pGEMT-Easy from Promega. The nucleotide sequences of all cloned cDNAs were determined automatically, as described previously [START_REF] Lario | Automated laser-induced fluorescence DNA sequencing: equalizing signal-to-noise ratios significantly enhances overall performance[END_REF], and compared against their respective records annotated in GenBank. Unexpectedly, after repeating PCR and isolating new clones, we found that the supplier had provided us with a mutant variant of the Atp8b1 gene (A1954G, which resulted in the amino acid mutation R652G). Since we were concerned about the possible undesired effects that this

Isolation of

Perifosine and NBD-phosphatidylethanolamine (NBD-PE) uptake are linked processes in KB cell lines

In our previous work we generated a perifosine-resistant line (designated as KB PER-R clone 10) from the perifosine-hypersensitive KB Wt cells (with an EC 50 of less than 1 µM for 24 h perifosine exposure). Along with a high resistance to perifosine (grows suitably at 20 µM perifosine), this cell line showed significantly decreased perifosine uptake, a low rate of NBD-aminophospholipid translocation through the plasma membrane, a relatively high exposure of endogenous phosphatidylserine in the cell surface, and lacked NEM-and orthovanadatesensitive, plasma-membrane-associated ATPase activity [START_REF] Muñoz-Martínez | The anti-tumor alkylphospholipid perifosine is internalized by an ATP-dependent translocase activity across the plasma membrane of human KB carcinoma cells[END_REF]. Similar perifosineresistant clones generated in our laboratory from KB Wt cells (clones 4 and 8)

showed intermediate levels of perifosine uptake with respect to KB cells and KB PER-R clone10 cells, along with intermediate rates of NBD-PE translocation through the plasma membrane, and were only moderately resistant to perifosine compared with clone 10 (Fig. 1). This correlation between perifosine uptake and NBD-PE translocation therefore suggested that, as proposed previously [START_REF] Muñoz-Martínez | The anti-tumor alkylphospholipid perifosine is internalized by an ATP-dependent translocase activity across the plasma membrane of human KB carcinoma cells[END_REF], a link could exist between both processes in KB cells. We also tested the uptake rates of BODIPYlactosylceramide (raft-dependent endocytosis marker) and Alexa Fluor 488-dextran 10,000 (fluid phase endocytosis marker) in a subset of fourteen cell lines, and found no significant correlation upon comparing their rates of perifosine uptake (Fig. 2).

A c c e p t e d M a n u s c r i p t 13

CDC50A is implicated in both perifosine and fluorescent aminophospholipid analogues uptake in KB cells

In our previous work we proposed that a putative APLT could be responsible for the internalization of both ALPs and aminophospholipids through the plasma membrane in human KB carcinoma cells. APLTs seems to be a functional complex of at least two proteins (alpha and beta subunits) [START_REF] Pérez-Victoria | Phospholipid translocation and miltefosine potency require both L. donovani miltefosine transporter and the new protein LdRos3 in Leishmania parasites[END_REF][START_REF] Sanchez-Cañete | Low plasma membrane expression of the miltefosine transport complex renders Leishmania braziliensis refractory to the drug[END_REF][START_REF] Saito | Cdc50p, a protein required for polarized growth, associates with the Drs2p Ptype ATPase implicated in phospholipid translocation in Saccharomyces cerevisiae[END_REF], and since fourteen different putative alpha but only three presumed beta subunits (CDC50A, B and C) have been described in humans, we decided to start our study of the potential implication of the latter proteins on perifosine uptake in mammalian cells. Since CDC50C transcripts have been found to terminate aberrantly in humans [START_REF] Osada | Aberrant termination of reproduction-related TMEM30C transcripts in the hominoids[END_REF], only CDC50A and B were likely to play a putative role in perifosine uptake. We therefore generated KB cells expressing CDC50A-V5 or CDC50B-V5 (Fig. 3A) and found that, whereas NBD-aminophospholipid uptake was increased by about 50% when either of these proteins was over-expressed, perifosine uptake was only boosted when CDC50A was over-expressed (Fig. 3B). No significant increase of NBD-PC uptake was detected in KB cells over-expressing CDC50A-V5 or CDC50B-V5. Interestingly, over-expression of CDC50A-V5 in KB PER-R clone10 cells did not increase the rate of perifosine or aminophospholipid uptake (data not shown). Moreover, knock-down of CDC50A and CDC50B reduced perifosine and NBD-aminophospholipid uptake to around 40% of the normal rate when CDC50A expression was lowered, whereas CDC50B silencing only affected NBD-aminophospholipid uptake (Fig. 4A). No significant reduction of NBD-PC was detected after the knock-down experiments.

However, since we had no suitable antibody for monitoring the expression levels of endogenous CDC50A in KB cells, we also performed the same experiments in

CDC50A-V5-transfected KB cells, thereby allowing us to monitor CDC50A-V5

A c c e p t e d M a n u s c r i p t 14 expression levels by Western blot (Fig. 4B). This study showed that lower CDC50A-V5 levels correlated with lower rates of perifosine and NBD-aminophospholipid uptake (Fig. 4C).

CDC50A subcellular localization in KB Wt vs. KB PER-R clone 10 cells

Cdc50-Lem3p proteins play a role as both molecular chaperones needed for the correct folding and export of the alpha subunit P4-ATPase from the endoplasmic reticulum (ER) [START_REF] Puts | Mechanism and significance of P4 ATPase-catalyzed lipid transport: lessons from a Na+/K+-pump[END_REF][START_REF] Paulusma | ATP8B1 requires an accessory protein for endoplasmic reticulum exit and plasma membrane lipid flippase activity[END_REF], and also as catalytically relevant, integral components of the phospholipid translocation machinery [START_REF] Lenoir | Cdc50p plays a vital role in the ATPase reaction cycle of the putative aminophospholipid transporter Drs2p[END_REF]. Since we have found previously that human CDC50A appears to play a role in perifosine uptake in KB cells, we decided to establish the cellular localization of this protein. Thus, we performed indirect immunofluorescence of V5-tagged CDC50A expressed in KB cells and then visualized the samples by confocal microscopy. It is apparent from Fig. 5, which shows the scattered intracellular distribution and co-localization with several markers of different intracellular membrane compartments, that CDC50A-V5 is a membrane protein that is distributed via the endocytic and exocytic pathways in KB cells.

However, we cannot conclude from this figure whether CDC50A-V5 reaches the plasma membrane or not. A comparison of KB and KB PER-R clone10 cells, both of which express similar levels of CDC50A-V5, as monitored by Western blot (not shown), highlighted a striking difference in the intracellular localization of CDC50A-V5. Thus, this protein was expressed well at the ER level, and was also exported to the CIS-Golgi and TRANS-Golgi network (TGN) in both cell lines, whereas its presence in early and late endosomes was significantly lower in KB PER-R cells than in KB cells, thereby suggesting that CDC50A-V5 distribution via the pathways As we were unable to ascertain the presence of CDC50A at the plasma membrane of KB cells by confocal microscopy, we biotinylated the cell surface proteins in CDC50A-V5-transfected KB and KB PER-R cells in parallel. Fig. 6 shows that most of the CDC50A-V5 is intracellular (unbound fraction to streptavidinagarose beads) in KB cells, as also determined by confocal microscopy, but that a detectable amount of CDC50A-V5 reaches the plasma membrane as well. In contrast, no CDC50A-V5 was observable in the fraction of biotinylated cell surface proteins in KB PER-R clone10 cells, although its intracellular levels were comparable to those for KB cells. We used α-tubulin as a reference to compare the amount of CDC50A-V5 expressed in both cell types and as a control of cell membrane integrity throughout the entire biotinylation process. Only cell surface proteins were biotinylated, since α-tubulin is only present in the unbound fraction.

The inability of CDC50A-V5 to reach the plasma membrane in KB PER-R cells is in agreement with its defective membrane trafficking, as suggested previously by confocal microscopy, and further supports our hypothesis that perifosine uptake is very low in this cell line due to the lack of a P4-ATPase APLT at the plasma membrane.

CDC50A and Atp8b1 associate to exit the ER and traffic to the plasma membrane in HeLa and HEK293T cells, where they form a functional perifosine and aminophospholipid translocator

In light of the functional involvement of CDC50A in perifosine uptake in KB cells, we decided to investigate whether CDC50A was also involved in perifosine uptake in A c c e p t e d M a n u s c r i p t other cell lines. Thus, we expressed CDC50A-V5 in HeLa and HEK-293T cells, although with disappointing results: no increase was observed in perifosine or NBDphospholipid uptake in either of these cell lines (see Fig. 8B). However, this finding does not allow us to state that CDC50A plays no role in perifosine uptake in this case, as CDC50A may lack the appropriate alpha subunit partner in these cell lines.

To test this hypothesis, we performed semi-quantitative RT-PCR to compare the mRNA expression levels of CDC50A and ATP8B1 (the best studied P4-ATPase [START_REF] Paulusma | The type 4 subfamily of P-type ATPases, putative aminophospholipid translocases with a role in human disease[END_REF],

which was likely to be responsible for aminophospholipid translocation in these cell lines). This study showed (Fig. 7) that the endogenous expression levels of CDC50A were similar in all three cell lines, and that ATP8B1 was expressed at similar levels in KB and HEK-293T cells, and at even higher levels in HeLa cells (Fig. 7). This finding suggests that a P4-ATPase other than ATP8B1 is likely to be the alpha subunit responsible for the APLT activity observed in KB, HeLa and HEK293T cells.

At this point, instead of simply determining the expression levels of the remaining thirteen putative human alpha subunits, we adopted a more straightforward approach involving the simultaneous transfection of CDC50A-V5 and murine Atp8b1 (a well-studied P4-ATPase, 95% homologous to its human counterpart [START_REF] Paulusma | Atp8b1 deficiency in mice reduces resistance of the canalicular membrane to hydrophobic bile salts and impairs bile salt transport[END_REF]) as a Cterminal myc-tagged protein in HeLa and HEK-293T cells (Fig. 8A). Interestingly, we found that the rates of perifosine uptake in HeLa and HEK-293T cells increased by 60% and 130%, respectively (Fig. 8B). The uptake of NBD-aminophospholipids (but not NBD-PC) was also similarly boosted in both cell lines (Fig. 8B). Confocal microscopy studies revealed that CDC50A-V5 was retained in the intracellular membrane compartments in HeLa cells when over-expressed alone, whereas when it was co-expressed with Atp8b1 in the same cells (Fig. 9), both proteins co-localized at the plasma membrane, in a similar manner to that reported for the pair ATP8B1 + A c c e p t e d M a n u s c r i p t CDC50A [START_REF] Paulusma | ATP8B1 requires an accessory protein for endoplasmic reticulum exit and plasma membrane lipid flippase activity[END_REF]. Co-localization of these two proteins was also observed when coexpressed in HEK-293T cells (data not shown).

DISCUSSION

A detailed knowledge of the molecular mechanism of action of anticancer drugs will help to identify those patients who are most likely to respond to a given treatment based on such drugs, to overcome drug resistance phenomena and to improve the final treatment outcome. Indeed, there seems to be a general consensus amongst researchers regarding the manner in which perifosine causes cell death in several different types of cancer cells. However, some discrepancies have arisen concerning how perifosine enters cancer cells, and what previous steps are required for the drug to reach its intracellular targets and trigger cell death. Thus, it has been proposed that ALP enters cells either by raft-dependent endocytosis [START_REF] Van Der Luit | A new class of anticancer alkylphospholipids uses lipid rafts as membrane gateways to induce apoptosis in lymphoma cells[END_REF][START_REF] Van Der Luit | Alkyllysophospholipid accumulates in lipid rafts and induces apoptosis via raftdependent endocytosis and inhibition of phosphatidylcholine synthesis[END_REF][START_REF] Van Der Luit | Different modes of internalization of apoptotic alkyl-lysophospholipid and cell-rescuing lysophosphatidylcholine[END_REF] or by transporter-mediated translocation of perifosine through the plasma membrane [START_REF] Muñoz-Martínez | The anti-tumor alkylphospholipid perifosine is internalized by an ATP-dependent translocase activity across the plasma membrane of human KB carcinoma cells[END_REF][START_REF] Vink | Lipid rafts and metabolic energy differentially determine uptake of anti-cancer alkylphospholipids in lymphoma versus carcinoma cells[END_REF][START_REF] Menez | Inward translocation of the phospholipid analogue miltefosine across Caco-2 cell membranes exhibits characteristics of a carrier-mediated process[END_REF].

We firmly support the latter hypothesis, and in this present work, which is a continuation of a previous study [START_REF] Muñoz-Martínez | The anti-tumor alkylphospholipid perifosine is internalized by an ATP-dependent translocase activity across the plasma membrane of human KB carcinoma cells[END_REF], we present strong evidence to suggest that perifosine enters cancer cells via APLT(s) at the plasma membrane, and that CDC50A plays a key role in this process. Thus, we have found a positive correlation between perifosine uptake and NBD-aminophospholipid translocation (a wellestablished marker for aminophospholipid translocase activity) in several KB clones obtained in our laboratory. Furthermore, we found no correlation whatsoever between perifosine uptake and endocytosis rates in a subset of fourteen cell lines, which is in stark contrast with previous works where the role of raft-dependent endocytosis in perifosine uptake has been proposed [START_REF] Van Der Luit | A new class of anticancer alkylphospholipids uses lipid rafts as membrane gateways to induce apoptosis in lymphoma cells[END_REF][START_REF] Van Der Luit | Alkyllysophospholipid accumulates in lipid rafts and induces apoptosis via raftdependent endocytosis and inhibition of phosphatidylcholine synthesis[END_REF][START_REF] Van Der Luit | Different modes of internalization of apoptotic alkyl-lysophospholipid and cell-rescuing lysophosphatidylcholine[END_REF]. This discrepancy may, however, only be apparent if we consider that only adherent cell lines, but no leukemia or lymphoma cells, were used in the present study, whereas the authors who proposed the raft-dependent hypothesis for perifosine uptake have recently published a further study [START_REF] Vink | Lipid rafts and metabolic energy differentially determine uptake of anti-cancer alkylphospholipids in lymphoma versus carcinoma cells[END_REF] where they discriminate between lymphoma and carcinoma cells on the basis of their prevalent perifosine-uptake mechanism (raftdependent endocytosis in the former, and a plasma membrane translocator in the latter). It is also possible that both mechanisms of perifosine uptake could co-exist in the same cell, with their relative importance varying with cell type. Thus, for example, KB cells accumulate perifosine very rapidly, and independently of rafts endocytosis, due to their high APLT activity [START_REF] Muñoz-Martínez | The anti-tumor alkylphospholipid perifosine is internalized by an ATP-dependent translocase activity across the plasma membrane of human KB carcinoma cells[END_REF][START_REF] Vink | Lipid rafts and metabolic energy differentially determine uptake of anti-cancer alkylphospholipids in lymphoma versus carcinoma cells[END_REF], whereas HeLa cells accumulate perifosine five times slower as they exhibit lower APLT activity. In this case, the slow drug uptake would depend mainly on endocytosis, as suggested previously [START_REF] Van Der Luit | Different modes of internalization of apoptotic alkyl-lysophospholipid and cell-rescuing lysophosphatidylcholine[END_REF].

We also found that modulating the CDC50A expression levels in KB cells markedly affected the rate of perifosine and aminophospholipid uptake, thereby suggesting that a putative APLT is responsible for translocation of both types of molecules through the plasma membrane. It also appears that the CDC50A beta subunit should be the limiting factor for the APLT complex in KB cells, where it is distributed throughout the endo-membranous compartments, trafficking from the ER to the plasma membrane of KB cells. In contrast, CDC50A exports from the ER only to the TGN in KB PER-R clone10 cells and does not reach the cell surface, thereby suggesting that these perifosine-resistant cells harbor a defect that results in the incorrect trafficking of CDC50A at this point of the exocytic pathway. Since CDC50A must couple with a P4-ATPase for correct trafficking of the full APLT complex to the plasma membrane [START_REF] Paulusma | ATP8B1 requires an accessory protein for endoplasmic reticulum exit and plasma membrane lipid flippase activity[END_REF], the observed traffic defect of CDC50A in KB PER-R cells could explain the inability of this cell line to translocate perifosine. In contrast to KB A c c e p t e d M a n u s c r i p t cells, over-expression of CDC50A in the human cell lines HeLa and HEK-293T did not modulate perifosine uptake, thus suggesting that, in this case, the alpha subunit rather than CDC50A may be the limiting component of the putative APLT. Human ATP8B1 is a P4-ATPase with APLT activity when complexed with CDC50A, although the expression levels of both CDC50A and ATP8B1 in HeLa and HEK-293T cells were found to be very similar to those in KB. Therefore, if we assume that no differential regulation of the APLT complex occurs in each cell line, and that the overall APLT activity depends mainly on the relative expression levels of its two components, it seems likely that ATP8B1 is not the alpha subunit responsible for the different APLT activities measured in these cell lines. However, we found that the uptake of perifosine and NBD-aminophospholipids (but not NBD-PC) increased remarkably when both HeLa and HEK-293T cells were co-transfected with Atp8b1myc and CDC50A-V5. Moreover, CDC50A-V5 remained intracellular when overexpressed alone, but when co-expressed with Atp8b1 both proteins exported from the ER and co-localized at the plasma membrane, as is the case with human ATP8B1 + CDC50A [START_REF] Paulusma | ATP8B1 requires an accessory protein for endoplasmic reticulum exit and plasma membrane lipid flippase activity[END_REF]. CDC50A may be coupled with a P4 ATPase to drive perifosine and aminophospholipids uptake and could contribute to the transport specificity of the complex.

Taken together, the results of this work strongly support the role of an APLT as the means of entry of the anticancer drug perifosine through the plasma membrane in adherent carcinoma cells. Such an APLT in humans is most likely formed by the beta subunit CDC50A plus a P4-ATPase alpha subunit(s) whose identity remains to be established.

A c c e p t e d M a n u s c r i p t 4 cells [ 14

 414 -16], have been described. A clearer insight into the molecular determinants of perifosine entry into cancer cells would therefore be of enormous importance for the identification of additional biomarkers for the efficacy of perifosine as an anticancer drug.

. 1 .

 1 Chemicals. Perifosine (octadecyl-[N,N-dimethyl-piperidinio-4-yl]-phosphate), and [ 14 C]-perifosine (30.9 mCi/mmol) were provided by Zentaris GmbH (Frankfurt, Germany). Lipofectamine 2000, Anti-V5 monoclonal antibody, Alexa Fluor 633conjugated anti-mouse IgG, pEF6-V5/His and pcDNA3.1-HYGRO vectors, FITCconcanavalin A, BODIPY-lactosylceramide, and Alexa Fluor 488-labeled dextranA c c e p t e d M a n u s c r i p t 5 10,000 were purchased from Invitrogen (Prat de Llobregat, Barcelona, Spain). Predesigned ON-TARGETplus siRNA for human CDC50A and CDC50B were purchased from Dharmacon (Lafayette, CO). Fluorescent (NBD-labeled) phospholipid analogues 1-palmitoyl-2-{6-[(7-nitro-2-(1,3-benzoxadiazol)-4yl)amino]hexanoyl}-sn-glycero-3-phosphoethanolamine (NBD-PE), -phosphoserine (NBD-PS) and -phosphocholine (NBD-PC) were purchased from Avanti Polar Lipids (Alabaster, AL). Anti-rabbit IgG monoclonal antibody, monoclonal anti-α-tubulin antibody, and the rabbit polyclonal antibodies anti-myc, anti-calnexin, anti-GM130, anti-TGN46, anti-Rab5 and anti-Rab7, and MTT [3-(4,5-Dimethylthiazol-2-yl)-2,5diphenyltetrazolium bromide] were purchased from Sigma-Aldrich (Madrid, Spain). Fugene HD, High Pure RNA isolation kit and Expand High Fidelity PCR System were purchased from Roche (San Cugat del Vallés, Barcelona, Spain), and the High-capacity cDNA Reverse transcription kit plus RNAse inhibitor from Applied Biosystems (Alcobendas, Madrid, Spain). The Perfectprep ® plasmid isolation kit from Eppendorf Ibérica (San Sebastian de los Reyes, Madrid, Spain) was used to isolate purified plasmids. QuickChange II XL Site-directed Mutagenesis kit was purchased from Stratagene-Cultek SLU (Madrid, Spain), and the pGEM-T Easy vector from Promega Biotech Ibérica (Alcobendas, Madrid, Spain). All other chemicals were purchased from Sigma-Aldrich (Madrid, Spain).

2. 2 .

 2 Cell lines and cultures. The mammalian cell lines described in the present work were grown at 37 °C in a humidified atmosphere and 5% CO 2 in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% (v/v) fetal bovine serum (FBS), 2 mM Glutamax I, 1:50 Micokill AB solution 50X (PAA Laboratories GmbH) and 250 U/ml penicillin G plus 250 μg/ml streptomycin sulphate (Invitrogen), unless indicated otherwise. The following cells were used: murine fibroblast embryo NIH A c c e p t e d M a n u s c r i p t 6 3T3 cells (ATCC number CCL-1658); chinese hamster ovary CHO-K1 cells (ATCC number CCL-61; cultured in F-12K medium); canine kidney MDCKII cells [22]; monkey epithelial kidney Vero and COS-7 cells (ATCC CRL-1586 and CRL-1651); human epidermal carcinoma KB cells (wild-type) (ATCC number CCL-17) and its perifosine-resistant counterparts KB PER-R clone 4, 8 and 10 [14]; human breast carcinoma MCF-7 and MDA-MB-23 cells [23, 24]; human cervix epithelial adenocarcinoma HeLa cells (ATCC CCL-2); human epithelial embryonic kidney HEK-293T cells (ATCC CRL-11268); human pancreas carcinoma AsPC1 and MiaPaca2 cells (ATCC CRL-1682 and CRL-1420); human 2008 ovary cells [25], kindly provided by Prof. Piet Borst (Division of Molecular Biology, The Netherlands Cancer Institute, Amsterdam, The Netherlands); and human osteosarcoma HOS cells (ATCC CRL-1543).

  CDC50A, CDC50B and Atp8b1, and DNA constructs. The human cDNAs of CDC50A (GenBank accession no. NM_018247, longer transcription variant 1) and CDC50B (GenBank accession no. NM_001017970) A c c e p t e d M a n u s c r i p t 7 genes, which code for the putative beta subunits of the human APLT(s) (P-type ATPase P4 subfamily), were isolated by RT-PCR with the High Capacity cDNA Reverse transcription kit plus RNAse inhibitor (Applied Biosystems), using total RNA obtained from KB cells with the High Pure RNA isolation kit (Roche) as template; [primers: CDC50A-F (5'-ACCATGGCGATGAACTATAACG), and CDC50A-R (5'-TTAAATGGTAATGTCAGCTGTATTAC); CDC50B-F (5'-ACCATGACCTGGAGCGCCACGGC) and CDC50B-R (5'-TCACTCCTCGTCGTCGTCCTGG)]. The resulting PCR products were cloned into the pEF6-V5/His TOPO TA expression vector from Invitrogen, in frame with the Cterminal V5 epitope and the polyhistidine (6 x His) tag. The cDNA of Atp8b1

A c c e p t e d M a n u s c r i p t 8 mutation 2 . 7 .

 827 could cause in the phenotypic function of the resulting protein, we decided to fix it by using the QuickChange II XL Site-directed Mutagenesis kit from Stratagene using the primers (fixing point-mutation is underlined): FAtp8b1-mut (5'-AAGTGAGACTCTCAGGACCCTGTGCCTGTG); RAtp8b1-mut (5'-CACAGGCACAGGGTCCTGAGAGTCTCACTTG). Once this annoying problem had been resolved, we re-amplified the Atp8b1 ORF with the primers FAtp8b1-NheI (5'-TTATTGCTAGCACCATGGGCACAGAAAGAGACTCG); RAtp8b1-BamHI (5'-TTATTGGATCCTCAAGATCTCTTCAGAGATGAG), and cloned the resulting PCR product as an NheI/BamHI fragment into the pcDNA3.1-HYGRO expression vector.2.5. Cell transfection.KB Wt, KB PER-R, HeLa and HEK-293T cells were stably transfected with pEF6-CDC50A-V5 and pEF6-CDC50B-V5 using Lipofectamine 2000 following the manufacturer's protocol. The stable transfectants were selected with 10 µg/ml blasticidin in only one week, and the expression levels of both proteins readily visualized by Western blotting or indirect immunofluorescence (IFI) under a microscope using monoclonal anti-V5 1:5000 from Invitrogen. Untransfected and CDC50A-V5-transfected KB Wt and KB PER-R cells were refractory to both stable and transient transfection with the pcDNA3.1-HYGRO-Atp8b1-myc construct using either Lipofectamine 2000 or Fugene HD, as judged by the fact that no Atp8b1-myc expression was observed, either by Western blot or by IFI using several anti-myc antibodies. In contrast, HeLa and HEK-293T cells could be transiently transfected with CDC50A-V5 alone or in combination with Atp8b1-myc, with the best results being obtained using Fugene HD. 2.6. siRNA-based CDC50A and CDC50B knockdown. Human CDC50A and CDC50B were targeted using SMARTpool siRNA oligonucleotides (Dharmacon),A c c e p t e d M a n u s c r i p t following a previously described procedure [14], at concentrations of 1, 10 and 100 nM of siRNA oligonucleotides in order to monitor the concentration dependency of the gene silencing. Biotinylation of cell surface proteins. KB Wt and PER-R cells were washed with PBS and incubated on ice with 1 mM EZ-Link Sulfo-NHS-SS-biotin (Thermo Scientific Pierce, Madrid, Spain) in 1 ml PBS for 2 h at 4 °C. The biotinylation reaction was quenched by washing the cells three times with 50 mM Tris-HCl, pH 7.4, then cells were incubated on ice for 30 min in 100 µl lysis buffer (50 mM Tris-HCl pH 7.4, 150 mM NaCl, 1% Nonidet P-40) plus a protease inhibitor cocktail (Sigma-Aldrich). The cell lysate was clarified by centrifugation at 21,000 x g for 10 min at 4 °C, and the supernatant incubated for 2 h at 4 ºC in continuous agitation with 100 µl of packed streptavidin-agarose beads (Thermo Scientific Pierce) pre-washed with lysis buffer. The agarose beads were then washed three times with 1 ml of PBS + 1 M NaCl and two times with 1 ml of PBS, and the biotinylated proteins eluted in 30 µl of 2 x Laemmli buffer. 2.8. Immunoblotting. Protein samples were fractionated by SDS-PAGE under standard conditions and electrotransferred onto Immobilon P membranes (Millipore, Bedford, MA). Immunodetection of the V5-tagged beta subunits and Atp8b1-myc was performed with 1:5000 monoclonal anti-V5 and 1:1000 anti-myc polyclonal antibody in buffer A (PBS containing 0.01% Tween 20 and 0.1% BSA). After washing, the membranes were incubated with horseradish peroxidase-conjugated secondary rabbit anti-mouse or goat anti-rabbit (1:5000) immunoglobulin G (Dako, Barcelona, Spain) in buffer A. Immunodetection of α-tubulin for sample normalization A c c e p t e d M a n u s c r i p t 10 was performed with monoclonal anti-α-tubulin antibody (Sigma-Aldrich) at 1:5000. Signals were detected using either the ECL or the ECL Plus chemiluminescent substrate (Thermo Scientific Pierce).2.9. Immunofluorescence and confocal microscopy. Cells were seeded onto cover slip circles (diameter: 12 mm; thickness: 0.13-0.16 mm, from Menzzel-Glässer) in 24-well plates (Thermo Scientific Nunc, Madrid, Spain) and cultured under standard conditions. On the day of the experiment, the cells were washed three times with PBS and fixed with 2% para-formaldehyde in PBS for 20 min at room temperature. After quenching with 10 mM Tris-HCl pH 7.4 and washing twice with the same buffer, the cells were permeabilized with a 0.1% solution of Triton X-100 detergent (Sigma-Aldrich) in PBS for 10 min at room temperature. They were then washed with PBS and incubated for 30 min with PBS 1% BSA to block the nonspecific binding of antibodies. The primary antibodies anti-V5 (1:5000), anti-myc (1:1000), anti-calnexin (1:200), anti-GM130 (1:2500), anti-TGN46 (1:200), anti-Rab5 (1:200) and anti-Rab7 (1:200) were incubated at 37 °C in a humidified atmosphere for 1-2 h in PBS 1% BSA. The cells were then rinsed three times with PBS and treated with the secondary antibodies FITC-conjugated anti-rabbit IgG antibody (1:1000-1:2000) from Sigma-Aldrich, and AlexaFluor 633-conjugated anti-mouse IgG antibody (1:500) from Invitrogen in PBS 1% BSA. After the labelling procedure was complete, the cells were washed three times with PBS and mounted onto glass slides using Prolong Gold mounting medium with DAPI (Invitrogen). Confocal images were acquired using a Leica TCS-SP5 confocal scanning laser microscope equipped with Ar and He-Ne lasers and connected to a Leica DMIRBE inverted epi-fluorescence microscope. The following filter settings were used: DAPI A c c e p t e d M a n u s c r i p t 11 (blue): excitation at 405 nm (405 diode laser, UV), emission at 409-469 nm; FITC (green): excitation at 458 nm (argon laser, 20% power), emission at 494-603 nm; Alexa Fluor 633 (red): excitation at 633nm (633 HeNe laser, 20% power), emission at 648-692 nm. Serial fluorescence images were obtained using the 63× objective and processed and assembled with LAS AF Leica confocal software. Image analysis and subcellular colocalization fluorograms were generated and analyzed using the Leica confocal software package. 2.10. Semi-quantitative RT-PCRs. Semi-quantitative RT-PCRs were performed using the total cDNAs synthesized above as templates to establish the endogenous expression of CDC50A and ATP8B1 genes in KB, HeLa and HEK293T cells. Standard PCRs were performed using serial dilutions of total cDNAs as templates and highly specific primers [ATP8B1 forward (5'-CAG ACT GCA TAC GAG GAT TGG); ATP8B1 reverse (5'-CCT TCA ACC GCT TGC GAT G); CDC50A forward (5'-CTA TAA CGC GAA GGA TGA AGT G); CDC50A reverse (5'-TCC ATA ATA CAT AAA CAC GTT GCC); β-actin forward (5'-CCT CAT GAA GAT CCT CAC CG); βactin reverse (5'-CAG TGA TCT CCT TCT GCA TCC)] to amplify internal 631-, 357and 660-bp sequences from the ATP8B1, CDC50A and β-actin cDNAs, respectively. The PCR conditions for each step were as follows: Initial denaturation (95 ºC, 5 min); 35 cycles (denaturation: 92 ºC, 30 s; annealing: 52 ºC, 45 s; extension: 72 ºC, 45 s); final extension: 72 ºC, 3 min. The resulting PCR products were run in 0.8% agarose gels stained with ethidium bromide and visualized in a UV transiluminator.

A

  c c e p t e d M a n u s c r i p t connecting TGN, early endosomes, the plasma membrane and late endosomes is hindered in KB PER-R cells.

FIGURE LEGENDS Figure 1 .

 LEGENDS1 FIGURE LEGENDS

Figure 2 .

 2 Figure 2. Comparative analysis of perifosine and endocytosis markers in different

Figure 3 .

 3 Figure 3. Uptake of perifosine and NBD-phospholipids in KB Wt cells over-

Figure 4 .

 4 Figure 4. Silencing of CDC50A and CDC50B in KB Wt cells. A) Uptake of NBD-

A c c e p t e d M a n u s c r i p t 27 The

 27 absolute uptake values measured in KB Wt cells with Dharmafect 1 in this set of experiments were: perifosine: 3657 ± 162 pmol/mg protein; NBD-phospholipids uptake expressed in flow cytometry mean channel values after substracting the cell autofluorescence: NBD-PE: 609 ± 71; NBD-PS: 505 ± 63; NBD-PC: 521 ± 66.

Figure 5 .

 5 Figure 5. CDC50A subcellular localization by confocal microscopy analysis. KB Wt

Figure 6 .

 6 Figure 6. Cell-surface expression of CDC50A. KB Wt and KB PER-R clone10 cells

Figure 7 .

 7 Figure 7. CDC50A and ATP8B1 mRNA levels in different cell lines. Semi-quantitative

Figure 8 .Figure 9 .

 89 Figure 8. Uptake of NBD-phospholipids and perifosine in cells over-expressing both

Figure 1 A

 1 Figure 1

Figure 2

 2 Figure 2

Figure 3 A

 3 Figure 3

Figure 6 Page

 6 Figure 6

Figure 7 A

 7 Figure 7

Figure 9 Page

 9 Figure 9

ACKNOWLEDGEMENTS

This work was supported by "Plan Andaluz de Investigación, Proyecto de Excelencia (CTS 187)", Junta de Andalucia. We acknowledge Aeterna Zentaris (Frankfurt, Germany) for the financial support and for providing the perifosine and [ 14 C]-perifosine used in this study. We want to thank Carlos Castillo Pérez for his help in the fixing and cloning of Atp8b1 cDNA and Dr. Luis Carvalho for his technical assistance in RT-PCR experiments.