Qin Zhang

Mitsuhiro Hayashibe

David Guiraud

Muscle Fatigue Tracking based on Stimulus Evoked EMG and Adaptive Torque Prediction

Functional electrical stimulation (FES) is effective to restore movement in spinal cord injured (SCI) subjects. Unfortunately, muscle fatigue constrains the application of FES so that output torque feedback is interesting for fatigue compensation. Whereas, inadequacy of torque sensors is another challenge for FES control. Torque estimation is thereby essential in fatigue tracking task for practical FES employment. In this work, the Hammstein cascade with electromyography (EMG) as input is applied to model the myoelectrical mechanical behavior of the stimulated muscle. Kalman filter with forgetting factor is presented to estimate the muscle model and track fatigue. Fatigue inducing protocol was conducted on three SCI subjects through surface electrical stimulation. Assessment in simulation and with experimental data reveals that the muscle model properly fits the muscle behavior well. Moreover, the time-varying parameters tracking performance in simulation is efficient such that real time tracking is feasible with Kalman filter. The fatigue tracking with experimental data further demonstrates that the proposed method is suitable for fatigue tracking as well as adaptive torque prediction at different prediction horizons.

I. INTRODUCTION

Functional electrical stimulation (FES) is one of the existing solution to partly restore the lost motor function in persons with spinal cord injury (SCI). The electrical stimulus can artificially generate action potential on the axons of the alpha motor neurons to drive the muscle contraction instead of central nervous system. To date, FES has been used over a wide range of FES-aided support of standing, gait, grasp, or foot drop correction and tremor compensation.

In these applications, the stimulation is required to provide strong, consistent muscle force. It is well known, however, that muscles get fatigued more rapidly when artificially stimulated than when excited by the central nervous system due to the way which the motor units are recruited: inverse size principle, synchronized activation of motor units and constant order of recruitment. As a result, successful implementation of FES paradigms for rehabilitation has been greatly constrained by the early occurrence of fatigue. In addition, since SCI individuals have also lost sensory pathways, they cannot perceive their muscle fatigue as stimulation proceeds. This leads to movement failures and sub-optimal FES parameters tuning in practical applications.

Several techniques have been developed for muscle fatigue tracking. A fitness function was introduced to describe the effect of muscle fatigue on shank motion prediction induced Q. Zhang, M. Hayashibe, and D. Guiraud are with DEMAR team, INRIA Sophia Antipolis, and LIRMM, UMR 5506 CNRS UM2, 161 Rue Ada, 34095 Montpellier Cedex 5, France. {qin.zhang, hayashibe, guiraud}@lirmm.fr by FES [START_REF] Riener | Biomechanical Model of the Human Knee Evaluated by Neuromuscular Stimulation[END_REF]. Muscle fatigue was modeled based on the metabolic profiles to represent force reduction and recovery in intermittent stimulation [START_REF] Mizrahi | Predicted and Measured Muscle Forces After Recoveries of Differing Durations Following Fatigue in Functional Electrical Stimulation[END_REF], [START_REF] Lim | Model of Activation Dynamics for an FES-Induced Muscle Fatigue[END_REF]. A fatigue model, coupled with the force model that predicts the fatigue induced by different stimulation patterns on different days during isometric contractions was reported [START_REF] Ding | A predictive model of fatigue in human skeletal muscles[END_REF]. However, due to high nonlinearity and dynamic complexities, these muscle models are difficult to implement in application and it remains difficult to estimate the parameters.

In order to take into account the effect of muscle fatigue for FES control, reliable sensors are required to produce desired feedback or feed-forward signals. However, torque measurement through external sensors might be affected by slight voluntary motion, recording missing or other disturbances, and the devices, such as torque measuring chair, are not convenient at all for patients' daily use. The implanted sensor is a possible technique [START_REF] Hoffer | Closed-Loop, Implanted-Sensor, Functional Electrical Stimulation System for Partial Restoration of Motor Functions[END_REF] but not available for practical use yet. As electromyography (EMG) signal was reported highly correlated with FES-induced torque in different muscle condition [START_REF] Chesler | Suface EMG as a Fatigue Indicator During FES-induced Isometric Muscle Contractions[END_REF] [START_REF] Erfanian | Excitation-Contraction Fatigue During Sustained Electrical Stimulation of Paralyzed Muscle[END_REF], it was demonstrated able to be used as fatigue indicator. A predictive model using EMG instead of electrical stimulation as model input was developed in [START_REF] Erfanian | Using Evoked EMG as a Synthetic Force Sensor of Isometric Electrically Stimulated Muscle[END_REF], and it was underlined evoked EMG could be used as a synthetic torque sensor. In [START_REF] Tepavac | Detection and Prediction of FES-Induced Fatigue[END_REF], it was suggested to combine time and frequency domain variables for better prediction of force. Embedded in these works, a fixed relationship between EMG and force/torque during sustained stimulation was assumed. Nevertheless, some researchers found nonlinear relationship and dissociation between myoelectrical and mechanical activity in different fatigue states [START_REF] Mizrahi | Myoelectric and Force Characteristics in Transcutanteous Isometric FES[END_REF] or in the recovery process [START_REF] Yu | Mechanical and Electromyographic Response to Stimulated Contractions in Paralyzed Tibialis Anterior Post Fatiguing Stimulations[END_REF]. The evidence figures out the limitation of fatigue prediction strategies based on models with fixed parameters. Therefore, it is essential to find a suitable model in order to properly represent various muscle behavior.

The present work focuses on developing a method to predict the FES-induced mechanical contraction behavior during different muscle fatigue states. Online estimation and fatigue prediction are presented, based on a time-varying myoelectrical mechanical model, where the parameters are estimated by Kalman filter (KF) with forgetting factor. In terms of the improvements in fatigue tracking, the proposed method can be directly applied for fatigue prediction when perfect sensor data are available. Furthermore, if the measurements of the sensors are affected by external disturbances, the method can bridge the disturbances and provide sufficient fatigue tracking accuracy.

II. MYOELECTRICAL MECHANICAL MODEL OF ELECTRICALLY STIMULATED MUSCLE

A. Model Structure

The discrete-time Hammerstein cascade is used to model electrically stimulated muscle, relating FES-evoked EMG to induced torque under isometric condition. This model structure has ever been used to represent some highly nonlinear systems, for example, biomechanical systems, such as stretch reflex [START_REF] Dempsey | Identification of Hammerstein Models With Cubic Spline Nonlinearities[END_REF] and electrically stimulated muscle, but relating stimulation to muscle force under isometric conditions [START_REF] Chia | Recursive Parameter Identificaiton of Constrained Systems: An Applicaiton to Electrically Stimulated Muscle[END_REF]. It consists of a memoryless nonlinear part followed by a linear dynamic part as shown in Fig. 1. In this work, the linear part is chosen as an autoregressive with external input (ARX) model, which has been shown experimentally to yield good prediction of output force/ torque in isometric situation [START_REF] Bernotas | A Discrete-Time Model of Electrically Stimulated Muscle[END_REF]. Given model order (l, m), it can be described in operator notation:

A(z -1 )y(t) = B(z -1 )h(t) (1)
where y(t) is the muscle torque at time instant t, h(t) is the activation level of muscle. A(z) and B(z) are the polynomials with respect to the backward shift operator z -1 , defined by the following equations:

A(z -1 ) = 1 + a 1 z -1 + a 2 z -2 + • • • + a l z -l B(z -1 ) = b 1 z -1 + b 2 z -2 + • • • + b m z -m (2)
The memoryless nonlinear part maps the system input, i.e., instantaneous EMG signal u(t), to h(t). Traditionally, it is modeled with a nth order polynomial, in which case,

h(t) = n i=1 γ i u i (t) (3)
Substituting ( 2) and (3) into (1) and expanding, the generated output at a given instant t can be parameterized as:

y(t, g) = l i=1 a i y(t -i) + m i=1 n j=1 μ ij (u(t -i)) j ( 4
)
where

μ ij = b i γ j , g = [a 1 , • • • , a l , μ 11 , • • • μ mn ] T is a
vector containing the muscle model coefficients. At a given instant t, the computation of k-step future output estimations using ( 4) is performed by assuming the system is stationary during the prediction horizon. The elements of g as well as the torque versus the EMG features are time-varying due to the effects of fatigue and the associated biochemistry. Moreover, model reidentification can effectively improve prediction performance in different muscle states [START_REF] Zhang | Torque Prediction Using Stimulus Evoked EMG and its Identification for Different Muscle Fatigue States in SCI Subjects[END_REF], which inspires us to employ Kalman filter to estimate the timevarying model parameters online, and improve torque prediction for fatigue tracking.

B. State-Space Representation of Time-Varying Model

For the use of KF algorithm, a model representation in state space form is required. Such a model consists of process equation and measurement equation. Generally, for a model (l, m, n), its state-space model can be written in a compact form:

1) process equation

x k = f (x k-1 , u k-1 ) = Ax k-1 + n i=1 B i (u k-1 ) i (5)
2) measurement equation

y k = Cx k (6) In (5), u is model input, x = x 1 , x 2 , • • • , x q T is called
state vector, q = max{l, m}. The q × q matrix A relates the state at the previous time step x k-1 , to the state at the current step x k . The q × 1 matrix B i relates the model input at the previous time step

(u k-1 ) i , i = 1, 2, • • • , n,
to the state at the current step x k . They are represented as

A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a 1 1 0 • • • 0 0 a 2 0 1 • • • 0 0 . . . . . . . . . . . . . . . . . . a q-1 0 0 • • • 0 1 a q 0 0 • • • 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , B i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ μ 1,i μ 2,i . . . μ q-1,i μ q,i ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
In practice, the matrices A and B i might change with each time step due to time-varying property of muscle fatigue. The y in ( 6) is the measurement of output torque. The 1 × q matrix C relates the state at the current step x k , to the measurement at the current step y k with the following expression:

C = 1 0 • • • 0 0

III. IDENTIFICATION OF THE MYOELECTRICAL MECHANICAL MUSCLE MODEL

In this study, a method to predict FES-induced torque only driven by EMG was developed, where KF with forgetting factor was employed to estimate the model states and parameters. The EMG and ankle torque will be used as model input and output for estimation.

A. Kalman Filter with Forgetting Factor

For online estimation and prediction, a KF was employed for the recursive estimation of the model parameters. The KF is an efficient recursive filter that estimates the internal states and parameters of a discrete-time system from a series of noisy measurements. It has some advantages in explicit modeling of the uncertainties associate with the proposed model and output measurements. Online estimation of the states in x and model parameters in g is investigated simultaneously. Parameter estimation with KF is performed considering the unknown parameters as elements of the state vector. In this way, the basic KF algorithm does not need to be changed except that the state vector x will be augmented with the unknown parameters. That is, the meta-state vector w has the expression w = x T g T T . The parameters in g are assumed time-invariant or slowly varying comparing to the process, that is, g k = g k-1 . The augmented system is then described by

w k = F(w k-1 , u k-1 ) = f (x k-1 , u k-1 ) g k , ( 7
)
y k = Hw k = [C 0 1×(l+m * n) ]w k (8)
The recursive estimation of the state-space model using KF consists of two phases, prediction and correction. The main equations are given as follows. More details can be found in [START_REF] Welch | An introduction to the kalman filter[END_REF].

1) prediction phase:

ŵ- k = F( ŵk-1 , u k-1 ) (9)
P - k = D k P k-1 D T k + Q k-1 (10)
2) correction phase:

K k = P - k H T k (H k P - k H T k + R k ) -1 (11) ŵk = ŵ- k + K k (y k -H ŵ- k ) (12) P k = (I -K k H k )P - k (13)
Equations ( 9), [START_REF] Mizrahi | Myoelectric and Force Characteristics in Transcutanteous Isometric FES[END_REF] project the state and error covariance estimates forward from time step k -1 to step k. The ŵk-1 and P k-1 are initial estimates for the state and for the measurement variance, respectively. R k is measurement noise covariance, while Q k is process noise covariance. D k is the Jacobian matrix of the plant transfer functions with respect to the involved variables at step k, with each element D [i,j] computed as:

D [i,j] = ∂F [i] ∂w [j] ( ŵk-1 , u k-1 )
In correction phase, equations ( 11) -( 13) adjust the projected estimates by an actual measurement at step k. K k is called Kalman gain.

Although KF is an effective way to estimate the state and parameters of a discrete-time controlled process. Its performance in estimating the time-varying parameters is degraded since it refers to the entire history of the past measurements [START_REF] Mack | Speech Parameter Estimation by Timeweighted-Error Kalman Filter[END_REF]. This is particularly troublesome since the myoelectrical and mechanical activity of stimulated muscle may vary with different fatigue states in prolonged or repetitive stimulation. Deliberately, forgetting factor λ is incorporated to allow tracking of time-varying parameters [START_REF] Mack | Speech Parameter Estimation by Timeweighted-Error Kalman Filter[END_REF]. Thus, equation ( 10) and ( 11) can be rewritten as

P - k = D k P k-1 D T k /λ (14) K k = P - k H T k (H k P - k H T k + λ) -1 (15)
Choosing λ ∈ [ 0, 1] is based on how much we expect the filter to forget the past measurements. In most application, λ is only slightly lower than 1 (for example, 0.85 ∼ 1). In our study, λ is chosen to be 0.97 in simulation and 0.997 in experimental validation. These values represent a compromise between the smoothness of tracking and lag in detecting the changes in model parameters or muscle condition.

B. Experimental Setup for Model Identification

In order to investigate the estimation and prediction method proposed in this paper, experiments were conducted on SCI subjects in PROPARA rehabilitation center, Montpellier, France. The experimental setup, stimulation and recording electrodes configuration are depicted in Fig. 2. The subjects were seated on the chair with their ankle at 90 o , while the foot was strapped on the pedal. Amplitude modulated electrical current pulses were sent to the the right triceps surae muscle via surface electrodes (10cm×3cm) and caused muscle contractions and ankle-joint motion in isometric condition. One electrode was placed 5cm beneath the popliteal cavity and the other beneath the insertion point of the medial and lateral gastrocnemius on the Achilles tendon. The muscle group was stimulated with amplitude modulation at a constant frequency (30Hz) and constant pulse-width (450μs), under isometric conditions by a portable stimulator (Cefar Physio 4, Cefar Medical AB, Lund, Sweden).

Evoked EMG activity of soleus were recorded, amplified (gain 1000) and sampled at 4KHz by an acquisition system (Biopac MP100, Biopac Systems Inc., Santa Barbara, CA, USA). Two bipolar AgCl EMG electrodes were positioned over muscle belly along the muscle fiber direction with 20mm interelectrode spacing. The reference electrode was placed on the patella of another leg. The skin under the electrodes was shaved to minimize the impedance. Isometric ankle plantar-flexion torque was measured using a calibrated dynamometer (Biodex, Biodex Medical Systems, Inc., New York, USA), interfaced with the acquisition system (Biopac MP100) and sampled at 4KHz. The detailed experimental protocol has been shown previously [START_REF] Zhang | Torque Prediction Using Stimulus Evoked EMG and its Identification for Different Muscle Fatigue States in SCI Subjects[END_REF]. The data used for this work include intermittent four series of five stimulation trials which were enforced to induce muscle fatigue. Each trial consists of 1s ramp, 2s plateau at intensity evoking 50% of maximum torque, 1s falling and 2s interval between adjacent trials. It means each trial began and ended without other stimulations. Hence, if desired, data from successive trials could be concatenated for model estimation.

This study was approved by the ethical committee of France. Three SCI subjects participated in the experiments and signed informed consent form.

IV. IDENTIFICATION RESULTS

The proposed model structure and KF estimation method were evaluated in simulations and with experimental data.

For a model (l, m, n), we need to estimate max(l, m) + (l+m * n) dimensional meta-state. As the stimulus amplitude is quite small at the beginning of stimulation, when initializing the filter, we choose the meta-state vector ŵ0 = 0 and the output estimate ŷ0 = 0. The initial error covariance is P 0 = I, where I is an identity matrix.

A. Time-Varying Parameters Tracking in Simulation

In simulation, invariant parameters tracking was evaluated with KF firstly to investigate the stability of the muscle model estimation. Secondly, since the muscle behavior represents gradual time-varying property with muscle fatigue, we slowly change the model parameters at different time instants to simulate the changes of muscle condition. The advantage of simulation is that the true parameters are known to be compared with the estimated ones. The drawback is that it is difficult to design suitable time-varying parameters in a high order model which will not cause too much variation in output. Thereby, we chose simulation model order (2, 2, 1) to reduce the model complexity and have fewer parameters.

Fig. 3 shows the plant simulator output as well as the corresponding simulator input. The output estimate of KF is also shown in this plot. In the upper graph, the solid line indicates the actual model output, while the dotted line indicates its estimate. The corresponding parameter estimates of the model are depicted in Fig. 4. The solid lines indicate the true parameters, while the dotted lines indicate parameters estimates. All the parameters converge steadily after 5s when the parameters are static. From 35s, 50s and 60s, the model parameters were changed gradually, the estimates track the changes well. Even at 35s, a 1 suddenly decreases, the estimate tracks the variance fast. The simulation results imply that the identification method is suitable for timevarying parameters tracking in muscle contraction model.

B. Fatigue Tracking Based on Experimental Data

The proposed estimation method was further evaluated with the experimental data on three SCI subjects. The data from successive trials were concatenated for estimation. The normalized MAV and normalized torque served as input and output, respectively. The signal was processed in four steps before estimation:

• Blanking window method is used to remove stimulation artifact from EMG raw signal and extract muscle response signal (Mwave); • Lowpass filter of measured ankle torque (6th-order, cutoff frequency 100Hz) and measured EMG (6th-order, cutoff frequency 300Hz); • The filtered EMG signal is divided into epochs with each epoch containing one Mwave, the MAV is calculated every five epochs, the average torque is calculated in the same time window; • The MAV and torque are normalized with respect to their maximum value. During 100s intermittent stimulation, the torque of each subject decreased to different levels due to muscle fatigue. The nonlinear relationship between MAV of EMG and torque was revealed in experiments [START_REF] Zhang | Torque Prediction Using Stimulus Evoked EMG and its Identification for Different Muscle Fatigue States in SCI Subjects[END_REF].

The selection of model order (l, m, n) is a key stage of model parameters estimation. We chose n = 3 for the polynomial function as proposed in [START_REF] Chia | Recursive Parameter Identificaiton of Constrained Systems: An Applicaiton to Electrically Stimulated Muscle[END_REF]. The model order of linear part is determined by comparing the Rissanen's minimum description length (MDL) [START_REF] Rissanen | A Universal Prior for Integers and Estimation by Minimum Description Length[END_REF], which provides a criterion for tradeoff between model simplicity and model's fitness for the data. Finally model order l = 3, m = 4 was chosen with relatively less MDL value and simpler model.

1) Torque Prediction Performance:

To investigate the torque prediction performance in time, different prediction horizons, 6s, 18s and 30s, are tested. The idea is illustrated in Fig. 5, featuring the estimated model at time instant t0. The torque predictions are computed using (4), considering a static system within the prediction horizon. The root mean square (RMS) error and peak prediction error at every torque sample for a given prediction horizon are computed and then averaged. The prediction errors are quantified in Table I. 2) Predictive Performance with Fatigue Tacking: In order to evaluate the tracking robustness against recording missing, slight voluntary motion or failures, torque prediction is integrated to the fatigue tracking task. The event, such as torque recording error or disturbances, is assumed to result in tracking failures. The estimated model at t0 and EMG signal are used to predict generated torque.

Fig. 6 reveals the prediction performance in subject1 when muscle torque decreases to 71% by the maximum torque during 100s-stimulation. Assuming the measurement is missing or interrupted from time t0 to t1, the estimated model at t0 was used to predict the torque until instant t1. When prediction is executed, the torque estimation is only driven by EMG signal, while the online model identification is suspended. To evaluate the prediction performance in different muscle fatigue states, this process is repeated until the end of stimulation, where every 18s the measurement is suspended. We can find the prediction fits the measured torque well, suggesting that even if the torque measurement is interrupted, the predicted torque based on EMG can be used to bridge such gap for long-time application of FES.

The estimated parameters are partly depicted in Fig. 7. This plot shows the parameters denoting the weights of the past torque. It is obvious that the parameters vary gradually representing different myoelectrical mechanical correlation during fatigue. Whereas, the proposed muscle model and KF algorithm with forgetting factor provide a good prediction performance as shown in Fig. 6.

C. Discussion for Fatigue Tracking Control

The prediction result presented in Fig. 5 reveals good performance of the predictive filter, successively acquiring the muscle mechanical behavior during FES. In addition, the error values given in Table I indicate that the RMS error and peak error do not increase significantly when the prediction horizon is extended. It indicates the proposed myoelectrical mechanical muscle model fits properly the muscle behavior. If the prediction quality degrades when an expanded prediction horizon is selected, it is probably due to the difference of muscle fatigue levels. This can occur when prolonged or repetitive stimulation is delivered to the muscle. In this application, 18s prediction horizon, is considered as the optimal choice.

The myoelectrical mechanical muscle model has been described in previous work [START_REF] Erfanian | Using Evoked EMG as a Synthetic Force Sensor of Isometric Electrically Stimulated Muscle[END_REF], with recursive least squares method for identification. However, there was no significant consideration for different fatigue condition in intermittent stimulation as illustrated in Fig. 6. Moreover, the proposed approach in this paper can be further expanded to use other measurement, like joint angle, or multiple measurements update with KF thus avoiding torque sensor.

In this paper, the proposed method represents feasible and effective torque prediction only based on EMG in isometric condition for fatigue tracking. It is also promising to be applied in dynamic condition by multiplying a force-length or force-velocity function to the proposed myoelectrical mechanical muscle model as proposed in [START_REF] Farahat | A Method Identification of Electrically Stimulated Muscle[END_REF].

V. CONCLUSIONS AND FUTURE WORKS

In this work, we propose a method to predict the FESinduced ankle torque based on a time-varying myoelectrical mechanical muscle model, where the model parameters are estimated by Kalman filter with forgetting factor. The proposed prediction method is applied in fatigue tracking task with satisfactory results. The future work will be extended to adaptive FES control for muscle fatigue compensation in SCI patients. It is also possible to introduce joint dynamics model and joint angle sensors along with the proposed method, then we will work on FES control with torque prediction based on EMG without the use of torque sensor.

VI. ACKNOWLEDGMENTS

The first author is supported in part by China Scholarship Council (CSC). We would also like to warmly thank patients for their cooperation, Maria Papaiordanidou and Dr. Charles Fattal for their help to do the experiments.

Fig. 1 .

 1 Fig. 1. Muscle model structure. The present work focus on the contraction dynamics model contained in the rectangle with dashed line. The memoryless nonlinearity models muscle activation level, h(t), as a nth-order polynomial function of the instantaneous EMG signal, u(t). The linear dynamics is chosen as an ARX model.

Fig. 2 .

 2 Fig. 2. Electrodes configuration (left) and experimental set-up (right).

Fig. 3 .

 3 Fig. 3. Data set for model identification. Model input is chosen as random, Gaussian signal.

Fig. 4 .

 4 Fig. 4. Convergence and tracking of model parameter estimates for plant simulator identification. Solid lines indicate true parameters. Dotted lines indicate parameter estimates. All the initial estimates are set at zero.

Fig. 5 .

 5 Fig.5. Measured and estimated torque at time instant t0. For evaluating prediction performance, the prediction errors for prediction horizon t1, t2, and t3 (6s, 18s and 30s respectively) are computed.

Fig. 6 .

 6 Fig.6. Fatigue tracking based on EMG-torque model and Kalman filter in subject1. The model estimated at instant t0 was used to predict the torque from instant t0 to t1, while the online identification and state estimation are switched off. This process is repeated until the end of stimulation. The blue solid line indicates the measured torque, while the red dotted lines indicate estimated or predicted torque (in the range with green background color).

Fig. 7 .

 7 Fig. 7. Model parameters represent time-varying feature during different muscle states. Black, red and blue lines respectively indicate the weights of past torque a 1 , a 2 and a 3 in equation (4).

TABLE I PREDICTION

 I FILTER PERFORMANCE ON EXPERIMENTAL DATA

	Subject	Averaged error	6s	Prediction horizon 18s	30s
	S1	RMS Peak	0.0278 0.0743	0.0314 0.0962	0.0366 0.1146
	S2	RMS Peak	0.0524 0.1208	0.0534 0.1523	0.0556 0.1697
	S3	RMS Peak	0.0387 0.1036	0.0418 0.1360	0.0437 0.1510


Dkw Chu

Llm Poon
email: llmpoon@hkucc.hku.hk

Y Guan

Jsm Peiris

Novel astroviruses in insectivorous bats

Keywords: astroviruses, bats, zoonosis, coronaviruses, wild-life, evolution

Bats are increasingly recognized to harbor a wide range of viruses and in most instances these viruses appear to establish long term persistence in these animals.

They are the reservoir of a number of human zoonotic diseases including Nipah, Ebola and SARS. We report the identification of novel groups of astroviruses in apparently healthy insectivorous bats found in Hong Kong, in particular bats belonging to the genera Miniopterus and Myotis. Astroviruses are important causes of diarrhea in many animal species including humans. Many of the bat astroviruses (BatAstV) form distinct phylogenetic clusters in the genus mammastrovirus within the family Astroviridae. Virus detection rates of 36%-100% and 50%-70% were found in Miniopterus magnater and Miniopterus pusillus bats respectively, captured within a single bat-habitat during four consecutive visits spanning one year. There was high genetic diversity of viruses in bats found within this single habitat. Some BatAstV may be phylogenetically related to human astroviruses and further studies in a wider range of bat species in different geographic locations are warranted. These findings are likely to provide new insights into the ecology and evolution of astroviruses and reinforce the role of bats as a reservoir of viruses with potential to pose a zoonotic threat to human health.

Introduction

The Astroviridae are a family of non-enveloped, positive sense single-stranded RNA viruses of approximately 28-30nm in size with a characteristic star-like surface structure. The genomes of these viruses range in sizes from 6.4 -7.3 kb and are polyadenlyated at their 3' ends [START_REF] Mendez | Astroviruses[END_REF]. The viral RNA is synthesized by its RNA dependent RNA polymerase (RdRp). The genome consists of 3 open reading frames (ORFs) designated as ORF1a, ORF1b and ORF2 (Figure 1). ORF1a encodes nonstructural polyprotein 1a while the entire ORF1 encodes polyprotein 1ab with a ribosomal frame shift at the ORF1a / 1b junction [START_REF] Jiang | RNA sequence of astrovirus: distinctive genomic organization and a putative retrovirus-like ribosomal frameshifting signal that directs the viral replicase synthesis[END_REF]. Efficiency of ORF1b translation is estimated as only 25 -28% of ORF1a [START_REF] Lewis | Astrovirus ribosomal frameshifting in an infection-transfection transient expression system[END_REF]. ORF2 encodes the viral structural polyprotein which is required for virion formation. Studies with human astroviruses showed that the structural polyprotein is intracellularly cleaved into functional units including a capsid forming unit and a host binding motif unit [START_REF] Geigenmuller | Studies on intracellular processing of the capsid protein of human astrovirus serotype 1 in infected cells[END_REF][START_REF] Krishna | Identification of structural domains involved in A C C E P T E D astrovirus capsid biology[END_REF][START_REF] Mendez | Proteolytic processing of a serotype 8 human astrovirus ORF2 polyprotein[END_REF].

Astroviruses have been identified from a variety of mammals (genus mamastrovirus) and birds (genus avastrovirus) including humans (HAstV), bovine (BAstV), pigs (PAstV), ovine (OAstV), mink (MAstV), dogs, cats, mice, chickens (CAstV) and turkeys (TAstV). In most species, these viruses are associated with gastroenteritis but some avian astroviruses have been associated with both intestinal and extraintestinal manifestations [reviewed in [START_REF] Mendez | Astroviruses[END_REF]]. Human astroviruses appear to cause milder disease than rotaviruses but are the second or third most common viral agent found in children with diarrhea [START_REF] Dalton | Astrovirus acute gastroenteritis among children in Madrid, Spain[END_REF][START_REF] Herrmann | Astroviruses as a cause of gastroenteritis in children[END_REF]. Astrovirus can also cause significant disease in the elderly [START_REF] Lewis | Outbreaks of astrovirus type 1 and rotavirus gastroenteritis in a geriatric in-patient population[END_REF] and in immunocompromised patients [START_REF] Gallimore | Use of a heminested reverse transcriptase PCR assay for detection of astrovirus in environmental swabs from an outbreak of gastroenteritis in a pediatric primary immunodeficiency unit[END_REF].

Most of the surveillance studies of astroviruses focused on humans and domesticated

A C C E P T E D

animals and relatively little is known about the prevalence of astroviruses in wildlife.

The role of bats as the reservoirs for zoonotic diseases including rabies, Hendra, Nipah and Ebola has been highlighted in recent years [reviewed in ref. [START_REF] Calisher | Bats: important reservoir hosts of emerging viruses[END_REF]].

Insectivorous bats have also been shown to harbor a range of novel coronaviruses including the precursor of SARS coronavirus [START_REF] Lau | Severe acute respiratory syndrome coronavirus-like virus in Chinese horseshoe bats[END_REF][START_REF] Poon | Identification of a novel coronavirus in bats[END_REF]. Some species of bats live in close proximity to human habitation and thus it is important to have a better understanding of the virus ecology found in bats. The range and diversity of coronaviruses found in bats has led to the hypothesis that bat coronaviruses may be the precursors of most other mammalian group 1 and 2 coronaviruses [START_REF] Vijaykrishna | Evolutionary insights into the ecology of coronaviruses[END_REF]. These findings highlight the importance of identifying novel viruses in wildlife in general and bats in particular. We used random primers to detect novel viruses in bat faecal specimens. Here, we report the discovery of novel astroviruses in bats in Hong Kong.

The remarkably high prevalence and genetic diversity of astroviruses in various bat species found within a relatively small geographic area highlights the need for study in other species of bats and in other geographic locations.

A C C E P T E D

Methods

Sample collection

The sampling of bats was carried out in two phases and the sampling methods have been described previously [START_REF] Chu | Coronaviruses in bent-winged bats (Miniopterus spp.)[END_REF][START_REF] Poon | Identification of a novel coronavirus in bats[END_REF]. Phase cDNA was generated from RNA using Superscript III reverse transcriptase (Invitrogen) in a 20 µl reaction containing 150 ng of random hexamers or 0.5 µM of a gene specific reverse primer 5'-TTTGGTCCNCCNCTCCAAA-3' targeting the 3' end

A C C E P T E D

of ORF1b, 10 mM of dithiothreitol, 0.5 mM of deoxynucleoside triphosphate mix, 1x First-Strand buffer (Invitrogen), and 200 U of reverse transcriptase. Reaction mixtures were incubated at 25°C for 5 min, followed by 50°C for 60 min and then the enzyme was inactivated by heating at 70°C for 15 min. Random hexamer-generated cDNA was screened for the presence of astrovirus using hemi-nested PCR targeting RdRp gene. A 50 µl PCR reaction was set up containing 1 U Accuprime Taq DNA polymerase in 1x reaction buffer (Invitrogen), 2 µM each of forward and reverse primers and 2 µl of cDNA or 1 µl of first PCR product as a template. First-round PCR was carried out with a mix of two forward primers 5'-GARTTYGATTGGRCKCGKTAYGA-3' and 5'-GARTTYGATTGGRCKAGGTAYGA-3', and reverse primer 5'-GGYTTKACCCACATNCCRAA-3'. After an initial incubation at 94°C for 1 min, 30 cycles of amplification were carried out consisting of denaturation at 94°C for 30 sec, annealing at 50°C for 30 sec and extension at 68°C for 30 sec. Hemi-nested PCR was carried out with a mix of two forward primers 5'-CGKTAYGATGGKACKATHCC-3' and 5'-AGGTAYGATGGKACKATHCC-3' and the same reverse primers used in the first-round PCR using the same thermocycling conditions as for the first-round PCR except that 40 cycles of amplification were performed. Water controls were included in each run of the RT-PCR assay. PCR products were analyzed by standard agarose gel electrophoresis.

The expected product size of the 2nd PCR is 422 base pairs. All positive results were verified by direct DNA sequencing of the PCR amplicons.

Cloning and sequencing of PCR products

A C C E P T E D

PCR products were purified by QIAquick PCR purification kit (QIAGEN) following the manufacturer's instruction. Long PCR products (product sizes > 1000 bp) were gel purified with QIAquick gel extraction kit (QIAGEN) and then cloned into pCR2.1-TOPO plasmids (Invitrogen) for DNA sequencing. Multiple clones of a PCR product were picked and sequenced by using BigDye Terminator v3.1 Cycle Sequencing kit (Applied Biosystems). Sequencing products were analyzed by PRISM 3700 DNA analyzer (Applied Biosystems).

One astrovirus positive bat specimen, bat astrovirus AFCD337, was chosen for sequencing of the viral genome directly from the original clinical specimen by using cDNA generated by random hexamers and consensus primer RT-PCR. ORF2 region and the 3' end of the virus sample were amplified using 3' Rapid Amplification of cDNA Ends (RACE) System (Invitrogen) following the protocol provided by the manufacturer (primers and PCR conditions available on request). ORF1a (partial) and ORF1b sequences were assembled from multiple overlapping sequences derived from PCR amplicons. Nine additional ORF1a (partial) / ORF1b astrovirus sequences from other samples were obtained using similar method. The deduced sequence of these samples has at least 3-fold sequence coverage.

Phylogenetic analysis

Sequence editing and sequence identities calculation were done using BioEdit version 7.0.4 (http://www.mbio.ncsu.edu/BioEdit/bioedit.html). Alignments of nucleotide sequences and amino acid sequences were done using Clustal W (29) with default parameters. Phylogenetic trees were constructed using Clustal X V2.0 [START_REF] Larkin | Clustal W and Clustal X version 2.0[END_REF] and Mega 4 (28) by the neighbour joining method with the nucleotide substitution model of

A C C E P T E D

maximum composite likelihood and default parameters. Bootstrap values of the phylogenetic tree constructed were generated by doing 1000 replicates.

Nucleotide sequence accession numbers.

The sequence of bat astrovirus AFCD337 reported in this paper were deposited in GenBank under accession numbers xxxxx, and RdRp and the partial ORF1b sequences of other strains of bat astroviruses were deposited under accession numbers xxxxx-xxxxx. For the genetic analysis in this paper, other astrovirus genomes were retrieved from GenBank including human astrovirus type 1 strain Oxford (L23513), human astrovirus type 1 strain Dresden (AY720892), human astrovirus type 4 strain Dresden (AY720891), human astrovirus type 4 isolate Goiania/GO/12/95/Brazil (DQ070852), human astrovirus type 4 isolate Guangzhou (DQ344027), human astrovirus type 5 isolate Goiania/GO/12/94/Brazil (DQ028633), human astrovirus type 8 (AF260508), mink astrovirus (NC_004579), ovine astrovirus (NC_002469) and turkey astrovirus type-1 (NC_002470).

A C C E P T E D

Results

Detection of astroviruses in bats

A total of 262 bats were captured and rectal and throat swabs sampled in the two phases sample collection in Hong Kong in 2004 1). On the other hand, only 19 (8%) positive throat swabs were found in 246 available samples tested. With the exception of two throat swabs collected from Miniopterus pusillus and Rhinolophus rouxi in 2004, all the throat swabs positive results came from bats with astrovirus detected in rectal samples. All the positive PCR products were sequenced to confirm the identity of the amplicon and since the sequences of the different PCR amplicons were largely non-identical, PCR cross-contamination can be excluded.

Astroviruses were detected in 7 out of the 9 species of bats screened, i.e. Miniopterus magnater, Miniopterus pusillus, Miniopterus schreibersii, Myotis chinensis, Myotis ricketti, Pipistrellus abramus and Rhinolophus rouxi (Table 1). The detection rates of astroviruses in these species of bats were remarkably high and ranged from 25% to 100%. However, astrovirus was not detectable in our samples collected from Cynopterus sphinx (n=11) and Hipposideros armiger (n=10). These same specimens had previously been tested for bat coronaviruses [START_REF] Chu | Genomic characterizations of bat coronaviruses (1A, 1B and HKU8) and evidence for co-infections in Miniopterus bats[END_REF][START_REF] Chu | Coronaviruses in bent-winged bats (Miniopterus spp.)[END_REF][START_REF] Poon | Identification of a novel coronavirus in bats[END_REF].

A C C E P T E D

While 6% of bats were co-infected with both a bat astrovirus and a bat coronavirus, such co-infection appears to be randomly distributed and there was no positive or negative statistical association between the presence of these two viruses (Chi square test with Yates correction, p=0.82).

Longitudinal study of astrovirus in bats at a single habitat

The bats listed in Table 1 include (ORF1a) and the complete ORF1b, ORF2 and 3' untranslated region followed by the poly-A tail at the 3' end (Figure 1). Amino acid sequences encoded by ORFs of the novel bat astrovirus were compared with that of other astrovirus genomes including human astrovirus (HAstV) types 1, 4, 5 and 8 [START_REF] Lewis | Analysis of astrovirus serotype 1 RNA, identification of the viral RNA-dependent RNA polymerase motif, and expression of a viral structural protein[END_REF][START_REF] Silva | Molecular characterization of human astroviruses isolated in Brazil, including the complete sequences of astrovirus genotypes 4 and 5[END_REF] as well as mink astrovirus (MAstV) [START_REF] Mittelholzer | Molecular characterization of a novel astrovirus associated with disease in mink[END_REF], ovine astrovirus (OAstV) and turkey astrovirus type-1 (TAstV-1) (11) (Table 2). The findings show that the identified bat astrovirus BatAstV AFCD337 is a novel mamastroviruses virus clearly distinct from other known astroviruses. It has <53% and <27% genetic similarity to other known astroviruses in the ORF1b and ORF2 regions, respectively.

The putative ORF1a (partial) and ORF1b of the bat astrovirus sequenced have sizes of 909 nt and 1572 nt respectively. A region at 5' end of ORF1a remains unsequenced so far. ORF1a and ORF1b encodes for non-structural proteins which are essential for virus replication. Characteristic features of BatAstV AFCD337 ORF1a and ORF1b includes the protease motif in ORF1a; an astrovirus "slippery sequence" (AAAAAC) at the junction between ORF1a and ORF1b which is required for inducing a ribosomal shift event; a RdRp motif in ORF1b; and a conserved sequence at the end of ORF1b of astroviruses [START_REF] Jonassen | Complete genomic sequences of astroviruses from sheep and turkey: comparison with related viruses[END_REF]. A conserved stem-loop structure that is predicted at the 3' end of the genomic RNA of human, ovine, porcine and turkey astroviruses type-1 was not found in BatAstV AFCD337 [START_REF] Jonassen | A common RNA motif in the 3' end of the genomes of astroviruses, avian infectious bronchitis virus and an equine rhinovirus[END_REF].

The putative ORF2 of the virus has a size of 2553 nt, which is the largest astrovirus capsid gene known. The N-terminal half of the ORF2 protein, which was previously

A C C E P T E D

shown to be more conserved among astroviruses and proposed to be the core assembly domain of the viral capsid [START_REF] Krishna | Identification of structural domains involved in A C C E P T E D astrovirus capsid biology[END_REF], was also found to be relatively conserved in this bat astrovirus. The amino acid sequence similarities of this N-terminal half of the bat astrovirus capsid protein to HAstV-1 Oxford, OAstV and MAstV are 36.3%, 45.0% and 39.5% respectively; compared with <27% similarity for the ORF2 region overall (Table 2). Thus the C-terminal half of this bat astrovirus protein was highly divergent when compared with other astroviruses. This observation supports the speculation that the C-terminal half of the protein is located on the surface of the viral particle and constitutes a region of the capsid that contributes to the species-specific tropism of the virus [START_REF] Krishna | Identification of structural domains involved in A C C E P T E D astrovirus capsid biology[END_REF].

Phylogenetic analysis on astrovirus in bats.

Phylogenetic analyses on ORF1a (partial) and ORF1b region and the ORF2 confirms that BatAstV AFCD337 is a novel distinct astrovirus (Figure 2). To better define the genetic diversity within the bat astroviruses by species, time and geographic location, 77 PCR RdRp amplicons (422 nucleotides) obtained from the rectal swab samples in the screening PCR assay were selected for genetic sequence analysis. Bat astrovirus gene sequences were aligned with that of other astroviruses including MAstV, OAstV, HAstV-1 Oxford, HAstV-1 Dresden, HAstV-2, HAstV-4 Dresden, HAstV-4 Goiania, HAstV-4 Guangzhou, HAstV-5 Goiania and HAstV-8. An avian astrovirus, Turkey astrovirus was included as an out-group. A phylogenetic tree was constructed from the sequence alignment (Figure 3). The 72 astroviruses detected in bats cluster together to form a novel group of viruses within the cluster of mamastroviruses (Figure 3). Within this group are found two subgroups of viruses, one which include the majority of astroviruses detected from Miniopterus magnater, Miniopterus pusillus and

A C C E P T E D

Miniopterus schreibersii (including BatAstV AFCD337) and another subgroup which includes most of the viruses detected in Myotis chinensis and Myotis ricketti. Other than these two major groups of bat astroviruses, a few astroviruses detected in Miniopterus magnater and Miniopterus pusillus and a virus detected in a Pipistrellus abramus appear to have an outgroup relationship to the others, albeit with weak levels of statistical confidence in this phylogenetic topology (Figure 3). To further investigate this, a 750 nt region of the RdRp and ORF1b (3'end) was sequenced of these and other representative bat astroviruses. A phylogenetic tree based on the aligned protein encoding sequences was shown (Figure 4). The phylogeny of the major group of viruses related to BatAstV AFCD337 is confirmed. Interestingly however, AFCD11 from Pipistrellus abramus and AFCD57 from Miniopteus magnater appear phylogenetically related to the human astroviruses although with modest bootstrap support. Their close relationship has been also confirmed by phylogenetic analysis using an alternative method, MrBayes (data not shown).

Occasionally the same virus strain can be found in multiple bats sampled at the same habitat in a single sampling trip (i.e. strains AFCD74 and AFCD79; strains AFCD175 and AFCD228). However, most viruses detected even at the same sampling occasion at a single habitat are genetically diverse and no dominant strain could be discerned, even though the detection rates of the virus was remarkably high.

A C C E P T E D

Discussion

We report the discovery of novel astroviruses in 7 out of 9 species of apparently healthy bats captured in Hong Kong. These astrovirus positive species include Miniopterus magnater, Miniopterus pusillus, Miniopterus schreibersii, Myotis chinensis, Myotis ricketti, Pipistrellus abramus and Rhinolophus rouxi. Attempts at viral culture were so far unsuccessful (unpublished data). Phylogenetic analysis revealed that 72 out of 77 BatAstV RdRp genes analyzed clustered together to form a novel group of astroviruses. This virus group can be divided into two sub-groups, subgroup A detected from Miniopterus species and subgroup B detected from Myotis species. In the longitudinal study carried out in the abandoned mine cave habitat, the subgroup A viruses appears to be circulating between Miniopterus magnater and Miniopterus pusillus bats without any evidence of species restriction. This is in marked contrast to the bat coronaviruses within the same habitat where Bat CoV 1A and 1B appear to have a marked host restriction to Miniopterus magnater and Miniopterus pusillus, respectively (4). Multiple clones of partial ORF1 sequence (approximately 1000 nt) for 10 representative samples were analyzed and no evidence of multiple infection was found although more systematic studies on this aspect are needed.

The diversity of astroviruses in bats is remarkable. There is no significant phylogenetic clustering of viruses found within a single sampling occasion. The values of astrovirus RdRp amino acid pair-wise similarity found within a single bat species, i.e. Miniopterus magnater captured in a single habitat (mine cave) ranged between 51.1% to 100%, with 97.9% of these values lower than 90%. In contrast, the pair-wise amino acid sequence similarities between the same gene region of human

A C C E P T E D

astroviruses from geographically diverse regions were estimated to range between 92.6% to 99.1%. It has been previously reported that the amino acid sequence identities of RdRp gene (covering 80% of the RdRp gene regions analysed in this report) between 4 groups of avastrovirus, i.e. turkey astrovirus type 1-like viruses, turkey astroviruses type 2-like viruses, avian nephritis virus-like viruses and chicken-origin astroviruses detected in different regions were also highly diverse ranging from 50.1% and 73.8% [START_REF] Pantin-Jackwood | Molecular characterization and typing of chicken and turkey astroviruses circulating in the United States: implications for diagnostics[END_REF]. The high virus detection rates of BatAstV in our surveillance taken together with the marked genetic diversity of viruses from bats within the same habitat are reminiscent of our observations with bat coronaviruses (4).

These findings may suggest that bats are persistently infected with astroviruses although mark-recapture studies are needed to confirm this contention. Other mammalian astrovirus infections tend to be short-lived in immunocompetent humans or other animals. However, type 3 HAstV has been associated with persistent gastroenteritis in immunocompetent children although the same virus serotype was not repeatedly demonstrated over the full period of clinical diarrhea [START_REF] Caballero | Persistent gastroenteritis in children infected with astrovirus: association with serotype-3 strains[END_REF].

Five BatAstV sequences from Miniopterus magnater, Miniopterus pusillus and Pipistrellus abramus failed to cluster with the subgroup A and B referred to above and some of these (AFCD11 from Pipistrellus abramus and AFCD57 from Miniopterus magnater) cluster rather with the human astroviruses although with weak statistical support (Figure 4). Whether these bat viruses are related to the precursor of HAstV is yet to be further investigated. Further sequence data from these strains may help elucidate this phylogenetic association.

Evidence of recombination between astroviruses and also that between coronaviruses

A C C E P T E D

is well documented [START_REF] Lai | The molecular biology of coronaviruses[END_REF][START_REF] Pantin-Jackwood | Phylogenetic analysis of Turkey astroviruses reveals evidence of A C C E P T E D recombination[END_REF][START_REF] Strain | Genomic Analysis of Closely Related Astroviruses[END_REF]. It has been reported that a stem-loop motif in 3' untranslated region (UTR) was found conserved in mamastroviruses, TAstV-1, and in avian infectious bronchitis virus which is a group 3 coronavirus [START_REF] Jonassen | A common RNA motif in the 3' end of the genomes of astroviruses, avian infectious bronchitis virus and an equine rhinovirus[END_REF]. However, this stem loop motif was not recognized in the BatAstV AFCD337. This 3' UTR stem-loop motif is also absent in TAstV-2. However a phylogenetic analysis of 3' UTRs did not indicate a close phylogenetic relationship between the two sequences of BatAstV AFCD337 (81 nt) and TAstV-2 (accession no.: NC_005790) (196 nt) which lack the 3'UTR stem-loop structure (data not shown). Recently a report on the identification of a novel coronavirus from liver tissue of a whale with pulmonary disease and terminal acute liver failure showed that the ORF6 of the novel coronavirus possessed significant amino acid similarity to human astrovirus capsid proteins [START_REF] Mihindukulasuriya | Identification of a Novel Coronavirus from a Beluga Whale Using a Pan-viral Microarray[END_REF]. The high rates of infection of bats in the same mine cave habitat with coronaviruses and astroviruses implies frequent co-infection with both viruses.

Therefore we searched for sequence of BatAstV AFCD337 similar to bat coronaviruses using BLAST program with algorithms allowing a word-size down to seven bases. However no sequences with significant similarity was detected between bat astroviruses and coronaviruses co-circulating within the same species within the same habitat.

The discovery of novel diverse astroviruses in bats and the genetic analysis of such viruses is likely to provide new insights into the ecology and evolution of astroviruses and reinforces the role of bats as a reservoir of viruses which sometimes pose a zoonotic threat to human health. More extensive surveillance for astroviruses in bats of different species and in different geographic areas is needed to further address these questions.

A C C E P T E D

A C C E P T E D

(HAstV), mink (MAstV), ovine (OAstV) and turkey type-1 (TAstV-1). These sequences were aligned based on the encoded amino acid sequences and reverse translated back to nucleotides for the phylogenetic analysis.

A C C E P T E D

A C C E P T E D

TAstV-1

  1 was carried out in 2004 and 2005 and during this phase, swab samples were collected from different species of bats captured in the wild in different habitats in Hong Kong. Phase 2 was carried out in 2005 and 2006 and swab samples were collected from bats of the genus Miniopterus captured on four sampling occasions (June, August and December 2005 and March 2006) in a single habitat, an abandoned mine cave in Lin Ma Hang, Hong Kong, near the border with Mainland China (4). In both phases, species of bats captured for sampling were healthy and identified by a bat taxonomist. Rectal and throat swabs, together with fresh fecal samples if available, were collected. Swabs were placed in viral transport medium (Earle's balanced salt solution, 0.2% sodium bicarbonate, 0.5% bovine serum albumin, 200 µg of vancomycin per liter, 18 µg of amikacin per liter, 160 U of nystatin per liter) in screw cap tubes and transported in a cool box to the laboratory for processing. Bats were released at the site after sampling. Viral nucleic acid extraction and RT-PCR RNA from 140 µl of sample in transport medium was extracted by QIAamp virus RNA mini kit (QIAGEN) following the protocol provided by the manufacturer. Purified RNA was eluted in 60 µl of Elution buffer provided in the extraction kit.

  157 Miniopterus magnater and Miniopterus pusillus bats which were captured at four separate visits over a two year period at one habitat, an abandoned mine cave in Hong Kong. M. magnater was found throughout the period while M. pusillus was mainly found in two visits carried out in Dec 2005 and Mar 2006 and only 1 found in a visit in Aug 2005. Sixty-two (54%) out of 115 rectal swabs and 2 (2%) out of 116 throat swabs collected from M. magnater and 18 (55%) out of 33 rectal swabs and 5 (15%) out of 33 throat swabs collected from M. pusillus were positive for astrovirus. The overall positive rate in individual bats (either rectal or throat swab or both positive) for M. magnater ranged from 36-100% at each of the four visits which spanned the winter, spring and summer seasons and for M. pusillus was 50% and 70% in the two instances when adequate bats were sampled. Genetic and Phylogenetic analysis of a novel astrovirus from Miniopterus pusillus A bat astrovirus (BatAstV) AFCD337 detected in a rectal specimen from a Miniopterus pusillus bat collected in March 2006 was chosen as a representative virus for more extensive genome sequencing. Approximately 74% of the genome of this novel astrovirus was obtained by direct RT-PCR amplification from a rectal swab sample. The partial 5067 nucleotide (excluding the poly-A tail) genome was constructed by aligning sequences from multiple overlapping regions. Amino acid A C C E P T E D sequences deduced from the viral genome include part of the open reading frame 1a

Figure 1 .

 1 Figure 1. Schematic diagrams of Mink astrovirus genome and BatAstV AFCD337

Figure 2 .

 2 Figure 2. Phylogenetic analysis of partial ORF1a (~800 nt), ORF1b and ORF2

Figure 3 .

 3 Figure 3. Phylogenetic tree constructed with RdRp gene sequences (422 nt) amplified

Figure 4 .

 4 Figure 4. Phylogenetic tree constructed with 750 nt sequences of RdRp gene and

  Figure 3

  to 2006. Bats sampled included nine species, Cynopterus sphinx, Hipposideros armiger, Miniopterus magnater, Miniopterus pusillus, Miniopterus schreibersii, Myotis chinensis, Myotis ricketti, Pipistrellus abramus and Rhinolophus rouxi. A hundred and sixteen positives were detected from 250 available rectal samples tested, representing a positive rate of 46% (Table

Table 1 .

 1 Detection of astrovirus in bats by RT-PCR

			Rectal samples			Throat samples			Bats
		No. tested No. positive Positive %	No. tested No. positive Positive %	No. tested No. positive Positive %
	Insectivorous bats Hipposideros armiger Miniopterus magnater Miniopterus pusillus Miniopterus schreibersii Myotis chinensis Myotis ricketti Pipistrellus abramus Rhinolophus rouxi Fruit bats	10 122 73 3 9 12 2 8	0 67 31 3 3 10 1 1 A C C E P T E D 0% 10 0 0% 10 55% 123 4 3% 132 42% 71 6 8% 74 100% 2 1 50% 33% 9 1 11% 83% 12 5 42% 12 50% 1 0 0% 13% 8 2 25%	0 67 32 3 3 10 1 2	0% 51% 43% 100% 33% 83% 33% 25%
	Cynopterus sphinx	11	0	0%	10	0	0%	11	0	0%
	Total	250	116	46%	246	19	8%	262	118	45%

Table 2 .

 2 Amino acid sequence similarities between prototype BatAstV AFCD337, subgroup A (BatAstV AFCD68), subgroup B (BatAstV WCF140) and other astroviruses.

	BatAstV AFCD337

Acknowledgements

This project is supported by National Institutes of Health (NIAID contract HHSN266200700005C) and a Research Excellence Award to JSMP from The University of Hong Kong. The study was approved and supported by the Department of Agriculture, Fisheries and Conservation, Hong Kong, Special Administrative Region, People's Republic of China. We thank K. S. Cheung, C. T. Shek and C. S. M. Chan of the Department of Agriculture, Fisheries and Conservation, Hong Kong for facilitating this study.

A C C E P T E D


Robert Te Winkel

An Michel Rodríguez

On the Detection of Absolute Motion Relative to a Preferred Frame

From the classical perspective, it is well known that the use of the Lorentz-Fitzgerald contraction of bodies in the direction of motion explains the null results reported by Michelson & Morley in their famous experiment (and all successive repetitions). In this paper, using a one-way interferometer and taking into account Lorentz-Fitzgerald contraction, we show that it is possible to detect changes in the interference pattern proportional to the velocity of the interferometer relative to a preferred frame in which light's speed is 'c'.

Introduction

Before Young's experiments in the 19th century, when he made interfere two light wave fronts, many experiments had already been done with the purpose of determining various properties of matter and light itself.

Given that many experiments demonstrated the wave-like nature of light, and therefore the belief that light propagated through a physical medium, during the 19 th and 20 th century many different experiments were made to find this medium, so called "luminiferous aether".

The postulation of Maxwell Equations of Electromagnetism (1865) encouraged further the search of the 'ether', because it was needed to find the medium over which light propagated with speed 'c', predicted in the derived electromagnetic wave equation.

As it is today generally accepted, such 'ether' has never been found. All experiments have been reported null or results obtained have been much less than expected (as for example the very famous Michelson and Morley Experiment in 1887). Among many, experiments claim to have reduced the existence of the 'ether' with wonderful precision [START_REF] Ch | Laboratory Test of the Isotropy of Light Propagation at the 10-17 Level[END_REF].

Before the Special Theory of Relativity (STR) was postulated by Einstein in 1905, the scientific community had theories that explained the null results obtained in the experiments in search for the 1 tewinkelrobert@yahoo.es 2 anmichel.rodriguez@gmail.com 2 'ether'. Michelson & Morley's (M&M) results were explained by Lorentz, introducing his well known transformations, predicting a contraction of bodies in the direction of their motion (or, equivalently, perpendicular expansion). This is known as Lorentz-Fitzgerald contraction.

In 1905, Einstein postulated the Special Theory of Relativity (STR), which had great success among the community of physicists, and that also explained the results of the M&M experiment, as well as all previous null experiments. Although analogous formulations of the Lorentz-Fitzgerald contraction can be obtained as a result of Einstein's postulates, the interpretation of the results differs greatly. In Einstein's perspective, the 'ether' doesn't exist, there is no relative motion between light and observer. The observer always measures the speed of light as the universal constant 'c'. In classical physicist's perspective, the 'ether' does exist as the medium over which light propagates. Hence, for an observer moving relative to the 'ether', the speed of light relative to the observer is different from 'c'.

On the next formulations, the Contraction of Bodies in the direction of motion is considered as ad-hoc principle. It is shown that changes in the interference pattern can be obtained when an interferometer as the one proposed is rotated.

Theory

It will be shown that changes in the interference pattern can be measured when a one-way interferometer is rotated. Classical pre-relativistic concepts such as the Lorentz-Fitzgerald contraction of bodies in the direction of movement are used. Velocities are considered relative to a preferred inertial referential frame where the speed of light is the same in all directions and has a constant value of 'c'. Consider an interferometer as the one shown schematically on Figure 1. When analyzed at rest, light emitted by the source located at L reaches beam splitter located at A . One beam travels the optical path between points ABCD , while the other travels AD . Both beams recombine at D , and the interference pattern can be seen on the screen S .

When the interferometer travels at a speed v with respect to the preferred frame, while the beam of light travels from A to B , the mirror at B moves to position ' B . In this fashion, one beam travels the distances between points ' ' ' AB C D . For both beams to recombine at ' D at the same time and interference can be seen, the second beam must be travel the path ' ' A D .

Without loss of generality, for all purposes it can be considered the screen S as located at ' D .

In order to obtain the time taken by each beam to reach ' D , we have to consider Lorentz-Fitzgerald contraction of bodies in the direction of their motion. We express this factor with the parameter

1 v L A B C C' D D' B' A' L 2 ' ' B C vt ' ' B C ct S
For the first orientation, we can find the value of the time ' ' B C t that takes light to travel the distance between points ' B y ' C using Pythagoras's theorem with the relation

( ) ( ) ( ) 2 2 2 ' ' 2 ' ' B C B C vt L ct + = (1)
We obtain,

2 2 ' ' 2 2 B C L L t c c v α = = - (2)
We also have that

1 1 1 1 ' ' ' 2 2 2 2 AB C D L L cL L t t c v c v c v c α α α α + = + = = - + - (3)
From equations ( 1), ( 2), (3), and using the fact that ' '

' ' B C A D t t =
, we have that the time taken by light to travel paths ' ' '

AB C D and ' '

A D is given by: Using Figure 2, when the interferometer is rotated 90º (second orientation), we have for the travel times:

1 1 2 1 2L L t c c α α = + (4) 1 2 2 L t c α = (5) 1 1 2 t t t ∆ = - (6)
2 1 2 1 2L L t c c v α α = + - (7) 2 2 2 L t c v α = - (8) 2 2 2 1 2 t t t ∆ = - (9)
From equations (6) and (9), we find that

v A B' A' C' D' S 2 1 ( ) 0 t t t ∆ = ∆ -∆ = (10)
The result obtained in equation ( 10) is the expected one, and is a direct consequence of considering Lorentz-Fitzgerald contraction α . This result means that, at ' D , there is no phase change as a consequence of rotating the interferometer. However, it is important to notice that, although

0 t ∆ = , 2 1 1 1 t t > (11)
And that

2 1 2 2 t t > (12)
That is, the time taken by light to travel the distance between points ' ' '

AB C D (or ' ' A D ) in the first
orientation is smaller than the time taken to travel the same paths in the second orientation of the interferometer.

Given that the radius r of a sphere of light can be calculated as r ct = we have from equations ( 11) and

(12) that the distance traveled by light in the first orientation is smaller than the distance traveled in the second orientation of the interferometer.

It is important to emphasize that the results obtained in equations ( 11) and ( 12) are not obtained when using a M&M type interferometer, because light travels each arm of the interferometer in a two-way fashion, thus obtaining 2 1

1 1 t t = and 2 1 2 2 t t = .
The results of equations ( 11) and (12) enables us to explain why the interferometer is sensible to anisotropies in the speed of light when the interferometer is rotated: equation (10) says that there is no phase change in the interference pattern that forms in ' D , however equations ( 11) and (12) say that the radii of the spheres of light from the source to ' D (the observer) is different in each orientation of the interferometer.

The change S ∆ in this distance is, ( ) ( )

2 1 2 1 2 2 1 1 2 1 S c t t c t t cL c v c α α   ∆ = -= -= -   -   (13)
For the proposed interferometer, it can be verified that S ∆ is identically equal to zero just in the case when 0 v = . For a M&M interferometer 0 S ∆ = for all v , because 2 1 2 2

t t = and 2 1 1 1

t t = : the radii of the spheres of light from the source to ' D remain constant when the interferometer is rotated.

As it was already said, for both beams to arrive simultaneously at ' D , we have to consider that one beams traverses A , and the other is reflected at '

A when the beam splitter located at A has moved to ' A .

We can see that It is a known fact the phase i Φ in a point i P on ' i D can be calculated as ( )

1 2 i i i k d d Φ = - (14)
Where

2 k π λ =
is the wave number, and λ is the wavelength of the light used.

Because in the second orientation the beams have to travel an additional distance S ∆ to reach the observer, it should be clear that 2 1 j j d d > .

Transmitted at A (beam 1)

Reflected at A' (beam 2)

1 1 r 1 2 r 2 1 r 2 2 r First Orientation Second Orientation ∆S 0 0 2 1 d 2 2 d 1 1 d 1 2 d 1 ' D 2 ' D 1 P 2 P Reflected at A' (beam 2) Transmitted at A (beam 1)
We define

2 1 ∆Φ = Φ -Φ (15)
As can be verified, in general, 0 ∆Φ ≠ . This states the fact that there is a phase change from the point of view of the observer. This doesn't contradict the result obtained in equation ( 10), which states that the interference pattern is the same in both orientations. a a a = -. The distances j a correspond to the modulus of the projection of i i d over a .

1 1 r 1 2 r 2 1 r 2 2 r First Orientation Second Orientation ∆S 0 0 2 1 d 2 2 d 1 1 d 1 2 d 1 P 2 P a 1 a 1 a 2 a 2 a a a 1 a 1 a 2 a 2 a a Reflected at A' ( beam 2)
Reflected at A' (

Transmitted at A (beam 1) Transmitted at A (beam 1)

1 ' D 2 ' D 10
The points i P correspond to the fixed position of the observer with respect to the origin 0 .

The distances i j d of the optical paths ' ' ' AB C D and ' '

A D are equal to

2 2 i i j j j d r a = +
, where

i i j j r ct = .
The phase change that an observer measures when the interferometer is rotated (and thus passing from

1 P to 2 P ) is given by ( ) ( ) 2 2 1 1 1 2 1 2 k d d d d ∆Φ = - - - (16)
Where λ is the wavelength of light. This is equivalent to a number of fringes 2 π ∆Φ = .

From equation ( 16) it can be seen that the more unaligned the interferometer, the greater fringe shift must be measured.

The effect described in this paper can be visualized in the next figure: In the case of an M&M type interferometer, in Figure 5 the distance between 1 P and 2 P is equal to zero.

Additional Results

Because the main objective of this paper is to show theoretically the results that should be obtained when a one-way interferometer is rotated if a preferred frame of reference exists, it escapes the scope of this paper to describe in detail the experimental procedure already done by the authors to verify the results. To illustrate the reader, and to encourage further and more precise verifications of our results, it is convenient to show some data plots of the data (raw) obtained with an interferometer as the one preciously described.

because of the tridimensional nature of earths velocity, which also affect the optical paths, and are not studied in the above equations.

Conclusions

It has been described and analyzed an interferometer sensitive to anisotropies of the speed of light. The most important difference of this interferometer with a M&M type interferometer is that light travels the distances between ' ' A D and ' '

B C in a one way fashion. When the interferometer is rotated (Figure 1 to Figure 2), the radii of spheres of light between source and observer increases, which explains the changes observed in the interference pattern. In a M&M type interferometer, each arm is traveled in a two way fashion: the radii of the spheres of light between source and observer are of the same magnitude on all orientations, thus no change in the interference pattern is seen.

In this paper we have considered a rotation of the interferometer from the first to a second orientation (0º to 90º). However, it can be verified that the maximum fringe shift occurs during a rotation from the second orientation to the "fourth" (90º to 270º).

It is important to notice that Lorentz-Fitzgerald contraction was introduced, from the classical perspective, to account for the null result reported by M&M in their famous experiment. However, taking into account this effect (and thus maintaining the null result reported by M&M), and from the classical perspective (the existing of a preferred frame relative to which light's speed is isotropical and equal to the constant value of 'c'), it has been shown that using the proposed interferometer it is possible to detect changes in the interference pattern upon rotation, when 0 v > .

We can also notice that because there is no relative motion between parts of the interferometer, according to the Special Theory of Relativity (STR) no changes in the interference pattern should be expected when the interferometer is rotated. Using classical (pre-relativistic) notions it is shown that although there is no relative motion between parts of the interferometer, motion with respect to a preferred frame can be detected as fringe shifts upon rotation of the interferometer.

It is interesting to mention that the Mach-Zehnder (MZ) interferometers, because of their one-way nature, are also sensitive to the detection of anisotropies in the speed of light. In fact, the interferometer described in this paper could be understood as a variation of the MZ type. Yet more interesting is to think a MM as a MZ 'folded back upon itself' [3], with the difference that the same beam splitter is used to divide and to recombine the beams of light. This is equivalent to bring together

Figure 1

 1 Figure 1 Interferometer by which calculations are made. Thin dotted lines represent the path light travels when the interferometer at rest (v=0), thick solid lines represent the path light travels when the interferometer moves with v>0, the dashed line help visualize a right triangle of sides vt B'C' and L 2 and hypotenuse c t B'C' .

  < , where c is the speed of light, and v is the speed of the interferometer, both relative to the preferred frame. As we will analyze the interferometer in two different orientations, for all subsequent formulations, the superscript 1, 2 i = denotes the first and second orientation, respectively; the subscript 1, 2 j = denotes the paths ' ' ' AB C D and ' ' A D , respectively.

L

Figure 2

 2 Figure 2 Second orientation of the same interferometer, after a 90º rotation.

1

 1

  value in both orientations. This expression means that the interferometer has moved with velocity v during the time 1 2L c α . For an observer at ' D , the interference pattern is produced by the interference of two beams: one beam transmitted at A and another one reflected at ' A , both arriving simultaneously at ' D . The difference between the radii of the spheres of light is

Figure 3 1 P and 2 P

 312 Figure 3 In each orientation, the interference pattern is formed by two beams: one that is transmitted at A, and another that reflects at A', arriving both beams simultaneously at D'. The difference between the radii of the spheres is greater than zero and of constant value in both orientations. In the second orientation, both beams have to travel an additional distance ∆ ∆ ∆ ∆S to arrive to D'. In reference to Figure 3, 1 P and 2 P represent a fixed distance from the origin 0 . They are the position of the observer with respect to the source of light for each orientation. The distances i i j j r ct = correspond to the perpendicular distances from the source to the beam splitter, for each orientation. The i j d represent the distance traveled by light from the source to the observer in each orientation.

Figure 4

 4 Figure 4 Schematic behavior of the interference pattern formed by the two beams. The points P1 y P2 represent the position of the observer in each orientation. In this figure, we consider also a misalignment, represented by the lateral (vectorial) distance a=a 1 -a 2 between the beams.

Figure 1 ,

 1 Figure1, Figure2and Figure3represent a perfectly aligned interferometer. In practice, in order to observe interference fringes, the interferometer must be slightly unaligned. That is, the angles of the mirrors and beam splitter don't have exact 45º inclinations. Figure4show an interferometer slightly unaligned. This misalignment is represented in the figure with the lateral distance between beams

Figure 5 P 2 ,

 52 Figure 5 Unaligned interferometric setup. This figure represents a snapshot of the spheres of light propagating away from the source. The two beams that form the interference pattern are one transmitted at A and the other reflected at A'. When the interferometer is rotated the observer measures the pattern at points P 1 and P 2 , thus measuring an apparent fringe shift. In this figure, the observer measures N≅ ≅ ≅ ≅1,5.

Figure 5 1 P and 2 P

 512 Figure 5 represents a snapshot of the spheres of light propagating away from the source. The two beams that form the interference pattern are one transmitted at A and the other reflected at A'. When the interferometer is rotated, the observer measures the pattern at points 1 P and 2 P , thus measuring an

Acknowledgements

We wish to thank very specially Dr. Eduardo Greaves (Nuclear Laboratory Physics at Universidad Simón Bolívar) for all the help, support, guidance and fruitful discussions since much before the beginning of this work. Also Dr. José Ruiz-Camacho (USB) for his support. Friends and colleagues of the USB for all the help provided during research and all the endless discussions. Also we thank Paul Falstad[4] for letting us use his wonderful software to render Figure 5. Finally, we thank the Universidad Simón Bolívar for gently providing the space and tools needed to work.

Instead of rotating the interferometer, it was left at rest in the laboratory, while the Earths rotates and orbits around the sun. Experiments have been going on for approximately four months, in different configurations, and the results are consistent among themselves and the theory. In Figure 6 the horizontal axis represents time (30s precision), while the vertical axis represents the vertical position (measured in pixels) of the interference pattern's maxima. The visible undulations correspond to the movement of the interference maxima from one side to the other on the viewing screen while the earth rotates. The vertical lines tagged with "East" and "West" denote the timestamps when the stellar object HIP54255 appears in the horizon on the earth's sky, as seen from the laboratory (located N 10º32', W 66º55'). Qualitatively, it can be seen an apparent correlation between fringe displacement and the position in the sky of HIP54255. The apparent position of the object in earth's sky of this constellation is approximately equal to the direction of the velocity of the local group of galaxies with respect to CMB's rest frame [START_REF] Kogut | Dipole Anisotropy in the COBE DMR First-Year Sky Maps[END_REF].

Given that the temperature was monitored (0,1K precision) so as the input voltage of the laser (0,1V precision) and no correlation was found with the observed fringe movement, the authors believe that the deviations from the ideal sinusoid are caused by natural laser (He-Ne, 632nm) instabilities, and the beam splitters of a MZ (realigning the mirrors correspondingly). While the beam splitters are brang together, the MZ interferometer transforms gradually into a MM interferometer, also gradually transforming its one-way nature into a two-way. Just before the beam splitters are in the exact same positions, the deformed MZ interferometer looks like a slightly unaligned MM interferometer, with vestiges of sensibility to anisotropies. As has already been told in this paper, the misalignment is proportional to the number of fringes that shift when the interferometer is rotated. It is the opinion of the authors that this is the main reason why experiments using MM type interferometers have not measured identically null results, but always report velocities with magnitudes much smaller than expected (unless, of course, the interferometer is perfectly aligned, in which case an identically null result is expected).