A controversy and the writing of a history

The discussion of "small oscillations" (1760-1860) from the standpoint of the controversy between Jordan and Kronecker (1874).

Frédéric Brechenmacher

Camille Jordan and Leopold Kronecker were having a great controversy throughout the whole year of 1874, a controversy originally caused by Jordan's ambition to reorganise the theory of bilinear forms through what he designated as the "simple" notion of "canonical form":

It is known that any bilinear polynomial P = A αβ x α y β (α = 1, 2, . . . , n; β = 1, 2, . . . , n) can be reduced to its canonical form x 1 y 1 + . . . + x m y m , by linear transformations applied to the two sets of variables x 1 , . . . , x n , y 1 , . . . , y n . We now consider the following questions: 1. To reduce a bilinear polynomial P to its canonical form by orthogonal substitutions applied to the two sets of variables x 1 , . . . , x n ; y 1 , . . . , y n . 2. To reduce a bilinear polynomial P to its canonical form by the use of the same substitution on the x's and the y's. 3. To reduce simultaneously two bilinear polynomials P and Q to a canonical form. [...] The third problem has already been solved by M. Weierstrass [...] the solutions given by the geometers from Berlin are, nevertheless, incomplete; we will therefore suggest an extremely simple new method that holds no exception [...]. We will show that the problem of the simultaneous reduction of two functions P and Q is identical to the problem of the reduction of a linear substitution to its canonical form. ([8], 7-11, translation F.B.) In Jordan's 1873 note to the Parisian academy quoted above, the first problem referred implicitly to geometry and Cauchy's results on the principal axis of conics and quadrics, the second one referred to the arithmetic of quadratic forms and the works of Gauss and Hermite, the third problem referred to analytical mechanic and the solution given by Lagrange to the systems P Y + QY = 0 of linear differential equations with constant coefficients. What was at stake in the 1874 controversy 1 was the organisation of the theory of bilinear form, a theory that was considered as giving a new "homogeneous" and "general" treatment to different problems referring to various theories developed throughout the 19th century. According to Kronecker, the main problem of the theory was the characterisation of the equivalence of pairs of forms (P, Q) (Jordan's third problem in the above quotation). Two theorems stated independently by Weierstrass and Jordan could be used to solve this problem and it was their opposition that generated the 1874 controversy. [START_REF] Brechenmacher | Histoire du théorème de Jordan de la décomposition matricielle (1870-1930[END_REF] On the one hand, Weierstrass had defined in 1868 a complete set of polynomial invariants of non singular pairs of bilinear forms that were computed from the determinant |P + sQ| and designated as the elementary divisors of (P, Q). Looking for invariants was a typical method of the theory of bilinear forms developed in the 1860's as a local field of research limited to a few geometers from Berlin such as Christoffel, Kronecker and Weierstrass. It was in the very different context of group theory that, on the other hand, Jordan had stated in 1870 that a linear substitution could be reduced to a simple canonical form. The two theorems were stated independently and belonged to distinct theories until 1873 when Jordan claimed that his theorem could be used in the theory of forms.

The connection between the two theorems was made thanks to a hundred-yearold mechanical problem. In 1766 Lagrange had devised a method for the integration of a system of n linear equations with constant coefficients for the purpose of describing the small oscillations of a swinging string loaded with an arbitrary number of weights. Lagrange's method was based on an implicit mechanical representation: considering that the oscillations of a swinging string loaded with n weights could be represented as the combination of independent oscillations of n strings loaded with a single weight, one should be able to represent the differential system as a combination of n independent equations dy i = σ k y i (i = 1, . . . , n) (σ k representing the periods of the independent oscillations). Lagrange proved that σ 1 , . . . , σ n were the roots of an algebraic equation of the n th degree computed from the differential system, his method therefore failed when multiple roots appeared. When the method was later applied to describe the "secular oscillations" of the planets on theirs orbits, the problem caused by the multiplicity of roots questioned the stability of the solar system and generated a century-long discussion around what was to be known as the "equation of secular variations." [START_REF]Regards croisés sur Camille Jordan[END_REF] The meeting point of Weierstrass's invariants and Jordan's canonical form was their capacity to give a complete resolution to the above problem, regardless of the multiplicity of roots. Weierstrass first proved in 1858 that the multiplicity of roots was irrelevant to the subject of mechanical stability. [START_REF] Cauchy | Sur l'équation à l'aide de laquelle on détermine les inégalités séculaires du mouvement des planètes[END_REF] His 1868 theorem then gave a complete characterisation of the systems of linear equations with constant coefficients. Jordan reached the same conclusions in the years 1871-1872, see [7]. They both referred to the different problems handled by Lagrange, Laplace, Cauchy, Hermite etc. as a single problem of transformation of pairs of forms: the mathematical notion of forms unified various problems that were in the past considered to belong to different theories such as celestial mechanic or geometry. The question at the heart of the controversy was then to decide which methods should be used to organise the theory of bilinear forms: reducing to canonical forms or computing invariants? On the one hand, Jordan argued that the main method of the theory should be the algebraic reduction of forms to their simplest canonical forms. On the other hand, Kronecker objected that the theory belonged to the arithmetic of forms in the tradition of Gauss and the main problem should thus be characterising equivalence classes by the computation of invariants. Kronecker argued that Weierstrass' invariants could be effectively computed because they were defined as g.c.d.s of the minors extracted from the polynomial determinant |P + sQ|4 ; on the contrary, Jordan's canonical form could not be computed effectively, in general, and should thus be considered a "formal notion" with no "objective meaning." With the aim of criticizing Jordan's canonical form, Kronecker wrote a critical history of the algebraic methods used during the discussion of small oscillations. These methods were blamed for their tendency to thinking in terms of the "general" case with little attention given to difficulties that might be caused by assigning specific values to the symbols. Weierstrass' invariants were considered by Kronecker to be exemplary of a "truly general" theorem as opposed to the "so called generality" that consisted in focusing on the "general case" where |P + sQ| = 0 had no multiple roots. Kronecker was thus the first to stress a history of what the historian T. Hawkins referred to as the "generic reasoning" in the 18 th -19 th centuries algebra ( [5], 122).

The 1874 controversy opposed two theorems that were two different ends given to a century-long mathematical discussion. The controversy therefore combined mathematical and historical arguments and gave rise to the writing of a history of the methods used by Lagrange, Laplace, Cauchy and Weierstrass. A historical study of how Kronecker and Jordan referred to the past sheds some light on some tacit ideals which led to an opposition on what should be the meanings of "generality," "simplicity," "algebra," 'arithmetic" and the notion of "form" in mathematics. Even long after the controversy, the tension between canonical forms and invariants played a major role in the development of the theory of bilinear forms and the creation of the theory of matrices. A careful study of the Jordan-Kronecker controversy was used as a preliminary to a wider historical understanding of the development of the theory of matrices in my doctoral thesis (see [1]).

Provided that (P, Q) is non singular that is |P + sQ| does not vanish identically and that the forms coefficients belong to an algebraically closed field.

For a history of this discussion of small oscillations, see[START_REF] Hawkins | Cauchy and the spectral Theory of Matrices[END_REF] and[START_REF]Regards croisés sur Camille Jordan[END_REF].

From a contemporary standpoint, the system Y = AY is stable if and if the matrix A is diagonalizable, that is if and if A's elementary divisors are simple.

From a contemporary standpoint, Kronecker defined the invariant factors of a matrix. These invariants can be defined for any matrix with coefficient belonging to a principal ring while Jordan's canonical form requires an algebraically closed field.


Maud Delattre

INFERENCE IN MIXED HIDDEN MARKOV MODELS AND APPLICATIONS TO MEDICAL STUDIES by

Keywords: 2000 Mathematics Subject Classification. -62-02. Key words and phrases.hidden Markov models, mixed-effects, longitudinal data, stochastic approximation EM, forward recursions, maximum a posteriori, Viterbi algorithm

The aim of the present paper is to document the need for adapting the definition of hidden Markov models (HMM) to population studies, which rigorous interpretation typically requires the use of mixed-effects models, as well as for corresponding learning methodologies. In this article, mixed hidden Markov models (MHMM) are introduced through a brief state of the art on hidden Markov models and related applications, especially focusing on disease related problems. Making the main assumption that a given pathology can be considered at different stages, hidden Markov models have for example already been used to study epileptic activity or migraine.

Mixed-effects hidden Markov models have been newly introduced in the statistical literature. The notion of mixed hidden Markov models is particularly relevant for modeling medical symptoms, but the data complexity generally requires specific care and the available methodology for MHMM is relatively poor. Our new approach can be briefly described as follows. First, we suggest to estimate the population parameters with the SAEM (Stochastic Approximation EM) algorithm, which has the property to converge quickly. The well-known forward recursions developed for HMM allow to compute easily the complete likelihood at each step of the MCMC procedure used within SAEM. Then, for dealing with the individuals, we suggest to estimate each set of individual parameters with the MAP (Maximum A Posteriori) of the parameter distributions. Finally, the hidden state sequences are decoded using the Viterbi algorithm. Some Monte-Carlo experiments are presented to illustrate the accuracy of our algorithms.

Résumé (Infrence dans les modèles de Markov cachés à effets mixtes et applications médicales)

Cet article veut montrer la nécessité d'étendre la définition des modèles de Markov cachés, ainsi que leurs méthodes d'estimation, au cadre des études de population, dont l'interprétation rigoureuse passe par l'utilisation de modèles à effets mixtes. Nous motivons alors les intérêts des modèles de Markov cachés à effets mixtes au travers d'un état de l'art succinct sur les modèles de Markov cachés et leurs nombreuses applications. Nous nous limiterons à des problématiques médicales. Dans ce cadre, les modèles de Markov cachés supposent que l'évolution des maladies peut s'interpréter à travers différents états. En effet, la distinction de divers stades évolutifs dans la maladie justifie l'application de modèles de Markov cachés à certaines pathologies, comme cela a déjà été le cas pour la migraine, ou encore l'épilepsie. La définition des modèles de Markov cachés à effets mixtes est très récente. Ces nouveaux modèles sont des candidats intéressants pour la modélisation de symptômes, celle-ci s'opérant généralement à partir de données complexes par leur structure, et nécessitant une démarche d'analyse particulière. Actuellement, les méthodes d'apprentissage pour les modèles de Markov cachés à effets mixtes restent peu nombreuses et ne se limitent qu'à l'estimation des paramètres de population. Or, de manière générale, les modèles mixtes s'interprètent à deux niveaux : au niveau de la population d'une part, puis au niveau des individus isolés d'autre part. Ce second point n'a pas encore été abordé dans le cadre de nos nouveaux modèles. Sur ce constat, nous proposons de compléter les approches déjà établies par la démarche d'apprentissage suivante. Pour commencer, l'algorithme SAEM sera adapté aux modèles de Markov cachés à effets mixtes pour en estimer les paramètres de population. A chacune de ses itérations, l'étape d'approximation stochastique mettra à profit la procédure forward, qui propose une méthode de calcul rapide de la vraisemblance dans les modèles de Markov cachés. Ensuite, cette première estimation à l'échelle populationnelle nous permettra d'établir les paramètres propres à chaque sujet par maximisation a posteriori de leur distribution. Enfin, la simple mise en oeuvre de l'algorithme de Viterbi dans les modèles de Markov cachés individuels nous donnera les séquences d'états les plus probables pour chaque individu. En dernier lieu, nous nous attacherons à illustrer les propriétés de nos algorithmes sur un court exemple. -Hidden Markov models first appeared in the statistical literature in the late 60's, in a series of articles from Baum and coworkers ([1], [START_REF] Baum | An inequality with applications to statisical estimation for probabilistic functions of markov processes and to a model for ecology[END_REF], [START_REF] Baum | A maximization technique occuring in the statistical analysis of probabilistic functions of markov chains[END_REF]).

Those new models weren't referred to as "hidden Markov models" yet. The expression "probabilistic functions of finite state Markov chains" was rather used, reflecting quite well the definition of hidden Markov models. HMM applications, especially to speech recognition, abounded in the 70's and were at the origin of some methodological developments for learning these new models. In 1989, Rabiner published a tutorial ( [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF]) in which hidden Markov models were rigorously defined. This paper also clarified the related methodologies and gave several application examples devoted to speech recognition problems.

In hidden Markov models, a double process is assumed, among which only one is observed. Consequently, two levels are separately considered in the definition of such models. First comes the definition of the hidden process, which is a Markov process, generally assumed to have memory one; second comes the definition of the observations' process conditional on the hidden one. In this section, we will consider a parametric framework and Ψ will denote the vector of all HMMs' parameters.

Let {Z j } j∈N be a stationary Markov chain with a discrete and finite state space S = {1, . . . , S}. In the following, let π be the stationary distribution of the chain and let A denote the S × S transition matrix of the Markov process, and ∀s, s ∈ S 2 , let also a s,s be the transition probability associated with the transition from state s to state s :

a s,s = P Ψ (Z j+1 = s |Z j = s) ; ∀j ≥ 1 where S s =1 a s,s = 1; ∀s ∈ S
In hidden Markov models, a second process {Y j } j∈N enables inference on the latent process. More precisely, the j th observation is assumed to be a probabilistic function of the j th hidden visited state.

For discrete observations, emission probabilities are introduced to specify how the observations are related to the state sequence. Let O be the observation space, and let b o,s be the probability associated to observation o ∈ O while current (unknown) state is s ∈ S:

b o,s = P Ψ (Y j = o|Z j = s) ; ∀j ∈ N
The emission probabilities are such that:

o∈O b o,s = 1; ∀s ∈ S
For example, when the observations are count data, Poisson distributions with parameters λ 1 , λ 2 , . . . , λ S are usually chosen to model emissions in states 1, 2, . . . , S respectively. The model for observations would thus read:

b o,s = P λs (Y j = o|Z j = s) = exp(-λ s ) λ o s o! ; ∀s ∈ S
At last, an additional assumption holds: conditionally to Z j , Y j is independent of Y 1 , ..., Y j-1 , Y j+1 , ... and Z 1 , ..., Z j-1 , Z j+1 , ....

From the above definitions and assumptions, the general expression of the likelihood in HMMs becomes:

(1)

L (y 1 , . . . , y J ; Ψ) = z1,...,z J ∈S J π (z 1 ) J j=1 b yj ,zj J-1 j=1 a zj ,zj+1
where J is the length of the observed sequence.

Remark: extension to continuous observations is straightforward by defining conditional emission densities, for example a Gaussian distribution with parameters (m s ; σ 2 s ) in state s ∈ S.

1.2. Learning methods. -Infering hidden Markov models is challenging, mostly due to the complex expression of the likelihood and to the non observable visited states. As a consequence, hidden Markov models are associated with three "basic problems". Quoting from [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF], (i) computing the likelihood, (ii) estimating the model parameters (emission probabilities, transition probabilities, and possibly the probability distribution of the states at time 1), and (iii) decoding the most probable state sequence for a given sequence of observations. Some algorithms are referenced and discussed in [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF], giving potential solutions to (i), (ii) & (iii). Among those presented in [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF], the forward procedure, the Baum-Welch algorithm and the Viterbi algorithm are the most relevant ones, applying respectively to (i), (ii) & (iii). The Baum-Welch algorithm is an EM-type algorithm which is expected to compute the maximum likelihood estimator (MLE). The question of the consistency of the MLE has been largely investigated. This is however a complex problem, and very strong assumptions were usually required to get the MLE consistency. Among other works, those of Leroux ([5]) and Douc and Matias ( [START_REF] Douc | Asymptotics of the maximum likelihood estimator for general hidden markov models[END_REF]) could be cited. More recently, Douc et al. demonstrated the consistency of the MLE under very weak assumptions in [START_REF] Douc | Consistency of the maximum likelihood estimator for general hidden markov models[END_REF]. Their result even holds in situations where the state space is not compact.

From hidden Markov models to mixed hidden Markov models

Many authors suggested applying hidden Markov models to deal with some biological problems. Hidden Markov models have become a very successful modeling tool in molecular biology, and applications to genetics abound in the statistical literature. Hidden Markov models have also been used in epidemiology, as an alternative to compartimental SIR (Susceptible-Infectious-Recover ) and SIS (Susceptible-Infectious-Susceptible) models, to study the propagation of epidemics or infections in specific populations ( [START_REF] Strat | Monitoring epidemiologic surveillance data using hidden markov models[END_REF], [START_REF] Cooper | The analysis of hospital infection data using hidden markov models[END_REF]). In this context, the observations consist of counts of infected patients at successive time points. As such phenomenons are characterized by patient to patient transmission, successive outcomes can't be considered as independent outcomes and assuming a Markov structure is of strong interest. Sometimes, the studied disease is asymptomatic on some patients. Hidden Markov models could thus be an adapted way to describe the infection latent process, through the estimation of important epidemiological parameters, or even through the distinction of periods with different transmission rates. For example, Le Strat and coworkers were interested in influenza-like illness (I.L.I.) and poliomyelitis in 1999 ( [START_REF] Strat | Monitoring epidemiologic surveillance data using hidden markov models[END_REF]). In [START_REF] Strat | Monitoring epidemiologic surveillance data using hidden markov models[END_REF], I.L.I. incidence rates were modeled as a two-state Gaussian hidden Markov model, leading to a clear distinction between epidemic and non epidemic periods. The number of monthly Poliomyelitis cases were analyzed in the same manner, with a two-state Poisson hidden Markov model. In a same way, Cooper and Lipsitch were interested in nosocomial infections ( [START_REF] Cooper | The analysis of hospital infection data using hidden markov models[END_REF]) and proved that models with a latent Markov process granted the best fit of their data. Then, although hidden Markov models have not been widely used in this field, some authors showed that HMMs could be particularly relevant for the study of chronic illnesses, such as migraine ( [START_REF] Anisimov | Analysis of responses in migraine modelling using hidden markov models[END_REF]), multiple sclerosis ( [START_REF] Albert | Time series for modelling counts from relapsing remitting disease : application to modelling disease activity in multiple sclerosis[END_REF], [START_REF] Altman | Application of hidden markov models to multiple sclerosis lesion count[END_REF]) or epilepsy ([13]). Here, the transitions between some unobserved states, whom total number is possibly predefined, are supposed to describe the evolution of disease symptoms, like the daily numbers of seizure counts in the context of epilepsy ( [START_REF] Albert | A two state markov mixture model for a time series of epileptic seizure counts[END_REF]), the monthly numbers of lesion counts in the context of multiple sclerosis ( [START_REF] Albert | Time series for modelling counts from relapsing remitting disease : application to modelling disease activity in multiple sclerosis[END_REF], [START_REF] Altman | Application of hidden markov models to multiple sclerosis lesion count[END_REF]), or headache scores when dealing with migraine ( [START_REF] Anisimov | Analysis of responses in migraine modelling using hidden markov models[END_REF]).

To study diseases, the interests of hidden Markov models are numerous. First, those particular models are quite easily interpretable, and appear to show up similarities with the biological process that governs the pathologies. The Markov states are thus associated with distinct stages or seriousness degrees for the studied illness, and the assumption is that patients alternate periods in those stages. As an example, patients with multiple sclerosis seem to undergo relapsing and remitting periods ( [START_REF] Albert | Time series for modelling counts from relapsing remitting disease : application to modelling disease activity in multiple sclerosis[END_REF], [START_REF] Altman | Application of hidden markov models to multiple sclerosis lesion count[END_REF]). In the same way, Albert ( [START_REF] Albert | A two state markov mixture model for a time series of epileptic seizure counts[END_REF]) assumes the epileptic patients to go through two distinct stages, namely a low and a high seizure susceptibility. In related clinical trials, markers are used to support disease diagnosis, and explaining their value by indirectly observed illness stages seems to have a biological meaning; at least this approach is widely used. Seeing this, the use of latent states, typically through mixture models, to study disease dynamics on specific patients is natural. It is also quite reasonable to assume that consecutive values of a biomarker for a given patient are interdependent. For example, past events or repeated past passages to acute forms of the disease could reinforce susceptibility to the illness. It is thus justified to enrich the mixture with a first-order memory, leading to hidden Markov models. Having a well-founded biological interpretation, we could also imagine hidden Markov models could improve the understanding of the process underlying some more obscure pathologies.

Therefore, hidden Markov models have conceptual validity in some disease studies. The estimated model parameters help thus to interpret the disease process at several levels. First, through the emission distributions, they give some idea of the way the biomarker's values are related to the hidden states. Second, the estimated transition probabilities help to see how state changes are frequent in the studied population. Interpretation could even be carried on by including covariates and regression variables in the parametrization of the model.

However, modeling disease using hidden Markov models is not straightforward. In particular, while hidden, the "design" of the underlying Markov process could be challenging. When enough knowledge on the disease of interest exists, the number of hidden states can be a priori fixed; but most often, the number of hidden illness stages is unknown. Several numbers of states have to be tried and adapted selection criteria are needed. As an example, hidden Markov models from two to six states are tried to model I.L.I. data in [START_REF] Strat | Monitoring epidemiologic surveillance data using hidden markov models[END_REF], and the BIC criteria is chosen to discriminate the most adapted model from the others, leading to a five-state Gaussian hidden Markov model. Le Strat stresses the lack of interpretability of such a result.

More specific modeling difficulties occur when modeling the outcomes of clinical studies, mainly due to their structure. Indeed, several patients are included, and are subject to repeated measurements.

As hidden Markov models are a possible way to analyze one particular sequence of data, the first approach consisted in considering as many hidden Markov models as included patients. Each individual set of parameters was therefore estimated independently of the others. Albert followed this approach to epileptics'seizure count data ( [START_REF] Albert | A two state markov mixture model for a time series of epileptic seizure counts[END_REF]) and to multiple sclerosis data ( [START_REF] Albert | Time series for modelling counts from relapsing remitting disease : application to modelling disease activity in multiple sclerosis[END_REF]). However, by continuing Albert's work on multiple sclerosis, Altman ( [START_REF] Altman | Application of hidden markov models to multiple sclerosis lesion count[END_REF]) underlined estimation inaccuracy, and noted the obtained estimates were always associated with large s.e.

Clearly, the individual fit approach to longitudinal data has the major drawback of incorrectly capturing the heterogeneity among patients. Indeed, the complete set of individual estimates only give a limited summary of the variation or heterogeneity of the individual parameters. The need for ". . . a model [that would] describe all patient's data simultaneously" was therefore argued for the first time in Altman's article ( [START_REF] Altman | Mixed hidden markov models : an extension of the hidden markov model to the longitudinal data setting[END_REF]). On the same idea as mixed models, the heterogeneity characterizing the data would be finely taken into account by including i.i.d. random effects in each patients' hidden Markov model parameters definition. It would also be a way to foresee possible correlations between parameters. This way, Altman supposed defining a hidden Markov model with random parameters would help to increase the precision of the estimates, and would best capture the potential variation among patients. Those remarks, dating back to 2005, are at the origin of mixed hidden Markov models.

3. Mixed hidden Markov models 3.1. Definition. -A rigorous definition of mixed hidden Markov models by Altman followed in 2007 in [START_REF] Altman | Mixed hidden markov models : an extension of the hidden markov model to the longitudinal data setting[END_REF]. Parallel to this work, Ip and coworkers also published an article on mixed hidden Markov models in 2007 ( [START_REF] Ip | Mixed effects hidden markov models[END_REF]). In both papers, mixed hidden Markov models appear as an extension of "classical" hidden Markov models to deal with the specific contexts using a population approach.

Mixed hidden Markov models include several levels of definition. Assume we have at our disposal data from n subjects. A hidden Markov model is used for each individual set of data, while the parameters for each individual model are assumed to be random with a common probability distribution. As for HMMs' definition, we will consider a parametric framework. Using the same notations, O is the common observation space, and S is the common state space.

3.1.1. Definition of n "distinct" hidden Markov models. -The first step of a MHMM's definition consists in specifying a hidden Markov model for the observations of each of the n subjects. More precisely, the distribution of the observations for each individual is based on a Markov chain, which sequence of visited states is unknown. Let us restrict to subject i (1 ≤ i ≤ n). Let n i be the number of observations for this subject, and let Y i = (y i1 . . . y ini )

T and Z i = (z i1 . . . z ini ) T be respectively the sequence of observations for individual i and his sequence of hidden states. As MHMMs' definition mainly goes through the specification of individual hidden Markov models, as many sets of parameters as subjects are needed instead of only one set of parameters Ψ for HMMs. Let Ψ i denote the vector of parameters for subject i. Typically, Ψ i is part of the definition of 1. the emission distributions, via a series of emission probabilities for discrete observations:

b (i) o,s = P Ψ i (y ij = o|z ij = s) ; ∀1 ≤ j ≤ n i ; ∀o ∈ O; ∀s ∈ S 2. the transition matrix a (i) s,s = P Ψ i (z i,j+1 = s |z ij = s) ; ∀j ≥ 1; ∀(s, s ) ∈ S 2
3.1.2. Model for the individual parameters. -The n vector of individual parameters Ψ i have a same probability distribution. The parameters θ of this population distribution are the so-called population parameters. We will consider a linear Gaussian model for the (transformed) individual parameters that can include covariates:

h (Ψ i ) = µ + C i β + D i η i η i i.i.d. N (0, Ω)
where h is a vector of link functions, C i and D i are known matrices of covariates for individual i, µ and β are unknown vectors of fixed effects, and Ω captures the variability of individual behaviors that the covariates can't explain themselves. Here, θ = (µ, β, Ω).

With such a hierarchical definition, a single statistical model describes the whole individuals'data simultaneously while taking into account the potential heterogeneity among patients.

The observed likelihood is given by

L (Y 1 , Y 2 , . . . , Y n ; θ) = n i=1 L (Y i , Ψ i ; θ) dΨ i (2) = n i=1 L (Y i |Ψ i ) L (Ψ i ; θ) dΨ i
where L (Y i |Ψ i ) has a similar expression as the observed likelihood of a "classical" HMM given in (1). This observed likelihood cannot be computed in a closed form and this complex expression makes the model inference directly intractable.

Inference in MHMM.

-Inference in mixed hidden Markov models is not straightforward, and only some partial methods have already been suggested. Here, we will finally put forward a different way to grasp MHMM's learning.

First steps.

-Mixed hidden Markov models are somehow "new" models. As a consequence, their usage is documented in a very limited way and the related methodologies are not well established. The maximum likelihood approach could be used for estimating the population parameters but the complex expression for the likelihood makes its maximum difficult to locate. Mixed hidden Markov models can be viewed as missing data models where the visited states and the individual parameters are the non observed data. As a consequence, the EM algorithm seems to be a natural parameter estimation method for such models ( [START_REF] Ip | Mixed effects hidden markov models[END_REF]) but the E-step cannot be performed in a closed-form. Altman suggested alternative methods, such as quasi-Newton methods or Gaussian quadrature methods, or even the MCEM algorithm ( [START_REF] Altman | Mixed hidden markov models : an extension of the hidden markov model to the longitudinal data setting[END_REF]). Nevertheless, these algorithms are time expansive. Several days may be required to estimate the model parameters when the number of random effects in the model exceeds three. This forces her to restrict her attention to models involving random effects on the emission distribution only. Considering the problem of predicting mastisis prevalence in cows, Detilleux suggested to estimate the model parameters of a mixed hidden Markov model using a Gibbs sampler ( [START_REF] Detilleux | The analysis of disease biomarker data using a mixed hidden markov model[END_REF]).

Our methodology.

-Performing parameter inference in mixed hidden Markov models has been underlined to be a complex problem. Knowledge of the population parameters is necessary to grasp the mean tendency, as well as its variability among individuals; but it is not enough when focusing on a particular individual. Using the population parameters alone could bias individual diagnosis. That's why it is important to divide the mixed hidden Markov models' problem into three main questions:

1. First naturally comes the question of estimating the population parameters.

2. Then the individual sets of parameters have to be estimated.

3. Estimating the most probable individual state sequences is the final issue to address.

In this paragraph, statistical methods dealing with the three points above are suggested. The use of the SAEM algorithm for estimating the population parameters is the most original part of our methodology, and is thus more detailed.

Population parameters' estimation.

-In models such as mixed hidden Markov models, the E-step of the EM algorithm is not directly tractable. Then, we propose to adapt the MCMC-SAEM algorithm ( [START_REF] Kuhn | Coupling a stochastic approximation version of em with an mcmc procedure[END_REF]) to the mixed hidden Markov model setting. Each iteration of the algorithm can be decomposed into three steps.

The non observed data are simulated (simulation step). These simulated data are used in a second step together with the observations to approximate the complete likelihood (stochastic approximation step). This likelihood is then be maximized to update the estimation of the parameters (maximization step).

In the context of mixed hidden Markov models, the first idea would be to consider the individual parameters (Ψ i ) and the Markov chains (Z i ) as the non observed data. Indeed, the conditional distribution of (Ψ i , Z i ) can easily be simulated by MCMC and the complete likelihood L (Y 1 , . . . , Y n , Z 1 , . . . , Z n , Ψ 1 , . . . , Ψ n ; θ) can easily be maximized.

Even if this first version of the algorithm can be implemented and gives good results, considering the Markov chain as a nuisance parameter of the model allows to propose a much more simple and efficient procedure. Note that a quick computation of the n individual likelihoods L (Y i , Ψ i ; θ) is the key of the algorithm. Indeed, the following decomposition allows many simplifications extremely useful for implementing the SAEM algorithm:

(3) L (Y i , Ψ i ; θ) = L (Y i |Ψ i ) L (Ψ i ; θ)
Computing L (Y i |Ψ i ) turns out to be easy by making use of the forward recursions that are part of the well-known Baum-Welch algorithm which allows computing the observed likelihood in hidden Markov models. Then, L (Ψ i ; θ) derive from the Gaussian distribution and is easy to compute and to maximize. Let us describe iteration k of the algorithm. Here, θ k denotes the current estimate of the population parameters.

Simulation

The k th iteration begins with drawing Ψ (k) i from the conditional distribution p (Ψ i |Y i ; θ k ) for all 1 ≤ i ≤ n. The Hasting-Metropolis algorithm used for this simulation step requires to compute L (Y i |Ψ i ; θ k ) in a closed form for evaluating each acceptance probabilities. As mentioned above, computing this conditional likelihood is straightforward thanks to the forward procedure.

Stochastic approximation

Follows a stochastic approximation of the log likelihood:

Q k (θ) = Q k-1 (θ) + γ k n i=1 log L Y i , Ψ (k)
i ; θ -Q k-1 (θ)
where (γ k ) k≥0 is decreasing to 0 over iterations. Q k (θ) can be written as the sum of two terms among which only one depends on parameter θ:

Q k (θ) = R k + T k (θ)
where

R k = R k-1 + γ k n i=1 log L Y i |Ψ (k) i -R k-1
and

T k (θ) = T k-1 (θ) + γ k n i=1 log L Ψ (k) i ; θ -T k-1 (θ)
Then, it is equivalent to maximize Q k (θ) or T k (θ) with respect to θ and our stochastic approximation step would just reduces in computing T k (θ).

Maximization

k th iteration ends in maximizing T k to update the estimation of θ:

θ k = argmax θ T k (θ)
Iterations of this procedure are repeated until numerical convergence of the sequence (θ k ) to some estimate θ is achieved.

Computing the standard errors (s.e.) of the estimated parameter θ requires computing the Fisher Information Matrix (F.I.M.). We propose to estimate the F.I.M. using the stochastic approximation procedure described in [START_REF] Kuhn | Coupling a stochastic approximation version of em with an mcmc procedure[END_REF] and based on the Louis formula.

Individual parameters' estimation.

-After estimating the population parameters with the SAEM algorithm, each individual parameter estimate Ψ i can be calculated through the MAP (Maximum A Posteriori) method: Ψi = argmax

Ψ i p Ψ i |Y i ; θ
Such maximization for each individual requires some optimization procedure.

Remark: An alternative would be to estimate the conditional mean E Ψ i |Y i ; θ with the MCMC procedure used within the SAEM algorithm. -Our simulations were inspired by the quite numerous studies on epileptic activity. Similarly to the works cited above, we assumed the existence of a hidden Markov chain, which would condition the intensity of the seizures in epileptic patients. The common intuition is the following. The first and the second states would respectively be associated with a low and a high epileptic activity. Periods in both states would thus alternate in epileptic patients. As in [START_REF] Albert | A two state markov mixture model for a time series of epileptic seizure counts[END_REF], the emission distributions are chosen to be Poisson distributions. This means that conditional to the state the number of daily seizures for a given epileptic patient is assumed to follow a Poisson distribution. Let λ [START_REF] Albert | Time series for modelling counts from relapsing remitting disease : application to modelling disease activity in multiple sclerosis[END_REF] and p (i) 21 be individual i's transition probabilities associated respectively with the transitions from state 1 to state 1 and from state 2 to state 1.

(i) 1 < λ (i) 2 . Let also p (i)
Our model is the following: logit(p

(i) 11 ) = γ 1 + η 1i (4) logit(p (i) 21 ) = γ 2 + η 2i (5) log(λ (i) 1 ) = log(λ 1 ) + η 3i (6) log(α (i) ) = log(α) + η 4i (7) λ (i) 2 = λ (i) 1 + α (i) (8)
The random effects are assumed to be independent and normally distributed:

η i = (η 1i , η 2i , η 3i , η 4i ) i.i.d.

N (0, Ω)

θ corresponds here to the concatenation of the fixed effects (γ 1 , γ 2 , λ 1 , α) and the elements of the variance-covariance matrix Ω.

4.2.

A first numerical experiment. -One dataset with 200 individuals and 100 observations per subject were simulated using the following values for the fixed effects: γ 1 = 1.4, γ 2 = -1.4, λ 1 = 0.8, α = 2.3. The random effects were simulated assuming a diagonal variance-covariance matrix Ω with the following diagonal elements: ω 2 γ1 = 0.1, ω 2 γ2 = 0.1, ω 2 λ1 = 0.2 and ω 2 α = 0.1. Table 1 displays the results of the SAEM algorithm used for estimating the population parameters and their standard errors. The true values of the population parameters θ , the initial values θ 0 and the estimates θ are given together with their estimated standard errors (s.e.) and relative standard errors (r.s.e.).

Table 1 here

Table 1 shows the SAEM estimates are similar to the true values. On this particular example, the relative estimation error is less than 15% on the whole, except for parameter ω 2 γ2 (36%). We also note that the (relative) standard errors for each parameter are low, which is very encouraging, except for variance parameters ω 2 γ1 and ω 2 γ2 (48% and 41% respectively).

Figure 1 here

Figure 1 shows the sequences of estimated parameters (θ k ). One clearly sees that SAEM converges in very few iterations to a neighborhood of the "true" value used for simulating the data, even with a poor initialization. Moreover, it took only 6 on a laptop for estimating both the population parameters and the Fisher information matrix with this dataset.

Then we have estimated the individual parameters (ψ i ; 1 ≤ i ≤ 200) by computing the MAP estimates for each subject.

Finally, each individual state sequence was estimated with the Viterbi algorithm. As dealing with simulated datasets, the "true" state sequences are known. Even if this information is omitted during the whole inference process, true and estimated states can be compared. Figure 2 presents the results obtained with three typical subjects. On each graph, the (simulated) observations (daily seizures) are represented as a function of time (number of days). The true unknown states are displayed in the left column. The second column depicts the raw data, i.e. the only information available in the practice for inference. The right column displays the estimated states. We can observe a very good agreement between the true and the decoded states. -This first result is encouraging but it was obtained on a particular dataset. Then a Monte Carlo study should confirm the good behavior of the proposed methodology. 100 dataset were simulated using the same design and the same population distribution. Then, the population parameters and their standard errors were estimated with SAEM for each dataset. For m = 1, 2, . . . , 100, let θm be the estimated vector of population parameters obtained with the mth simulated dataset and let rse m be their respective estimated standard-errors. For each model parameter, we have computed the mean estimated parameter θ, the mean estimated relative standard error rse and the relative standard deviation of the estimated In Monolix, it is possible to choose the number of Markov chains used for the SAEM algorithm. Here, two Markov chains were used instead of only one chain. That allowed to slightly improve the convergence of the algorithm by reducing its stochastic behavior.

Each initial guess for the SAEM algorithm was randomly chosen for each Monte Carlo run.

The initial probability distribution of the hidden Markov chain was not estimated. It was assumed that π

(i) 1 = π (i) 2 = 1 2 .

Conclusion and perspectives

A brief state of the art on HMMs' applications to disease progression data shows that hidden Markov models are a reasonable modeling tool in this context. However, longitudinal data need models able to take account of the existing heterogeneity between individuals. This remark recently lead to the use of mixed hidden Markov models. However, related algorithms initially first tackled the population parameter estimation only including a small number of random effects. We suggested a new and complete inference methodology. The originality of our work consists in the use of the SAEM algorithm for estimating the model population parameters. A Monte Carlo study showed its good practical properties. More precisely, the SAEM algorithm converges to a neighborhood of the good parameter values in very few iterations even when the initial guess is poor. The estimation process is fast, even with large datasets. The impact of the data size and the theoretical properties of the algorithm keep to be rigorously studied in future works. The main problems to tackle would be a precise analysis of the convergence of the algorithm as well as the statistical properties of the maximum likelihood estimate in MHMM.

From a practical point of view those new models seem to offer very promising statistical applications. More precisely, mixed-effects hidden Markov models could help for a more finely analysis of clinical trials, when the collected data often consist of longitudinal count data and when one suspect several hidden states. Having a measure of interest, the most popular approach consists of mean comparisons, between groups of patients or treatment periods. However, the classical comparison methods could sometimes lead to improper conclusions. Assume a finite set of hidden stages give a plausible interpretation for the dynamics of the studied pathology, then an inappropriate choice for the statistical model (i.e. a model ignoring the transitions between distinct stages) does not catch enough information on the phenomenon observed. We could imagine a treatment effect occurs at transition level and a variation of the time spent in one specific hidden state once entering this state could constitute the only difference between treatment groups or between treatment periods, without modifying the observations' distribution in any state. Then, a simple comparison between means of observed outcomes would either fail in bringing to light a significant treatment effect or show an overestimated or underestimated treatment induced change. Hidden Markov models including the treatment group or the treatment dose as covariates in the definition of both the transition model and the observation model, would thus help to capture the true treatment effect and adequately locating it.

Contents 1 .

 1 Hidden Markov models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2. From hidden Markov models to mixed hidden Markov models . . . 4 3. Mixed hidden Markov models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 4. Application. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5. Conclusion and perspectives. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1. Hidden Markov models 1.1. Definition.

3. 2 . 2 . 3 .

 223 Most likely state sequences' decoding. -Once the individual parameters (Ψ i ) are estimated, each individual model can be considered separately and the optimal individual state sequences can be decoded using the Viterbi algorithm: Ẑi = argmax Zi p Z i |Y i , Ψi 4. Application 4.1. The model.

  i's Poisson parameters in state 1 and in state 2, with λ

Figure 2 here 4 . 2 . 1 .

 421 Figure 2 here 4.2.1. Monte Carlo study.-This first result is encouraging but it was obtained on a particular dataset. Then a Monte Carlo study should confirm the good behavior of the proposed methodology. 100 dataset were simulated using the same design and the same population distribution. Then, the population parameters and their standard errors were estimated with SAEM for each dataset. For m = 1, 2, . . . , 100, let θm be the estimated vector of population parameters obtained with the mth simu-

4. 3 .

 3 Technical remarks. -The proposed methodology for MHMM has been implemented in the Monolix software (http://software.monolix.org). All the numerical examples were performed with Monolix 3.1.

Figure 1 .

 1 Figure 1. Estimation of the population parameters: convergence of the SAEM algorithm.

Figure 2 .

 2 Figure 2. State sequences of three typical subjects. Left: the observed data and the true unknown states, center: the observed data without the states, right: the estimated states.

Table 1 .

 1 Estimation of the population parameters: the true values, the initial values, the estimations, their standard errors and relative standard errors.

Figure 3 .

 3 Figure 3. Estimation of the fixed effects: empirical distribution of the relative errors of estimations.

Figure 4 .

 4 Figure 4. Estimation of the variance parameters: empirical distribution of the relative errors of estimations.

Table 2 .

 2 Estimation of the population parameters: the true values, the means and the relative standard deviations of the estimated parameters, the mean estimated relative standard errors.

	γ 1	1.4 1.394	4	4
	γ 2	-1.4 -1.414	4	4
	λ 1	0.8 0.778	3	4
	α	2.3 2.292	3	3
	ω 2 γ1 ω 2 γ2 ω 2 λ1 ω 2 α	0.1 0.112 0.1 0.121 0.2 0.200 0.1 0.098	44 51 14 14	54 58 14 13

θ θ rsd( θ) (%) rse (%)

Here, θ denotes the true values of the population parameters.

Figure 3 here Figure 4 here Except for ω 2 γ1 and ω 2 γ2 , one can observe that the estimates show very little bias and small ranges, (table 2. The median REE for the fixed effects γ 1 , γ 2 , λ 1 and α remains between -5% and 5% (figures 3). The estimated variances ω λ1 and ω α are also very well estimated. The variances ω γ1 and ω γ2 are more difficult to estimate accurately. Indeed, the REE boxplots suggest quite important relative RMSE (root mean square errors) for those two parameters (44% and 51%) but this apparent estimation difficulty is in accordance with the estimated relative standard errors for those two parameters (54% and 58%).

More generally, one can remark the very good agreement between the estimated standard errors and the empirical standard deviations. The empirical standard deviations obtained from simulated data allow to evaluate the uncertainty of the estimated parameters. Of course, these empirical standard deviations cannot be computed in the practice when only one dataset is available and when the true population parameters are unknown. Nevertheless, one can have confidence with the estimated s.e. provided by the algorithm for evaluating the uncertainty of the estimated parameters.

These numerical results suggest that our algorithm produces unbiased and consistent population parameter estimates and standard errors in large databases. A theoretical study of the statistical property of the maximum likelihood population estimates is beyond the scope of this paper, but we will consider this issue in future works. On the other hand, more exhaustive studies should be led considering more difficult and more realistic contexts than the one here.


Ehab Al-Moubarak

Claire Simons
email: simonsc@cardiff.ac.uk

A homology model for Clostridium difficile methionyl tRNA synthetase: active site analysis and docking interactions

Keywords:

Treatment of C. difficile infection is one of the most difficult biomedical challenges. To develop novel antibacterials, researchers have been targeting bacte rial molecular functions that are essential for its growth. The methionyl tRNA synthetase (MetRS) is strictly required for protein biosynthesis and success was reported in developing antibacterials to inhibit this enzyme. The present study was aimed at building and analyzing a homology model for C. difficile MetRS in the context of drug design.

A homology model of C. difficile MetRS was constructed using Molecular Operating Environment (MOE) software. A. aeolicus MetRS was the main template while the query zinc binding domain was modeled using T. thermophilus MetRS. The model has been assessed and compared to its main template (Ramachandran, ERRAT and ProSA). The active site of the query protein has been predicted from its sequence using a detailed conservat ion analysis (ClustalW2). Using MOE software, suitable ligands were docked in the constructed model, including a C. difficile MetRS inhibitor REP3123 and the enzyme natural substrate, and the key active site residues and interactions were identified. These docking studies have validated the active site conformation in the constructed model and identified binding interactions.

Introduction

Clostridium difficile (C. difficile) is the most common cause of healthcare related diarrhoea in hospitals and health care institutions. Hospitalized patients on broad spectrum antibiotic treatment are considered those who are particularly at risk from developing the disease, an example is patients suffering from spinal cord injury who are subject to frequent antibiotic exposure and are hospitalized or using hospital facilities for long periods [START_REF] Richmond | Clostridium difficile infection in patients with SCI[END_REF]. It is estimated there has been a 40 percent increase in the C. difficile infection (CDI) incidence in the past two years [START_REF] Gerding | Deaths involving Clostridium difficile: England and Wales 2004-2008[END_REF] and news sources have shown that in 2008 8,324 C. difficile infection is currently treated with a range of therapies with meronidazole, routinely used as initial therapy for CDI [START_REF] Gerding | Deaths involving Clostridium difficile: England and Wales 2004-2008[END_REF]. Discontinuation of antibiotic treatment can restore the gut normal flora and stop the C. difficile overgrowth that results in the disease incident [START_REF] Gerding | Deaths involving Clostridium difficile: England and Wales 2004-2008[END_REF]. Nevertheless this is not always feasible especially when the patient has been treated for a serious condition, moreover, in the best situations only about 20-25 percent benefit from this intervention. Vancomycin is used in more severe cases and in subsequent episodes; combinations of antimicrobials have also been tried in some cases and treatment guidelines and algorithms are available to help in the management of CDI cases [START_REF] Gerding | Deaths involving Clostridium difficile: England and Wales 2004-2008[END_REF].

A potential target in the development of drugs for the treatment of C. difficile infections is aminoacyl-tRNA synthetase (AaRS). Aminoacyl-tRNA synthetases, can exist as either homodimers or monomers and catalyze the specific aminoacylation of aminoacid transfer RNAs [START_REF] Deniziak | Methionyl-tRNA synthetase[END_REF]. Class I aminoacylsynthetases consist of our domains: a catalytic Rossmann fold, a connective peptide (CP) inserted between the two halves of the Rossmann fold, a KMSK domain and a C-terminal -helix bundle domain, responsible for the recognition of the CAU anticodon of aminoacyl tRNAs [START_REF] Hountondji | Sequence similarities among the family of aminoacyl-tRNA synthetases[END_REF], whereas class II synthetases are characterized by a catalytic core of antiparallel -sheets and three consensus motifs [START_REF] Eriani | Partition of transfer-RNA synthetases into 2 classes based on mutually exclusive sets of sequence motifs[END_REF].

Although the natural substrates for AaRS enzymes are the same for all living organisms, the mode of binding of these amino acid substrates can be different depending on several factors including the active site conformation. In recent years increasing numbers of compounds have been reported to inhibit AaRS enzymes, many of them however inhibit human enzymes [START_REF] Hurdle | Prospects for aminoacyl-tRNA synthetase inhibitors as new antimicrobial agents[END_REF], which can cause toxicological issues. Success in identifying molecules that specifically inhibit bacterial AaRS enzymes has recently been achieved [START_REF] Hurdle | Prospects for aminoacyl-tRNA synthetase inhibitors as new antimicrobial agents[END_REF]. Success in clinical applications started with Mupirocin (marketed as Bactroban), which is the first AaRS inhibitor antibiotic that selectively inhibits bacterial isoleucyl-tRNA synthetase [START_REF] Hurdle | Prospects for aminoacyl-tRNA synthetase inhibitors as new antimicrobial agents[END_REF][START_REF] Boyce | MRSA patients: proven methods to treat colonization and infection[END_REF]. Mupirocin is now the first choice topical antibiotic for the treatment of methicillin-resistant Staphylococcus aureus (MRSA).

Several potent nanomolar S. aureus methionyl tRNA synthetase (MetRS) inhibitors were developed and found to have excellent antibacterial activity against Staphylococcal and Enterococcal pathogens [START_REF] Jarvest | Nanomolar inhibitors of Staphylococcus aureus methionyl tRNA synthetase with potent antibacterial activity against gram-positive pathogens[END_REF][START_REF] Finn | Discovery of a potent and selective series of pyrazole bacterial methionyl-tRNA synthetase inhibitors[END_REF][START_REF] Tandon | Potent and selective inhibitors of bacterial methionyl tRNA synthetase derived from an oxazolone-dipeptide scaffold[END_REF][START_REF] Finn | Identification of novel inhibitors of methionyl-tRNA synthetase (MetRS) by virtual screening[END_REF]. One of the quinolinone derivatives tested selective ly inhibited S. aureus MetRS and was found to be effective in vivo in rat [START_REF] Jarvest | Nanomolar inhibitors of Staphylococcus aureus methionyl tRNA synthetase with potent antibacterial activity against gram-positive pathogens[END_REF]. Very interesting is the generation of REP3123, which is a potent MetRS inhibitor of C. difficile [START_REF] Critchley | Spectrum of activity and mod of action of REP3123, a new antibiotic to treat Clostridium difficile[END_REF]. REP3123 was developed as an anti-microbial agent for treatment of C. difficile associated diseases (CDAD) and showed potent antimicrobial activity against C. difficile and Gram-positive aerobic bacteria; however it does not affect the normal gut flora content. REP3123 prevents C. difficile sporulation and therefore inhibits toxin production. The compound is an effective inhibitor of both C. difficile and S. aureus MetRS and does not inhibit human MetRS [START_REF] Critchley | Spectrum of activity and mod of action of REP3123, a new antibiotic to treat Clostridium difficile[END_REF].

The success of REP3123 as an inhibitor of C. difficile MetRS provides a useful ligand, in addition to known substrates of MetRS, to develop and evaluate a homology model.

There are no crystal structures of the AaRS enzymes from C. difficile, however comparative structure modeling can be used to construct a putative structure for the desired C. difficile MetRS enzyme that may facilitate inhibitor development. The construction and subsequent evaluation of a C. difficile MetRS homology model is described.

Methods

Construction of the C. difficile MetRS model

Homology searching

The C. difficile MetRS sequence was obtained from the ExPASy proteomics server [START_REF] Gasteiger | ExPASy: the proteomics server for in-depth protein knowledge and analysis[END_REF]. The enzyme sequence has the Uniprot identifier Q181D9 and was obtained from C. difficile strain 630, it has 645 residues [START_REF] Sebaihia | The multidrug-resistant human pathogen Clostridium difficile has a highly mobile, mosaic genome[END_REF]. A homology search was performed using NCBI BLAST2 service [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF] accessible from the ExPASy server, which was used to align the query against sequences in the Protein Data Bank (PDB) [START_REF]RCSB Protein Data Bank[END_REF] and thus the close homologous proteins were identified. The alignment parameters and thresholds used for screening candidate homologues were used with their default values and BLOSUM62 comparison matrix.

The Phylogeny server [18,[START_REF] Dereeper | Phylogeny.fr: robust phylogenetic analysis for the non-specialist[END_REF] was used to build a phylogenetic tree for these identified homologous proteins, the query sequence and other MetRS enzymes selected from different organisms.

Multiple sequence and structure alignment

The secondary structure of C. difficile MetRS together with the templates were determined using PSIPRED [START_REF] Jones | Protein secondary structure prediction based on position-specific scoring matrices[END_REF]. The query sequence was then aligned against the most homologues templates, Aquifex aeolicus MetRS, Thermus thermophilus MetRS, Pyrococcus Abyssi MetRS, Escherichia Coli MetRS using Clustal W2 [START_REF] Larkin | ClustalW and ClustalX version 2[END_REF] for identification of specifichelices, -sheets, coils and loops, and the common features and motifs: 'HIGH region', KMSKS motifs, zinc metal and ATP binding sites. Further evaluation of fold recognition and identification of template structures by structural alignment was performed using mGenTHREADER [START_REF] Lobley | pGenTHREADER and pDomTHREADER: New Methods For Improved Protein Fold Recognition and Superfamily Discrimination[END_REF].

3D Model Building

The molecular modeling experiments were performed on a RM Intel Core2 Duo CPU 2.2GHz running Windows XP using Molecular Operating Environment (MOE) 2008.09 molecular modeling software [23]. Homology models were built using MOE-Homology using AMBER99 forcefield [START_REF] Weiner | An all atom forcefield for simulations of proteins and nucleic acids[END_REF], which uses a dictionary to set the partial charges of atoms in amino acids. The final homology model was constructed using all the A. aeolicus MetRS (2CSX) crystal structure and the T. thermophilus MetRS zinc finger residues (from 1A8H). MOE-align was used to identify the zinc binding domain residues in the second template (T. thermophilus MetRS). The identified residues were selected to override the main template for building the model, the zinc and the Met-tRNA atoms were selected and included in the model as environment for induced fit, all other setting options were set to the default. Ten intermediate models were generated and the final model was taken as the Cartesian average of all the intermediate models. All the minimizations were performed with MOE until RMSD gradient of 0.05 kcal mol -1 Å -1 with the forcefield specified and the partial charges were automatically calculated.

Model validation

Stereochemical quality of the polypeptide backbone and side chains was evaluated using Ramachandran plots obtained from the RAMPAGE server [START_REF]RAMPAGE[END_REF]. Amino acid environment was evaluated using ERRAT plots [26]. ERRAT assesses the distribution of different types of atoms with respect to one another in the protein model. The ProSA tool [START_REF] Wiederstein | ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins[END_REF] was used to check the overall model protein structure for potential errors. Validation data for the templates A. aeolicus MetRS (2CSX) and T. thermophilus MetRS (1A8H) was used as the baseline to assess the respective models.

Docking

Ligand structures were obtained from the relevant complex crystal structures, REP3123 was built using MOE-builder [23] and then the ligand was energy minimized. These ligands were docked using the MOE-Dock with setting: Placement: Triangle Matcher; Rescoring 1: London dG, 30 poses were constructed for each compound and the best scoring model-ligand complexes were selected; the ligand interactions within these complexes were visualized using the MOE ligand interactions simulation.

Results and discussion

Homology models and validation

Initial screening for possible templates performed by running a BLAST analysis [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF] of the amino acid sequence of C. difficile MetRS, obtained from the ExPASy proteomics server [START_REF] Gasteiger | ExPASy: the proteomics server for in-depth protein knowledge and analysis[END_REF], against the PDB resolved structures, identified five structures for further consideration as possible templates (Table 1). For a structure to be considered it should ideally be wild-type, rather than mutant or engineered, have reasonable sequence identity with the C. difficile MetRS and preferably have the same function. The first ten native hits are bacterial MetRS enzymes, with the MetRS enzyme of Aquifex aeolicus (A. aeolicus) being the most appropriate template due to the high sequence identity (48%). The Thermus thermophilus (T. thermophilus) enzyme also had a good sequence identity of 39%. In addition Escherichia Coli (E. coli) MetRS (23% identity) and Pyrococcus Abyssi (P. Abyssi) MetRS (28% identity) can still be considered to have a reasonable sequence similarity. These MetRS enzymes const itute the dataset of suggested templates.

<Table 1>

The MetRS enzymes have common features and motifs: 'HIGH region', KMSKS motif, zinc metal and ATP binding sites (Table 2). Class I AaRS enzymes are characterized by a catalytic center built around a Rossmann fold where the two signature sequences HIGH and KMSKS are found [START_REF] Hountondji | Sequence similarities among the family of aminoacyl-tRNA synthetases[END_REF][START_REF] Landes | A structure-based multiple sequence alignment of all class I aminoacyl-tRNA synthetases[END_REF]. The HIGH motif is located in a loop between 1 and 1 elements, with the two histidine residues stacked upon each other [START_REF] Landes | A structure-based multiple sequence alignment of all class I aminoacyl-tRNA synthetases[END_REF]. The KMSKS peptide in located in a ----structure, named the stem contact (SC) fold [29,[START_REF] Kamijo | Analysis of conditions for KMSS loop in Tyrosyl-tRNA synthetase by building a mutant library[END_REF]. The HIGH region has been implicated in amino acid activation and the KMSKS in docking the acceptor stem of tRNA.

The positions of these features in different species are varied or shifted by a few residues but are still present in all of them [START_REF] Hountondji | Sequence similarities among the family of aminoacyl-tRNA synthetases[END_REF].

<Table 2>

To obtain more information regarding the best possible template, the Phylogeny server [18,[START_REF] Dereeper | Phylogeny.fr: robust phylogenetic analysis for the non-specialist[END_REF] was used to construct a phylogenetic tree using sixteen more MetRS enzyme sequences in addition to the query and the above four suggested MetRS templates. A diverse range of species was considered including mammals, fish, yeast and bacteria so as to appreciate the relative distances between the suggested four templates and the query sequence (Fig. 1).

<Figure 1>

The enzyme sequences of S. aureus and C. acetobutylicum were the closest to the query sequence in the phylogenetic tree. However, by considering the next closest branch, it was interesting to observe that the enzyme sequences of A. aeolicus and T. thermophilus were the next closest. The other two suggested templates of E. coli and P. Abyssi were diverting in branches relatively far from the query sequence.

Multiple sequence and structural alignments

ClustalW2 [START_REF] Larkin | ClustalW and ClustalX version 2[END_REF] was used to align the suggested template sequences and the query sequence of C. difficile MetRS (Fig. 2).

<Figure 2>

Very few gaps in alignment were observed when considering the C. difficile, T. thermophilus and A. aeolicus MetRS sequences alone until approximately residue 500. However, consistent gaps were observed between these three sequences in respect to the other two suggested templates of E. coli and Pyrococcus Abyssi MetRS sequences. The difference between the two groups was in agreement with the relations described by the phylogenetic tree (Fig. 1).

HIGH and KMSKS motifs are identified in all suggested templates and query sequence (Fig.

2).

To compare structural alignment, the C. difficile MetRS protein sequence was used as a template for the mGenTHREADER program [START_REF] Lobley | pGenTHREADER and pDomTHREADER: New Methods For Improved Protein Fold Recognition and Superfamily Discrimination[END_REF], to search the Protein Da taBank (PDB) for proteins that are likely to contain similar structural motifs. The mGenTHREADER program is a recommended method for fold recognition and identification of distant homologues and makes use of profile-profile alignments and predicted secondary structure (using PSIPRED) as inputs. The GenTHREADER method also runs PSI-BLAST on the target sequence before the search. A. aeolicus (45% identity) and T. thermophilus (38% identity) are likely the closest proteins structurally (Table 3).

<Table 3>

Template analysis therefore suggested A. aeolicus MetRS as the main template, however the A. aeolicus MetRS crystal structure (2CSX) lacks the zinc binding domain involving residues 124 to 158. Employing the missing sequence folds from T. thermophilus MetRS rectified this problem. This multi-template approach is quite appropriate here because the missing residues are found clustering in the structure just around the zinc. Visualization of both templates in MOE by superimposition has shown that the T. thermophilus (1A8H) and A. aeolicus (2CSX) MetRS structures overlap well with the zinc binding domain of 1A8H clearly identified (Fig. 3).

<Figure 3>

Both the T. thermophilus and A. aeolicus MetRS sequences are too short to align with the C. difficile MetRS sequence C terminus. The C terminal structure of P. Abyssi MetRS (1MKH) seems to be the closest in terms of sequence identity and secondary structure; however a considerable gap in alignment of the C terminal sequence would prevent using it as a template for modeling the C terminus alone. It has been observed that many crystal structures of MetRS have an incomplete C terminal structure. For example, the structure of P. Abyssi MetRS (1RQG) lacks the whole C-terminus. The C-terminal domain is less well conserved overall among the MetRS enzymes, which makes modeling this domain challenging [START_REF] Deniziak | Methionyl-tRNA synthetase[END_REF][START_REF] Spencer | Characterization of the human mitochondrial methionyl-tRNA synthetase[END_REF].

Studies of the C-terminal domains of P. Abyssi MetRS [START_REF] Crepin | Three-dimensional structure of methionyl-tRNA synthetase from Pyrococcus abyssi[END_REF] and E. coli MetRS [START_REF] Mechulam | Crystal structure of Escherichia coli methionyl-tRNA synthetase highlights species-specific features[END_REF] have indicated that the enhancement of tRNA-binding affinity conferred by the C-terminal domain is only expressed when the enzyme is in the dimeric state and, for E. coli MetRS, the enzyme remains active in monomeric form after removal of a C-terminal region required for dimerization [START_REF] Mechulam | Crystal structure of Escherichia coli methionyl-tRNA synthetase highlights species-specific features[END_REF]. By weighing the risks and assessing the need for modeling the C terminus of the query, it was decided not to incorporate an unsatisfactory C-terminus crystal structure in the multiple template combination, which may adversely affect the rest of the model. PSIPRED [START_REF] Jones | Protein secondary structure prediction based on position-specific scoring matrices[END_REF] secondary structure prediction from the mGenTHREADER program [START_REF] Lobley | pGenTHREADER and pDomTHREADER: New Methods For Improved Protein Fold Recognition and Superfamily Discrimination[END_REF] for C. difficile MetRS clearly showed the high helix content predicted throughout the sequence.

The C-terminus had a higher degree of coils and strands predicted than for the rest of the protein sequence, this in addition to the poor alignment observed for the C-terminus would make the C terminus a challenging sequence to model. The HIGH region motif (Table 2) is a sequence of about ten residues and was found in most of the homologous enzymes at the sequence positions just before the 22 nd residue in the N termini of these enzymes. This particular area of the query sequence was expected by PSIPRED to fold in coils and strands, this holds for A. aeolicus where the resolved structure shows residues 6-15 form a strand very close to where the HIGH region motif was suggested. The KMSKS motif is another known motif among the MetRS enzymes described. The motif contains the ATP binding site and has been found in the query sequence aligned to the T. thermophilus and A. aeolicus sequences with less than a three position differences in both cases (Table 2). The A. aeolicus 3D structure contains strands and turns from position 283 until at least 301. This is very much in agreement with the predictions of PSIPRED for the query sequence.

Zinc binding

Zinc is an absolute requirement for enzyme activity [START_REF] Xu | Identification of the metal ligands and characterization of a putative zinc finger in methionyl-tRNA synthetase[END_REF]. There are two steps associated with the enzyme activity catalyzing the transfer of the Met amino acid into the cognate tRNA: tRNAMet aminoacylation and ATP-pyrophosphate exchange activities. It is certain now that the zinc should be in the actual structure of the enzyme and that it is most likely bound to four cysteine residues in the query protein and hence the intended query model. This was supported by the alignment of C. difficile MetRS with the suggested templates where several conserved Cys residues and a histidine residue (E. coli) are identified in the likely region of this domain in the query sequence (Fig. 2).

3D Homology model

Homology models, constructed using all the A. aeolicus MetRS (2CSX) crystal structure and the T. thermophilus MetRS zinc finger residues (1A8H), were generated using Molecular Operating Environment (MOE) [23] software as described in the experimental section. Using combined templates to generate the C. difficile MetRS model allowed identification of the cysteine residues binding the zinc ion. These cysteine residues were all near the position of the zinc, at least three of which had the sulfur directly pointing towards the zinc. For the fourth cysteine residue side chain manipulation using MOE-rotamer explorer was employed to improve the zinc binding. The sulfur atoms binding the zinc at the constructed model were from Cys128, Cys131, Cys146 and Cys149 (Fig. 4), also noted in close proximity was Asp148 that could also chelate with the zinc ion.

<Figure 4>

The Met tRNA sequence of C. difficile, available from the National Centre for Biotechnology Information [START_REF] Xu | Identification of the metal ligands and characterization of a putative zinc finger in methionyl-tRNA synthetase[END_REF]35] consists of 73 bases, similar to that of A. aeolicus. Moreover, the CAU anticodon positions are also the same. RNAfold [START_REF] Mathews | Expanded sequence dependence of thermodynamic parameters provides robust prediction of RNA secondary structure[END_REF], a web server for predicting the RNA secondary structure, indicated that the A. aeolicus Met tRNA secondary structure is quite similar to the predicted C. difficile Met tRNA structure using the RNA fold (Fig. 5) therefore the tRNA chain of A. aeolicus MetRS was used in the query MetRS homology modeling. <Figure 5> Asn356, Arg360 and Trp430 of C. difficile MetRS are fully involved in binding the anticodon part of Met-tRNA (Fig. 6). The corresponding residues in the MetRS enzymes of A. aeolicus, T. thermophilus and E. coli play the same role. Asn356 and Arg360 form hydrogen bonds with 35A and 34C. As seen in the A. aeolicus MetRS complex the Trp residue (422 in A. aeolicus, Trp430 in the query sequence) is not involved in direct interaction with the anticodon 36U but interacts by hydrogen bond with a Phe residue (492 in A. aeolicus and Phe508 in the query sequence). While these findings are predictable by literature mining along with sequence alignments, having found that the model met these e xpectations can give some confidence to additional observations: Lys367 and Lys434 of C. difficile MetRS interact with 35A and 34C while Leu507 and Arg510 interact with 36U (Fig. 7).

<Figure 6>

The final C. difficile MetRS model (Fig. 7A) contained 513 residues from the 645 residues of the sequence; this was the result of disabling the C-terminus. The secondary structure of these residues involved twenty -helices and 14 -sheets (Table 4). This model also contained a zinc metal and A. aeolicus tRNA.

<Figure 7> <Table 4>

The Rossmann fold domain is formed by two polypeptide segments (residues 1-115 and 229-289, colored orange in Fig. 7B). Consistent with other AaRS enzymes [29,[START_REF] Cavarelli | L-arginine recognition by yeast arginyl-tRNA synthetase[END_REF][START_REF] Nureki | Architectures of class-defining and specific domains of glutamyl-tRNA synthetases[END_REF], the last residue of the Rossmann fold, His289 in the 12 strand, forms a hydrogen bond the backbone N-H group of Val259 in the 11 strand and is followed by a glycine residue. The HIGH motif of C. difficile MetRS (His20-Ile21-Gly22-His23) is located between 2 and 1, with the two histidine residues stacked upon each other.

Inserted between the two Rossmann fold polypeptide segments is the connective polypeptide (CP) domain, composed of an antiparallel sheet core (cyan color, Fig. 7B), a MetRSspecific -strand insertion (green, Fig. 7B), zinc finger (yellow, Fig. 7B) and antiparallel 2helical structure (blue, Fig. 7B). The CP core consists of 5, 6, 7, 8 and 9, between the second and third strands ( 6 and 7) is the zinc finger with one Zn 2+ ion coordinated by sulfur atoms (Cys128, Cys131, Cys146 and Cys149, Fig. 4). In T. thermophilus MetRS the zinc binding domain consists of an antiparallel -sheet insertion involving 3 -strands [29]. In this C. difficile MetRS model, the zinc binding domain consists of a loop between 5 and 6 with an -helix ( 5) positioned in the middle of the four zinc-coordinating cysteine residues.

Between the second and third -strands ( 6 and 7), two antiparallel -helices ( 6 and 7)

are inserted (Fig. 7B), and between the last antiparallel -sheet domain and the C-terminal Rossmann fold domain, a -strand ( 10, Val225-Ile224, green in Fig. 7B) is inserted in an antiparallel manner. This -strand insertion is considered to be a MetRS specific element [29].

The region connecting the N-terminal Rossmann fold (orange, Fig. 7B) with the C-terminal -helix bundle domain (magenta, Fig. 7B) is the SC-fold domain (red, Fig. 7B). The SC-fold domain, which is the RNA binding domain, has a ----( 13-12-13-14-14)

topology. This homologous topology is well conserved and characteristic of the monomeric class Ia and Ib synthetase enzymes [29,[START_REF] Cavarelli | L-arginine recognition by yeast arginyl-tRNA synthetase[END_REF][START_REF] Nureki | Architectures of class-defining and specific domains of glutamyl-tRNA synthetases[END_REF]. The KMSKS (Lys298-Met299-Ser300-Lys291-Ser292) motif is situated in the SC-fold domain in a loop between 13 and 12.

The anticodon binding -helix bundle domain (magenta, Fig. 7B) consists of five antiparallel -helices [START_REF] Sebaihia | The multidrug-resistant human pathogen Clostridium difficile has a highly mobile, mosaic genome[END_REF][START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF][START_REF]RCSB Protein Data Bank[END_REF](18)[START_REF] Dereeper | Phylogeny.fr: robust phylogenetic analysis for the non-specialist[END_REF]. The amino acid residues involved in anticodon recognition are clustered on one face of the -helix bundle domain, with Asn356, Arg360 and Lys367 in the 14 helix and Trp430 and Lys434 in the 16 helix.

Model validation

To evaluate the quality of the modeled structures, the lowest energy model generated was subjected to a number of checks. Stereochemical quality was assessed us ing Ramachandran plots, using the Cambridge RAMPAGE server [START_REF]RAMPAGE[END_REF], and amino acid environment was assessed using ERRAT [26]. Overall protein structure was assessed using ProSA [START_REF] Wiederstein | ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins[END_REF].

Validation results would suggest that the model performed well in terms of mainchain stereochemistry and amino acid environment. In the Ramachandran plot, a total of 97.9% of the residues were in the allowed region, which was comparable with the main template 2CSX (97.1%), indicating that the backbone dihedral angles and in the model were reasonably accurate. There were 11 outliers in the model however they were far from the active site and do not significantly contribute to its function. The ERRAT results showed that the overall quality factor of the protein was good with a score of 86.139 compared with the main template 2CSX (84.862). It is surprising for the model to perform better than the template, however the template was missing the Zn-binding domain residues. The Zn-binding template 1A8H scored better in both Ramachandran (99.4 % of residues in the allowed region, 3 outliers) and ERRAT (98.374) evaluation, therefore the combined template approach has incorporated some improvements in stereochemical quality compared with the 2CSX template alone.

<Figure 8>

ProSA analysis provides output in two plots: the first plot (Fig. 8A)

Active site analysis and docking

By using the ClustalW2 multiple sequence alignment [START_REF] Larkin | ClustalW and ClustalX version 2[END_REF] and the MOE alignment service Docking the natural substrate, methionine, identified key residues with His53, Asp51 and Lys56 consistently forming H-bonding interactions with the carboxylic acid moiety of methionine. Ile12 was also found to hydrogen bond with this moiety in different dockings and the residues Ile224, Val226 and Ala230 forming a hydrophobic pocket (Table 5).

<Table 5>

The MSP ligand was docked in the C. difficile MetRS model. It was observed that the residues around the docked ligand remained almost the same in all the dockings (Table 5).

The following residues interact with the ligand in different docking attempts: Ile12, His20, 5).

Asp51
2FM was found to dock in two sites within the active site, the first was the site involving interaction with Asp51 and Lys56 seen previously with methionine, MSP and MF3 and the second site involved interaction with Glu55, Ser133 and Tyr225. Many conserved residues in the N terminal and in positions from 224 to 240 and from 255 to 268 are also found in the 2FM dockings (Table 5). Docking of 1-amino-3-methylsulfanyl-propyl)-phosphonic acid (MPH), found in the resolved complex structure of E. coli MetRS (1P7P) [START_REF] Crepin | Use of analogues of methionine and methionyl adenylate to sample conformational changes during catalysis in Escheric hia coli methionyl-tRNA synthetase[END_REF] resulted in diverse ligand interactions but comparable with those found when docking 2FM, that is MPH docked in two distinct binding sites within the active site. Again the first was the site involving interaction with Asp51 and Lys56 seen previously with MSP and MF3 and the second site involved interaction with Glu55, Ser133 and Tyr225 (Table 5). REP 3123 is the only ligand found in the literature that has been described as a selective C.

difficile MetRS inhibitor [START_REF] Critchley | Spectrum of activity and mod of action of REP3123, a new antibiotic to treat Clostridium difficile[END_REF]. The antimicrobial effect is described and the chemical structure of the inhibitor. It was confirmed in their study that REP 3123 is a selective and potent competitive inhibitor of C. difficile MetRS, which implies that the ligand has features that enable it to strongly bind at the active site of C. difficile MetRS. REP 3123 is not a competitive inhibitor of ATP and the presence of ATP facilitated inhibition [START_REF] Critchley | Spectrum of activity and mod of action of REP3123, a new antibiotic to treat Clostridium difficile[END_REF], therefore it could be deduced that the ligand binds by a similar mechanism as methionine. REP 3123 does not inhibit E. coli MetRS however it does inhibit S. aureus MetRS [START_REF] Critchley | Spectrum of activity and mod of action of REP3123, a new antibiotic to treat Clostridium difficile[END_REF], which may be explained by the close position of C. difficile and S. aureus MetRS indicated by the phylogenic tree (Fig. 1).

<Figure 9>

Docking of REP3123 indicated that this ligand docks across both binding sites, the first pocket is formed by Asp51, Ile12, Ala230 and Trp227 and the second pocket, linked by Lys56 and Tyr14, is formed by Glu55, Ser133, Tyr225 and Gln139 (Fig. 9). Interestingly a pharmacophore virtual screening study by Finn et al [START_REF] Finn | Discovery of a potent and selective series of pyrazole bacterial methionyl-tRNA synthetase inhibitors[END_REF], employed to develop S. aureus

MetRS inhibitors, identified Asp51 as a key residue with optimal inhibitory activity observed on formation of two hydrogen bond interactions between the ligand and the critical Asp51 group. Another study investigating the effect of resistance mutations on S. aureus MetRS activity found that Gly54 (G54S) and Gly223 (G223C) mutations were the most deleterious single mutants with markedly impaired enzymatic function [START_REF] Green | Inhibition of methionyl-tRNA synthetase by REP8839 and effects of resistance mutations on enzyme activity[END_REF].

  figure shows a four times increase in the deaths associated with CDI compared with 2001 [3].

  shows the local model quality by plotting energies as a function of amino acid sequence position, in general, positive values correspond to problematic or erroneous parts of the input structure. The second plot (Fig. 8B) shows overall model quality from which the z-score is calculated, its value is displayed in a plot that contains the z-scores of all experimentally determined protein chains in the current PDB, a negative score indicates a good model whereas a positive score would indicate errors. The z-score of the model was -10.57 compared with 2CSX (-11.79) and 1A8H (-12.72). Furthermore, superimposition of the model with the main template 2CSX showed a very low RMSD (0.816 Å) suggesting a high similarity between them, superimposition with 1A8H resulted in a higher RMSD (2.198 Å), however as only the zinc finger residues were used from 1A8H this was not expected to have a detrimental effect on the model. Overall based on Ramachandran plot, ERRAT and ProSA a good model, in terms of quality of the local geometry of the backbone and side-chain conformations, for C. difficile MetRS was generated. Further validation of the active site architecture and key binding interaction were achieved by substrate and ligand docking experiments.

[

  23], the positions of the active site conserved residues in the query sequence, A. aeolicus, E. coli and T. thermophilus MetRS can be identified. Docking suitable ligands into the constructed model was used to validate the predicted active site. The docked ligands include the MetRS natural substrate (methionine), several E. coli MetRS inhibitors 5'-O-[N-(Lmethionyl)-sulfamoyl] (MSP), difluoromethionine (2FM) and trifluoromethionine (MF3), (1amino-3-methylsulfanyl-propyl)-phosphonic acid (MPH) and a C. difficile MetRS inhibitor (REP3123), no site restriction was specified in the MOE docking procedures.

  , Lys56, Val226, Trp227 and Trp291 with H-bonding interactions most frequently observed with Asp51, Lys56 and Ile12. Difluoromethionine (2FM) and trifluoromethionine (MF3) are only available bound to the E. coli MetRS [39]. MOE dockings of MF3 showed that Asp51 of the C. difficile MetRS model is consistently interacting with the ligand. Ile12, Asp51, His267, Val226, Trp227, Leu233 and Ile263 are in the model binding pocket surrounding the ligand and are very close to the expectations derived from E. coli MetRS complex with MF3 (Table

  The binding mode of ligands in the C. difficile MetRS model was comparable with the binding of identical ligands in E.coli crystal structures (PDB identifier: 1F4L (methionine); 1PFY (MSP); 1PFV (2FM); 1PFW (MF3); 1P7P (MPH). In these crystal structures the key H-bonding interactions involved the conserved residues Leu13 and Asp52, with the residues Ala256 (Ala230 in C.difficile), Pro257, Tyr260, Ile297 (Ile263 in C.difficile) and His301 (His267 in C. difficile) forming hydrophobic interactions. The overlap in binding residues between identical ligands and the C. difficile MetRS model and E. coli MetRS crystal structures demonstrates that the modeled structure is of sufficient quality to correctly reproduce ligand binding poses in docking simulations.

  Conf = Description of confidence level, CERT: p -value < 0.0001; HIGH: p-value < 0.001 Net Score = Raw score fro m SVM; p-Value = Probability of false positive; PairE = Pairwise energy for model; So lv E = So lv ation energy for model; Aln Score = Sequence align ment score; Aln Len = Length of align ment; Str Len = Length of PDB entry; Seq Len = Length of target sequence; % ID = percentage identical residues between target protein and PDB protein; PDB_ID = PDB identifier (+ chain code + do main code in CATH format)

Fig. 2 Fig. 3 Fig. 4 Fig. 5 AFig. 7 (Fig. 8 Fig. 9

 2345789 Figure captions

Table 1

 1 The high sequence identity (50%) between S. aureus MetRS and C. difficile MetRS may suggest that the proximity of The five optimal crystal templates identified in the C. difficile MetRS BLAST resultsThe BLAST score for an align ment is calcu lated by summing the scores for each aligned position and the scores for gaps b (Number of identical residues)/ (length of sequence fragment identified by BLAST)

	Glu55 to Gly54 in the binding pocket would indicate the importance of Glu55 in binding
	interactions and inhibition of enzyme function. The ability of REP3123 to interact with both
	binding pockets and form H-bonding interactions with several residues including Asp51,
	Ile12, Lys56, Glu55, Ser133 and Tyr225 may account for the high selectivity and inhibitory
	activity against C. difficile MetRS.

a

Table 2

 2 Key residues of the suggested templates

	PDB	HIGH	KMSKS	Zinc Binding Residues	ATP binding
		region motif	motif
	2CSX-A	14-18	295-299	129, 132, 147, 150	298
	1A8H-A	12-22	297-301	127, 130, 144, 147	300
	1RQG-A	11-21	344-348	143, 146, 156, 159	347
	1F4L-A	15-25	333-337	146, 149, 159, 162	336

Table 3

 3 mGenTHREADER results for C.difficile MetRS query

	Conf Net	p-Value E Pair	E Solv	Aln	Aln	Str	Seq	%	PDB	Description
	Score				Score	Len	Len	Len	Iden	code
	CERT 178.758 5e-17	-650.4	-36.9 853.3 462	465	645	45.2	2ct8A0	Aquifex aeolicus MetRS
	CERT 174.608 1e-16	-578.6	-27.4 857.0 504	544	645	21.0	3h9cA0	Escherichia coli MetRS
	CERT 172.541 2e-16	-447.3	-26.9 867.0 545	606	645	22.4	1rqgA0	Pyrococcus abyssi MetRS
	CERT 171.741 2e-16	-518.2	-39.1 825.0 490	500	645	37.8	2d5bA0	Thermus thermophilus MetRS
	CERT 135.154 1e-12	-419.6	-14.0 642.0 561	862	645	17.7	1gaxA0	Thermus thermophilus ValRS
	CERT 106.519 9e-10	-322.7	-15.3 469.0 521	821	645	14.0	1ileA0	Thermus thermophilus IleRS
	CERT 103.727 2e-09	-256.6	-6.5	476.0 580	917	645	13.2	1ffyA0	Staphylococcus aureus IleRS
	CERT 97.544	7e-09	535.5	13.1	389.0 362	399	645	9.8	3c8zA0	Mycobacterium smegmatis CysRS
	CERT 91.009	3e-08	-151.6	3.2	457.0 444	628	645	9.6	2zueA0	Pyrococcus horikoshii ArgRS
	CERT 90.077	4e-08	-318.0	-3.9	396.0 440	810	645	10.9	2v0cA0	Thermus thermophilus LeuRS

Table 4

 4 Secondary structure of C.difficile MetRS model

	Residues	2 o structure	Do main
	PHE6-TYR9	1
	ILE12-TYR13	2
	ILE21-ARG38	1
	VA L44-THR50	3	Rossmann fold
	GLU55-ILE62	2
	GLU68-M ET86	3
	ASP92-A RG95	4
	LYS99-GLN115	4
	ILE118-TRY127	5	Antiparallel -sheet
	GLU137-GLN139	5	Zinc-finger
	LEU155-PHE163	6	Antiparallel -sheet
	GLU169-LYS177	6	2 helical structure
	ARG187-LYS198	7	CP
	LEU203-SER204	7
	THR206-A RG207		Antiparallel -sheet
	LYS215-VA L216	9
	VA L223-ILE224	10	MetRS-specific insertion
	ASP229-THR236	8
	GLU246-LYS250	9
	GLN257-GLY260	11
	LYS261-VA L264	10	Rossmann fold
	ILE270-A LA277	11
	VA L286-HIS289	12
	ILE292-LEU293	13
	PRO309-A RG315	12
	THR320-LEU326	13	Stem contact (SC) fold
	GLY335-SER336	14
	HIS339-TYR268	14
	LEU386-ASP402	15
	PHE406-TYR423	16
	PRO429-LYS434	17	-Helix bundle
	ASP443-ILE460	18
	ASN465-ASP481

Table 5

 5 Representative dockings into C. difficile MetRS model

	Substrate	Docking	H-bonds
	Methionine		Asp51,
			His53, Lys56

Conclusions