Xavière Panhard

Adeline Samson

Extension of the SAEM algorithm for nonlinear mixed models with two levels of random effects

Keywords: Multilevel nonlinear mixed effects models, SAEM algorithm, Multiple periods, Cross-over trial, Bioequivalence trials

This article focuses on parameter estimation of multi-levels nonlinear mixed effects models (MNLMEMs). These models are used to analyze data presenting multiple hierarchical levels of grouping (cluster data, clinical trials with several observation periods,...). The variability of the individual parameters of the regression function is thus decomposed as a between-subject variability and higher levels of variability (for example within-subject variability). We propose maximum likelihood estimates of parameters of those MNLMEMs with two levels of random effects, using an extension of the SAEM-MCMC algorithm. The extended SAEM algorithm is split into an explicit direct EM algorithm and a stochastic EM part. Compared to the original algorithm, additional sufficient statistics have to be approximated by relying on the conditional distribution of the second level of random effects. This estimation method is evaluated on pharmacokinetic cross-over simulated trials, mimicking theophyllin concentration data. Results obtained on those datasets with either the SAEM algorithm or the FOCE algorithm (implemented in the nlme function of R software) are compared: biases and RMSEs of almost all the SAEM estimates are smaller than the FOCE ones. Finally, we apply the extended SAEM algorithm to analyze the pharmacokinetic interaction of tenofovir on atazanavir, a novel protease inhibitor, from the ANRS 1 107-Puzzle 2 study. A significant decrease of the area under the curve of atazanavir is found in patients receiving both treatments.

Introduction

The use of non-linear mixed effects models (NLMEMs) increases in several fields such as agronomy, forestry, clinical trials, population pharmacokinetics (PK) and pharmacodynamics or viral dynamics to model longitudinal data. In some settings, data can present multiple hierarchical levels of grouping, leading to multiple nested levels of variability. For instance, we may study patients that are grouped in medical services that are themselves grouped into hospitals. In this article, we call multilevel non-linear mixed effects models (MNLMEMs) the models that describe such data. MNLMEMs represent a natural extension of models with single variability level, and they have recently been subject to a great deal of attention in statistical literature. In the field of forestry, Hall and Clutter [START_REF] Hall | Multivariate multilevel nonlinear mixed effects models for timber yield predictions[END_REF] analyze longitudinal measures of yield and growth that are measured on each tree within a plot. In the field of agronomy, Rekaya et al. [START_REF] Rekaya | Hierarchical nonlinear model for persistency of milk yield in the first three lactations of holsteins[END_REF] consider milk yield data where each cow is observed longitudinally during its first three lactations. Jaffrézic et al.

(3) perform genetic analyses of growth measurements in beef cattle acknowledging the fact that several cows come from the same sire. Another example is population PK, where concentration measurements may be taken with several patients over several distinct time intervals, that are often named periods or occasions. That grouping pattern is used for instance in PK cross-over trials.

In NLMEMs with only one level of variability, often corresponding to between-subject variability, the analysis results in the estimation of the fixed effects parameters and of the betweensubject variability of the parameters, also called inter-subject variability. When there is more than one level of grouping, the higher levels of variability can be estimated. In the specific case where the second level of grouping corresponds to multiple periods of measurement, this variability is called within-subject variability (or intra-subject variability, or inter-occasion variability), and corresponds to the variation of the individual parameters across the different study periods or units. In the context of pharmacokinetics, Karlsson et al. (4) demonstrate the importance of modeling this second level of variability in two-levels NLMEMs. They show that neglecting it resulted in biased estimates for the fixed effects.

The parameter estimation of NLMEMs is not trivial because the likelihood of NLMEMs cannot be expressed in a closed form due to the non-linearity of the regression function in the random effects. Therefore, several estimation methods have been proposed. The First Order Conditional Estimates (FOCE) algorithm performs a first order linearization of the regression function with respect to the random effects (5; 6). The implementation of the FOCE algorithm in the NONMEM software and in the nlme function of Splus and R enables the estimation of both between-and within-subject variabilities. From our practice, the main drawback of this method is however that it does not always converge when one estimates simultaneously the betweenand the within-subject variabilities on several parameters. Furthermore, this linearization-based method cannot be considered as fully established in theory. For instance, Vonesh [START_REF] Vonesh | A note on the use of Laplace's approximation for nonlinear mixed-effects models[END_REF] and Ge et al. [START_REF] Ge | An approximate likelihood approach to nonlinear mixed effects models via spline approximation[END_REF] give examples of specific designs resulting in inconsistent estimates, such as when the number of observations per subject does not increase faster than the number of subjects or when the variability of random effects is too large.

Several estimation methods have been proposed as alternatives to linearization algorithms.

A common method to handle numerical integration is the adaptative Gaussian quadrature (AGQ) method. An estimation algorithm of NLMEM parameters based on this classical AGQ method has been proposed by Pinheiro and Bates [START_REF] Pinheiro | Approximations to the log-likelihood function in the non-linear mixedeffect models[END_REF] and is implemented in the SAS procedure NLMIXED [START_REF] Wolfinger | Laplace's approximations for non-linear mixed-effect models[END_REF]. However, the AGQ method requires a sufficiently large number of quadrature points implying an often slow convergence with very high computational time. Furthermore, two-levels NLMEM can be implemented in the NLMIXED procedure, but to our knowledge, the convergence is difficult to obtain in practice [START_REF] Jaffrezic | Genetic analysis of growth curves using the SAEM algorithm[END_REF]. Improvements upon this method are thus needed. The second alternative to linearization is the use of the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] in order to estimate models with missing or non-observed data such as random effects. Because of the nonlinearity of the model, stochastic versions of the EM algorithm have been proposed. Wei et al. [START_REF] Wei | Calculating the content and boundary of the highest posterior density region via data augmentation[END_REF]; Walker [START_REF] Walker | An EM algorithm for non-linear random effects models[END_REF] and Wu (14) propose MCEM algorithms, with a Monte-Carlo approximation of the E-step. However the MCEM algorithm may have computational problems (i.e slow or even non convergence). As an alternative to address computational problems, a stochastic approximation version of EM (SAEM) has been proposed in (15; 16), which requires the simulation of only one realization of the missing data for each iteration, sub-stantially reducing the computation time. Kuhn and Lavielle (16) propose to combine the SAEM algorithm with a Monte-Carlo Markov Chain (MCMC) procedure adapted to the NLMEMs, and prove that the resulting estimates are convergent and consistent.

To date, none of the EM-based algorithms are directly applicable to the case of multilevel NLMEMs and have to be adapted. The objective of this paper is to extend the SAEM algorithm to MNLMEMs with two levels of variability: both E and M steps need to be adapted to integrate higher levels of random effects. We also propose estimates of the likelihood and of the Fisher information matrix. We evaluate this algorithm on a PK example, more precisely a two-periods one-sequence cross-over trials simulated mimicking theophyllin concentration data [START_REF] Pinheiro | Approximations to the log-likelihood function in the non-linear mixedeffect models[END_REF]. We also apply the SAEM algorithm to the PK interaction of two HIV molecules (tenofovir and atazanavir) from a PK substudy of the ANRS 107-Puzzle 2 trial.

After describing the model and notations (section 2), section 3 describes the SAEM algorithm. Section 4 reports the simulation study and its results. In Section 5, we study the PK interaction of tenofovir on atazanavir in HIV patients. Section 6 concludes the article with some discussion.

Models

Let us denote y ijk the observation in unit k (k = 1, . . . , K) for subject i (i = 1, • • • , n) and at time t ijk (j = 1, • • • , n ik ). For instance, the different units can be the different periods in the case of PK trials, or the different parents in the case of genetic analyses. We assume, as a known fact, two nonlinear functions f and g such that the two-levels non-linear mixed effects model linking observations to sampling times can be written as:

y ijk = f (t ijk , φ ik ) + g(t ijk , φ ik )ε ijk , ε ijk ∼ N (0, σ 2 ),
where φ ik is the p-vector of the parameters of subject i for unit k and ε ijk is the measurement error. We hypothesize that the errors ε ijk given φ ik are mutually independent. We assume that the individual parameters φ ik are random vectors and that for each unit k, φ ik can be broken into:

φ ik = µ + β k + b i + c ik , (1)
b i ∼ N (0, Ω), c ik ∼ N (0, Ψ),
where µ + β k is the mean value for unit k, b i is the random effect of size p of subject i, and c ik is the random effect of size p of subject i and unit k. To ensure the identifiability of the parameters, we assume that β 1 = 0, ie µ is the mean of the first unit and β k represents the difference (or effect) of the kth unit in comparison to this first unit. The random effects (b i )

and (c ik ) are assumed to be mutually independent. The total variance of the parameters is thus broken into a between-subject variance Ω and a within-subject variance Ψ. Finally, the individual parameters pK-vector φ i = (φ i1 , . . . , φ iK ) of subject i is distributed with a Gaussian distribution with mean vector (µ, µ + β 2 , . . . , µ + β K ) and a pK × pK covariance matrix Γ equal to

Γ =          Ω + Ψ Ω . . . Ω Ω Ω + Ψ . . . . . . . . . . . . . . . Ω Ω . . . Ω Ω + Ψ         
.

Let θ = (µ, β, Ω, Ψ, σ 2 ), the vector of all the parameters of the model where β denotes the vector of unit effect β = (β 1 , . . . , β K ). The aim of this paper is to propose an estimation of θ by maximizing the likelihood of the observations y = (y ijk ) ijk .

Let us denote bi := µ + b i . The likelihood of y can be written as:

p(y; θ) = p(y, φ, b; θ)d(φ, b)
where p(y, φ, b; θ) is the likelihood of the complete data (y, φ, b), with φ = (φ ik ) i=1,...,n,k=1,...,K and b = ( b1 , . . . , bn ). Because of the nonlinearity of the regression function f with respect to the random effects φ ik , the likelihood has no closed form. Therefore, the maximization of the likelihood in θ, θ ∈ Θ, is a complex problem. We propose to use a stochastic version of the EM algorithm, which is presented in detail in the next section.

3 Estimation algorithm Let us show that the function Q can be reduced in the case of a MNLMEM. First, let us note that as p( b|y, φ; θ) = p( b|φ; θ), by application of the Bayes theorem we have:

The SAEM algorithm

p(φ, b|y; θ) = p( b|φ; θ)p(φ|y; θ) = n i=1 p( bi |φ i ; θ)p(φ i |y i ; θ). (2)
Second, through the linearity of the individual parameters model in equation ( 1), the posterior distribution p( bi |φ i ; θ) of the ith subject is explicit: p( bi |φ i ; θ) is a Gaussian distribution N (m(φ i , θ), V (θ)) of mean and variance equal to:

m(φ i , θ) = V (θ) Ψ -1 K k=1 (φ ik -β k ) + Ω -1 µ , (3)
V (θ) = (Ω -1 + KΨ -1 ) -1 .
Due to the factorization given in equation ( 2), function Q can be rewritten as: Let us detail the explicit computation of the first integral, denoted by R(y, φ, θ, θ ′ )

Q(θ|θ ′ ) = L c (y, φ, b; θ)p( b|φ; θ ′ )d b p(φ|y; θ ′ )dφ.
R(y, φ, θ, θ ′ ) = L c (y, φ, b; θ)p( b|φ; θ ′ )d b.
This integral has an analytical form. Indeed, the complete log likelihood

L c (y, φ, b; θ) is equal to L c (y, φ, b; θ) = - 1 2 K k=1 n i=1 nik j=1 log(2πσ 2 g 2 (t ijk , φ ik )) - 1 2 K k=1 n i=1 nik j=1 (y ijk -f (t ijk , φ ik )) 2 σ 2 g 2 (t ijk , φ ik ) - nK 2 log(2π det Ψ) - 1 2 K k=1 n i=1 (φ ik -bi -β k ) t Ψ -1 (φ ik -bi -β k ) - n 2 log(2π det Ω) - 1 2 n i=1 ( bi -µ) t Ω -1 ( bi -µ).
As the posterior distribution p( b|φ; θ) is known (equation 3), R(y, φ, θ, θ ′ ) is equal to

R(y, φ, θ, θ ′ ) = - 1 2 i,j,k log(2πσ 2 g 2 (t ijk , φ ik )) 1 2 i,j,k (y ijk -f (t ijk , φ ik )) 2 σ 2 g 2 (t ijk , φ ik ) (4)
- nK 2 log(2π det Ψ) - nK 2 Ψ -1/2 V (θ ′ )Ψ -1/2 - 1 2 i,k (φ ik -m(φ i , θ ′ ) -β k ) t Ψ -1 (φ ik -m(φ i , θ ′ ) -β k ) - n 2 log(2π det Ω) - n 2 Ω -1/2 V (θ ′ )Ω -1/2 - 1 2 n i=1 (m(φ i , θ ′ ) -µ) t Ω -1 (m(φ i , θ ′ ) -µ).
Therefore Q is reduced to the computation of the second integral under the posterior distribution p(φ|y; θ) as follows:

Q(θ|θ ′ ) = R(y, φ, θ, θ ′ )p(φ|y; θ ′ )dφ. (5)
Given the non-linearity of function f with respect to φ, the posterior distribution p(φ|y; θ ′ ) has no closed form and the function Q defined by ( 5 Let us detail this SAEM algorithm in the case of two-levels NLMEMs. Let us note that the quantity R(y, φ, θ, θ ′ ) belongs to the regular curved exponential family, i.e, it can be written as

R(y, φ, θ, θ ′ ) = -Λ(θ) + S(y, φ, θ ′ ), Φ(θ) , (6)
where ., . is the scalar product, Λ and Φ are two functions twice continuously differentiable on Θ and S(y, φ, θ ′ ) is known as the minimal sufficient statistics of the complete model. Those statistics are detailed later. In this case, the Q function is reduced to

Q(θ|θ ′ ) = -Λ(θ) + S(y, φ, θ ′ )p(φ|y; θ ′ )dφ , Φ(θ) ,
In this case, at the ℓth iteration, the SAEM algorithm proceeds as follows:

• Simulation step: simulation of the missing data (φ

i ) i under the conditional distribution p(φ|y; θ ℓ )

• Stochastic approximation step: computation of a stochastic approximation s ℓ+1 of E S(y, φ, θ ℓ )|y; θ ℓ = S(y, φ, θ ℓ )p(φ|y; θ ℓ )dφ, using (γ ℓ ) ℓ≥0 , a sequence of positive numbers decreasing to 0:

s ℓ+1 = s ℓ + γ ℓ (S(y, φ (ℓ) , θ ℓ ) -s ℓ ).
• Maximization step: update of the estimate θ ℓ+1 :

θ ℓ+1 = arg max θ∈Θ (-Λ(θ)) + s ℓ+1 , Φ(θ)) ) .
The sufficient statistics of the complete model ( 4) evaluated during the SA step of the SAEM algorithm are as follows:

s 1,i,ℓ+1 = s 1,i,ℓ + γ ℓ K k=1 φ (ℓ) ik -s 1,i,ℓ , i = 1, . . . , N, s 2,k,ℓ+1 = s 2,k,ℓ + γ ℓ n i=1 φ (ℓ) ik -s 2,k,ℓ , k = 1, . . . , K, s 3,ℓ+1 = s 3,ℓ + γ ℓ n i=1 m(φ (ℓ) i , θ ℓ ) t m(φ (ℓ) i , θ ℓ ) -s 3,ℓ , s 4,ℓ+1 = s 4,ℓ + γ ℓ K k=1 n i=1 φ (ℓ) ik -m(φ (ℓ) i , θ ℓ ) t φ (ℓ) ik -m(φ (ℓ) i , θ ℓ ) -s 4,ℓ , s 5,ℓ+1 = s 5,ℓ + γ ℓ   i,j,k y ijk -f (t ijk , φ (ℓ) ik ) g(t ijk , φ (ℓ) ik ) 2 -s 5,ℓ   ,
The expression of the M step is obtained by derivation of equation ( 4). The parameter estimates are as follows:

µ ℓ+1 = V ( θ ℓ ) Ψ -1 ℓ 1 n n i=1 s 1,i,ℓ+1 - K k=1 β k,ℓ + V ( θ ℓ ) Ω -1 ℓ µ ℓ , β k,ℓ+1 = s 2,k,ℓ+1 n -µ ℓ+1 , for k = 2, . . . , K, Ω ℓ+1 = V ( θ ℓ ) + s 3,ℓ+1 n -( µ ℓ+1 ) t µ ℓ+1 , Ψ ℓ+1 = V ( θ ℓ ) + s 4,ℓ+1 nK - 1 K K k=1 ( β k,ℓ+1 ) t β k,ℓ+1 , σ 2 ℓ+1 = s 5,ℓ+1 n i=1 K k=1 n ik .
Comparing with the classic SAEM algorithm for single-level NLMEMs, the extension of SAEM to the two-levels NLMEMs is finally split into an explicit EM algorithm and a stochastic EM part. Furthermore, it requires the computation of two intermediate quantities (the conditional expectations m(φ i , θ) and variance V (θ) of the between-subject random effects parameters b i ) as well as two additional sufficient statistics (S 3 and S 4 ), functions of m(φ i , θ). The M-step differs from the one of the classic SAEM for single-level NLMEMs, especially for the estimation of the variance matrix Ω and Ψ which uses the additional quantity V (θ).

As proved by (15; 17), the convergence of the SAEM algorithm is ensured under the following assumption:

Assumption (A1):

1. Functions Λ and Φ are twice continuously differentiable on Θ.

2. The log-likelihood log p(y; θ) is d times differentiable on Θ, where d is the dimension of S(y, φ, θ ′ ).

Function s defined as

s(θ, θ ′ ) = S(y, φ, θ ′ )p(φ|y; θ)dφ
is continuously differentiable on Θ with respect to its first variable.

4. For all ℓ in N,

γ ℓ ∈ [0, 1], ∞ ℓ=1 γ ℓ = ∞ and ∞ ℓ=1 γ 2 ℓ < ∞.
For a convenient step sizes sequence γ ℓ , the assumption (A1) is trivially checked in our model.

A choice of step sizes sequence γ ℓ is presented in Section 4.

However, the simulation step of the SAEM algorithm, which performs the simulation of the non-observed vector φ under the posterior distribution p(φ|y; θ) cannot be directly performed because the posterior distribution is only known up to a constant. In this case, Kuhn and Lavielle [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] propose to combine the SAEM algorithm with a Markov Chain Monte Carlo (MCMC) procedure for the simulation step. This version of the SAEM-MCMC algorithm can be used for the estimation of MNLMEMs. The MCMC procedure used in this case is detailed in section 3.2.

As proved by [START_REF] Kuhn | Coupling a stochastic approximation version of EM with a MCMC procedure[END_REF], the convergence of the SAEM-MCMC algorithm is ensured under assumption (A1) and the following additional assumption:

Assumption (A2):
For any θ in Θ, we assume that the conditional distribution p(.|y; θ) is the unique limiting distribution of a transistion probability Π θ , that has the following properties:

1. For any compact subset V of Θ, there exists a real constant L such that for any (θ, θ ′ ) in

V 2 sup {φ,φ ′ }∈E |Π θ (φ ′ |φ) -Π θ ′ (φ ′ |φ)| ≤ L θ -θ ′ .
2. The transition probability Π θ supplies an uniformly ergodic chain whose invariant probability is the conditional distribution p(φ|y; θ), i.e.

∃K θ ∈ R + , ∃ρ θ ∈]0, 1[ | ∀ℓ ∈ N Π ℓ θ (•|φ) -p(•, •|y; θ) T V ≤ C θ ρ ℓ θ ,
where • T V is the total variation norm. Furthermore,

C = sup θ∈Θ C θ < ∞ and ρ = sup θ∈Θ ρ θ < 1.
3. Function S is bound on E.

At iteration ℓ, the S-step of the SAEM-MCMC algorithm consists thus in simulating φ (ℓ) with the transistion probability Π θℓ φ (ℓ-1) |dφ (ℓ) .

The assumption (A2.2) is the most delicate to check, and it depends on the choice of the MCMC algorithm. This is detailed in the next subsection after presenting the MCMC procedure.

In practice, the SAEM algorithm being a stochastic algorithm, there exists no deterministic convergence criterion which could be used to stop the iterations of the algorithm as soon as the convergence is reached. Therefore, we recommend to implement the SAEM algorithm with a sufficiently large number of iterations and to graphically check the convergence by plotting the values of the SAEM estimates obtained along iterations versus the iterations. Such a figure is described in Section 4.

MCMC algorithm for the simulation step

Let us detail the simulation step of the SAEM-MCMC algorithm, which performs the simulation of the missing data φ through a Markov chain which has p(φ|y; θ) as unique stationary distribution. For subject i, by Bayes formula, this conditional distribution is proportional to

p(φ i |y i ; θ) ∝ K k=1 p(y ik |φ ik ; θ)p(φ i ; θ).
We propose to use a Metropolis-Hastings (M-H) algorithm to simulate this Markov chain.

Let us recall the principle of this algorithm. At iteration r of the M-H algorithm, given the current value φ (r) i of the Markov Chain, the M-H algorithm proceeds as follows:

1. Simulate φ c i with a proposal distribution q(•, φ

i )

2. Compute the acceptance probability

α(φ c i , φ (r) i ) = K k=1 p(y ik |φ c ik ; θ)p(φ c i ; θ) K k=1 p(y ik |φ (r) ik ; θ)p(φ (r) i ; θ) q(φ c i , φ (r) i ) q(φ (r) i , φ c i ) 11
3. Simulate u with a uniform distribution U[0, 1]

4. Update the Markov chain with

φ (r+1) i =      φ c i if u ≤ α(φ c i , φ (r) i ) φ (r) i else
The convergence of the M-H algorithm strongly depends on the choice of the proposal distribution q. The convergence is ensured for some proposal distributions such as independent (q(•, φ [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: a general framework[END_REF]. These proposals are detailed below. Given the dimension of φ, we also consider a Metropolis-Hastings-Within-Gibbs algorithm, combining both Gibbs algorithm and M-H procedure. The advantage of the Gibbs algorithm is to reduce the multi-dimensional simulation problem to the successive simulations of one-dimension vectors. Finally, at iteration ℓ of the SAEM algorithm, given the current estimate θ ℓ , we combine the three following proposal transitions:

(r) i ) independent of φ (r) i ) or symmetrical (q(•, φ (r) i ) = q(φ (r) i , •)) proposals
1. the prior distribution of φ i , that is the Gaussian distribution N ( µ ℓ + β ℓ , Γ ℓ ), corresponding to an independent M-H algorithm, 2. the multidimensional random walks N (φ where Γ ℓ is equal to

(ℓ-1) i , ρ Γ ℓ ) (symmetric proposal),
Γ ℓ =          Ω ℓ + Ψ ℓ Ω ℓ . . . Ω ℓ Ω ℓ Ω ℓ + Ψ ℓ . . . . . . . . . . . . . . . Ω ℓ Ω ℓ . . . Ω ℓ Ω ℓ + Ψ ℓ         
.

Given the proposal distributions, as previously detailed, and using the theoretical convergence results proposed in [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: a general framework[END_REF], this hybrid Gibbs algorithm converges and generates an uniformly ergodic chain with p(φ|y; θ) as the stationary distribution. Consequently, by applying the convergence theorem of Kuhn and Lavielle ( 16) and under assumptions (A1) and (A2), we prove that the estimate sequence ( θ ℓ ) ℓ≥0 produced by the extended SAEM algorithm converges towards a (local) maximum of the likelihood p(y; θ).

In practice, the convergence of the MCMC algorithm is difficult to verify. As in Bayesian inference, the only convergence criteria existing for MCMC procedure are graphical criteria. We have to check if the estimate sequence explores a sufficiency large domain of the Markov chain.

A convergence figure is presented and commented in section 4.

Estimation of the Fisher information matrix and the likelihood

i ; θ) = 1 T T t=1 p(y i , φ (t) i ; θ) q i (φ (t) i )
where (φ (t) i ) t=1,...,T are simulated using the individual instrumental distribution q i . As an individual instrumental distribution q i , we propose a Gaussian approximation of the individual conditional posterior distribution p(φ i |y i ; θ).

4 Simulation study: a PK example

Simulation settings

The objective of this simulation study is to illustrate the main statistical properties of the extended SAEM algorithm (bias, root mean square errors, group comparison tests) and to compare them to the FOCE algorithm. We do not use the AGQ algorithm, since procedure NLMIXED does not succeed, in practice, in estimating such complex variance models.

We use the PK data of orally administered theophyllin to define the population model for the simulation study. These data are classical ones in population pharmacokinetics, often used for software evaluation [START_REF] Pinheiro | Mixed-effect models in S and Splus[END_REF]. We assume that concentrations can be described by a one-compartment model with first order absorption and first order elimination:

f (t, φ) = DK a V K a -Cl e -Cl V t -e -Kat
where D is the dose, V is the volume of distribution, K a is the absorption rate constant and Cl is the clearance of the drug elimination. These parameters are positive and distributed according to a log-normal distribution. Thus, φ has the following components: φ = (log V, log K a , log AU C),

with AU C = D/Cl. We assume identical sampling times for all subjects: for all i in 1, . . . , n, k = 1, 2, t ijk = t j for j = 1, . . . , J. Additive Gaussian random effects are assumed for each parameter with a diagonal covariance matrix Ω and a a diagonal covariance matrix Ψ. Let

ω 2 = (ω 2 V , ω 2 Ka , ω 2 AU C ) and ψ 2 = (ψ 2 V , ψ 2 Ka , ψ 2 AU C
) denote the vector of the variances of the random effects. A combined error model is assumed by setting g(t, φ) = 1 + f (t, φ).

We set the dose for all subjects to the value of 4 mg. For all the parameters, the values are those proposed by Panhard and Mentré (22): log V = -0.73, log K a = 0.39 and log AU C = 4.61,

ω 2 V = 0.01, ω 2 Ka = 0.04, ω 2 AU C = 0.04, ψ 2 V = 0.0025, ψ 2 Ka = 0.01, ψ 2
AU C = 0.01 and σ 2 = 0.01. We generate n = 24 and n = 40 total numbers of subjects with J = 10 blood samples per subject, taken at 15 minutes, 30 minutes, 1, 2, 3.5, 5, 7, 9, 12 and 24 hours after dosing. The individual data of one simulated trial are displayed in Figure 1.

[Figure 1 about here.]

Evaluation of estimates

Our aim is to evaluate and compare the estimates produced by the extended SAEM algorithm with those produced by the nlme function of the R software. We fit the simulation model and compute the relative bias and relative root mean square error (RMSEs) for each component of θ from 1000 replications of the two trials described below (n = 24 and n = 40 total number of subjects).

The simulation model includes a treatment effect on all components of θ. We test the null hypothesis {β log AU C = 0} using the Wald test. We also fit the model where the treatment effect on log AU C is not estimated, and test the same null hypothesis using the Likelihood Ratio Test (LRT).

The SAEM algorithm is implemented with 500 iterations. During the first 200 iterations, a constant step size γ ℓ = 1 is chosen, in order to let the Markov chain explore the parameters domain. Then during the last 300 iterations, the stochastic approximation scheme is implemented with a step size equal to γ ℓ = 1 ℓ-200 at iteration ℓ. This choice of size sequence verifies convergence assumption (A1.1). The evolution of each SAEM parameter estimates is plotted against iterations (logarithmic scale) on Figure 2. During the first iterations of the SAEM algorithm, the estimate sequences explore randomly some neighborhoods of the initial values, through the Markov chain simulation. In particular, these behaviors are clearly visible for the fixed effect parameters (µ and β). After 200 iterations, the estimates converge then rapidly to a neighborhood of the maximum likelihood, due to the stochastic approximation scheme. In this example, the iteration number has been chosen such that the convergence is clearly attained before the last iteration.

[Figure 2 about here.]

The relative bias and RMSEs obtained on the 1000 datasets with n = 24 and n = 40 subjects are presented in Table 1.

[Table 1

about here.]

The bias and the RMSEs of the fixed effects (µ) are small with the SAEM algorithm and are almost half of those obtained with nlme (especially the RMSEs). The bias and RMSEs of the unit effect (β) are small and on the same order with both methods. For the between-subject variability parameters (Ω), the bias are reduced with SAEM, while the RMSEs are of the same order with both methods. For the within-subject variability parameters (Ψ), the bias and the RMSEs are satisfactory, and on the same order with both methods. The bias and RMSE for σ 2 are small and satisfactory for both methods.

The type I error of the Wald test and of the LRT are evaluated on the same 1000 datasets.

For n = 24, the type I errors are 6.0% and 6.5% for SAEM and nlme, respectively, for the Wald test, and 4.6% and 5.6% for SAEM and nlme, respectively, for the LRT. For n = 40, the type I error are 5.6% and 5.4% for SAEM and nlme, respectively, for the Wald test, and 5.8% and 5.2% for SAEM and nlme, respectively, for the LRT.

5 Application to the population pharmacokinetics of atazanavir with tenofovir

Study population: ANRS 107 -Puzzle 2 study

The Puzzle 2 -ANRS 107 trial was a randomized open-label, multiple-dose study supported by the French Agence Nationale de Recherche sur le Sida (ANRS) with HIV-infected patients in treatment failure with their antiretroviral therapy. Patients were randomized to receive for the first two weeks either their unchanged treatment with PIs and nucleoside reverse transcriptase inhibitors (NRTIs) (group 1) or unchanged treatment with NRTIs in combination with atazanavir (300 mg QD) plus ritonavir (100 mg QD) as a substitute for the failing PI therapy (group 2). From week 3 (day 15) to week 26, patients from either group switched to atazanavir (300 mg QD) plus ritonavir (100 mg QD) plus tenofovir DF at 300 mg QD and NRTIs selected according to the baseline reverse transcriptase genotype of the HIV isolated in each patient.

In this paper, we analyze concentration data obtained from 10 patients from group 2 who were included and measured at 1, 2, 3, 5, 8, and 24 h after administering drug during each treatment period. Those exact dosing intervals were recorded. The objective of the substudy was to measure the pharmacokinetic parameters of atazanavir (administered with ritonavir) either before (day 14 [week 2]) or after (day 42 [week 6]) initiation of tenofovir DF in HIVinfected patients in order to detect pharmacokinetic interactions of tenofovir on atazanavir.

Data of this substudy were analyzed using a nonlinear mixed effect model by Panhard et al [START_REF] Panhard | Impact of modelling intra-subject variability on tests based on non-linear mixed-effects models in cross-over pharmacokinetic trials with application to the interaction of tenofovir on atazanavir in hiv patients[END_REF] and a significant effect of the co-administration of tenofovir on the pharmacokinetic parameters of atazanavir was found using the nlme function of the Splus software.

The aim of the present analysis is to evaluate the effect of tenofovir on the PK parameters of atazanavir using the SAEM algorithm and the Wald test described in section 4.

We use the one-compartment model with zero-order absorption proposed by Panhard et al. [START_REF] Panhard | Impact of modelling intra-subject variability on tests based on non-linear mixed-effects models in cross-over pharmacokinetic trials with application to the interaction of tenofovir on atazanavir in hiv patients[END_REF] to describe atazanavir concentrations:

f (t, φ) = F D T a Cl (1 -e -Cl V t ) t<Ta + e -Cl V τ t<Ta (1 -e -Cl V Ta )e -Cl V (t-Ta) (1 -e -Cl V τ )
with F the bioavailability, V the volume of distribution of atazanavir, (T a ) the absorption duration, Cl the elimination clearance and τ the dosing interval (24 hours until the PK visit).

The vector of the logarithm of the identifiable parameters is φ = (log(V /F ), log(T a ), log(AU C)).

Data of both treatment periods are simultaneously analyzed using a NLMEM with two levels of variability (the between-patient and within-patient variabilities) on each PK parameter. A treatment effect is also estimated for each PK parameter, and a homoscedastic error model is used.

Results

Concentrations versus time are displayed in Figure 3. The SAEM algorithm succeeds in the estimation of all the parameters. The resulting parameters estimates are displayed in table 2. The SAEM algorithm estimates the AU C between-patient variability and the V /F and T a within-patient variabilities to 0.48, 0.69 and 0.19, respectively. The three other variances are estimated to 0.

[Table 2

about here.]

A significant effect of co-medication with tenofovir is found on log(AU C) (p=0.00015) with the Wald test based on the SAEM algorithm.

The individual prediction curves for the two periods are overlaid on the concentration data on The main original element of this study is the development of the SAEM algorithm for twolevels non-linear mixed effects models. We extend the SAEM algorithm developed by Kuhn and Lavielle [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF], which was not yet adapted to the case of MNLMEMs with two levels of random effects. This algorithm will be implemented in the 3.1 version of the monolix software, freely available on the following website: http://monolix.org. The two levels of random effects are the between-subject variance and the within-subject (or between-unit) variance, with N subjects and K units, with no restriction on N or K. We show that the SAEM algorithm is split into two parts: an explicit EM algorithm and a stochastic EM part. The integration of the term p(b|φ; θ)

in the likelihood results in the derivation of two additional sufficient statistics compared to the original algorithm. Furthermore it uses two intermediate quantities, the conditional expectations and variance of the between-subject random effects parameters b. The addition of higher levels of variability would therefore require other extensions of the algorithm.

The convergence of the algorithm is monitored from a graphical criterion, as shown in Fig1.

An automatic implementation of that stopping criterion to optimize both the number of iterations and the stochastic approximation step should be considered in future work and extension of the Monolix software.

The simulation study illustrates the accuracy of our approach. We show that the bias and RMSEs obtained by the extended SAEM algorithm are satisfactory for all parameters. The bias are reduced compared to those obtained with the FOCE algorithm implemented in the nlme function of the R software. The bias are especially divided by two for the fixed effects parameters with SAEM. Furthermore, whereas the nlme implementation of the FOCE algorithm does not always converge with both between-and within-patient variability on all parameters, the extended SAEM algorithm does. We develop the tests for a difference between the units, and we obtain type I errors close to the expected 5% for the Wald test and the LRT.

The analysis of the pharmacokinetics of atazanavir with tenofovir in the Puzzle 2 -ANRS 107 trial also demonstrates the ability of the extended SAEM algorithm to detect treatment interaction on a real data set. When testing for an interaction of tenofovir on the PK of atazanavir, the impact of tenofovir on the absorption of atazanavir is confirmed; more precisely, a decrease of the AUC of atazanavir as shown by Panhard et al. [START_REF] Panhard | Impact of modelling intra-subject variability on tests based on non-linear mixed-effects models in cross-over pharmacokinetic trials with application to the interaction of tenofovir on atazanavir in hiv patients[END_REF] is found.

We compare the extension of SAEM to the FOCE algorithm, that is the most popular method in population pharmacokinetics, which is one of the largest application fields of NLMEM. We try to use the NLMIXED procedure of SAS implementing Gaussian quadrature. However, procedure NLMIXED does not succeed, in practice, in estimating such complex variance models, on our simulated data, neither on the atazanavir dataset. The next step is a comparison with a Bayesian estimation of the parameters using Winbugs, which is beyond the scope of this paper.

The next ambitious development would be an extension of the algorithm to the case of MNLMEMs with more than two levels of random effects, in order to analyze, for instance, genetic data where more than one generation of parents are taken into account. However, it would be difficult to develop a general estimation method since it strongly depends on the relation (linear or not) of the different levels of random effects.

To conclude, the extended SAEM algorithm combines the statistical properties of an exact method together with a high computational efficiency. We thus recommend the use of this method in MNLMEMs.

A Index of notations

Model notations

i (i = 1, • • • , n): index of subject j (j = 1, • • • , n ik ): index of
i ): proposal distribution of the Metropolis-Hastings (M-H) algorithm

φ c i : candidate simulated using q(φ c i , φ (ℓ) i ) α(φ c i , φ (ℓ)

The

  EM algorithm introduced by Dempster et al. (11) is a classical approach to estimate parameters of models with non-observed or incomplete data. In two-levels NLMEMs, the non-observed data are the individual parameters (φ, b) and the complete data of the model is (y, φ, b). Let us denote L c (y, φ, b; θ) = log p(y, φ, b; θ) the log-likelihood of the complete data. The principle of the iterative EM algorithm is to maximize the function Q(θ|θ ′ ) = E(L c (y, φ, b; θ)|y; θ ′ ) where the expectation is the conditional expectation under the posterior distribution p(φ, b|y; θ ′ ), the maximization of Q being often easier than the direct maximization of the observed data loglikelihood. Each iteration of the EM algorithm is computed through two steps: the Expectation step (E-step) and the Maximization step (M-step). At the ℓ th iteration of the algorithm, the E step is the evaluation of Q(θ | θ ℓ ), while the M step updates θ ℓ by maximizing Q(θ | θ ℓ ).

  Because of the explicit posterior distribution of random effects b given in equation (3), the computation of this conditional expectation can be split into two parts : the computation of the integral with respect to the posterior distribution of b which has an analytical form, and the computation of the integral with respect to the posterior distribution of φ which has no analytical form. Therefore the EM algorithm is split into an explicit direct EM algorithm for the computation of the first integral and the use of a stochastic version of the EM algorithm for the computation of the second integral.

  ) is intractable. Thus we propose to use the stochastic version SAEM of the EM algorithm proposed by Delyon et al. (15) which evaluates the integral Q by a stochastic approximation procedure.

  where ρ is a scaling value chosen to ensure a satisfactory acceptation rate, namely around 30% as proposed in (19), 3. a succession of Kp unidimensional Gaussian random walks (symmetric proposal), i.e each component of φ i is successively updated, leading to a Metropolis-Hastings-Within-Gibbs algorithm,

Figure 3

 3 Figure 3 for 10 patients. The goodness-of-fit plots (population and individual predicted concentrations versus observed concentrations; standardized residuals versus predicted concentrations and versus time) are judged satisfactory, and are displayed in Figure 4.

  measurement of subject i for unit k k (k = 1, . . . , K): index of unit t ijk : measurement time in unit k for subject i and measurement j y ijk : observation in unit k for subject i at timet ijk y = (y ijk ) ijk : vector of the observations in the K units for all the n subject f and g: non-linear functions linking observations to sampling times φ ik : p-vector of the parameters of subject i for unit k φ i = (φ i1 , . . . , φ iK ): pK-vector of individual parameter of subject i φ = (φ ik ) i=1,...,n,k=1,...,K µ: p-vector of the mean of the individual parameters for k = 1 β k : effect of the kth unit in comparison to this first unit β = (β 1 , . . . , β K ): vector of the unit effects b i : random effect of size p of subject i bi := µ + b i b = ( b1 , . . . , bn ) c ik : random effect of size p of subject i and unit k ε ijk : measurement error Ω: p × p between-subject covariance matrix Ψ: p × p within-subject covariance matrix Γ: pK × pK covariance matrix of the individual parameters φ i (i = 1, • • • , n) σ 2 : variance of the measurement error θ = (µ, β, Ω, Ψ, σ 2 ): vector of all the parameters Algorithm notations p(y, φ, b; θ): likelihood of the complete data L c (y, φ, b; θ): log-likelihood of the complete data p(φ, b|y; θ ′ ): posterior distribution of (φ, b) given (y; θ ′ ) p(φ|y; θ ′ ): posterior distribution of φ given (y; θ ′ ) p( b|φ; θ): posterior distribution of b given (φ; θ) m(φ, θ): mean of p( b|φ; θ) V (θ): variance of p( b|φ; θ) Q(θ|θ ′ ) := E(L c (y, φ, b; θ)|y; θ ′ )) R(y, φ, θ, θ ′ ) := L c (y, φ, b; θ)p( b|φ; θ ′ )d b Λ(θ, θ ′ ) and Φ(θ, θ ′ ): functions of Θ × Θ ℓ: iteration number of the SAEM algorithm φ (ℓ) i : missing data simulated at iteration ℓ S(y, φ): minimal sufficient statistics of the complete model of dimension d s ℓ+1 : stochastic approximation of E S(y, φ)| θ ℓ (γ ℓ ) ℓ≥0 : step sizes sequence Π θ : transition probability q(φ c i , φ

iFigure 1 :Figure 2 :

 12 Figure 1: Simulated theophyllin concentration data for 24 subjects during the first period (plain line) and during the second period (dotted line)

Figure 3 :

 3 Figure 3: Individual concentrations and individual predicted curves for the pharmacokinetics of atazanavir in 10 subjects: x and * , observations with and without tenofovir, respectively; dotted and plain line, individual predictions of the atazanavir pharmacokinetics with and without tenofovir, respectively.

Figure 4 :

 4 Figure 4: Goodness-of-fit plots for atazanavir final population PK model: population (a) and individual (b) predicted concentrations (in ng/mL) versus observed concentrations (in ng/mL), standardized residuals versus predicted concentrations (in ng/mL) (c) and versus time (in hours) (d).

Table 1 :

 1 Relative biases (%) and relative root mean square errors (RMSE) (%) of the estimated parameters evaluated by the extended SAEM algorithm and the FOCE algorithm (nlme function) from 1000 simulated trials.

		n=24 subjects		n=40 subjects
		Bias		RMSE	Bias		RMSE
		SAEM nlme SAEM nlme	SAEM nlme SAEM nlme
	V	0.01	0.53	3.9	6.4	-0.06	0.54	2.91	5.00
	k a	0.48 -1.48	14.4	24.4	0.02 -3.07	10.79	18.4
	AU C	-0.08 -0.20	1.0	1.5	-0.11 -0.24	0.79	1.13
	β V	-0.00 -0.01	3.6	3.6	-0.05 -0.06	2.83	2.80
	β ka	-0.73 -0.76	14.2	18.8	0.24	0.27	10.73 10.60
	β AU C ω 2 V ω 2 ka ω 2 AU C ψ 2 V ψ 2 ka ψ 2 AU C σ 2	0.02 -0.02 -5.13 -5.92 -3.99 -7.07 -4.88 -7.29 -8.67 -7.29 -10.94 -9.09 -5.37 -5.37 -0.33 -0.10	0.7 38.7 42.4 34.5 69.4 73.5 43.6 7.7	1.1 38.4 41.5 34.2 68.5 72.0 42.6 7.7	0.00 -3.45 -4.28 1.14 -3.23 -4.63 -1.51 -3.80 -5.93 -4.91 -7.06 -5.17 -4.92 -5.79 0.28 0.67	0.57 30.30 30.37 1.14 33.49 33.32 27.41 27.02 58.78 57.81 62.00 60.60 33.31 32.47 6.03 6.09

Table 2 :

 2 Pharmacokinetic parameters of atazanavir (estimate and SE (%)) estimated with the SAEM

	Parameters	Estimate SE (%)
	log(V /F ) (L)	4.01	5.79
	log(T a ) (h)	1.36	6.72
	log(AU C) (ng.mL -1 .h)	10.67	1.61
	β V /F	0.12	267.43
	β Ta	0.33	45.03
	β AU C	-0.38	25.31
	ω V /F	0	-
	ω Ta	0	-
	ω AU C	0.48	25.48
	ψ V /F	0.55	28.30
	ψ Ta	0.16	35.76
	ψ AU C	0	-
	σ (ng.mL -1 )	732.29	8.40

Acknowledgments

The authors are grateful to Marc Lavielle for his supportive advice and help. The authors would like to thank the principal investigator (Dr Christophe Piketti), the pharmacology coordinator (Dr Anne-Marie Taburet) and the scientific committee of the Puzzle 2 -ANRS 107 trial for giving us access to the concentration measurements of the patients.

During this work, Xavière Panhard was supported by a grant from the Agence Nationale de la Recherche (France).


Omer-Landry Nguena-Timo

Pierre-Alain Reynier
email: pierre-alain.reynier@lif.univ-mrs.fr

On Characteristic Formulae for Event-Recording Automata

A standard bridge between automata theory and logic is provided by the notion of characteristic formula. This paper investigates this problem for the class of event-recording automata (ERA), a subclass of timed automata in which clocks are associated with actions and that enjoys very good closure properties (complementation, determinization...). We first study the problem of expressing characteristic formulae for ERA in Event-Recording Logic (ERL), a logic introduced by Sorea to express event-based timed specifications. We prove that the construction proposed by Sorea for ERA without invariants is false. More generally, we prove that bisimulation can not be expressed in ERL for the class of ERA, even without invariants. Then, we introduce the logic WTµ, a new logic for event-based timed specifications, closer to the timed logic Lν . We prove that it is strictly more expressive than ERL, and that its model-checking problem against ERA is EXPTIME-complete. Finally, we provide constructions for characterizing ERA up to timed (bi)similarity and study the complexity issues.

Introduction

In the untimed setting, automata and logics are central tools for the formal verification of reactive systems. While the system is usually modelled as an automaton, the specification may be described both as a formula of a logic or as an automaton. In the first case the correctness of the system reduces to a model checking problem, whereas in the second case it requires to compare the two automata, and different relations can be envisaged, such as bisimulation [START_REF] Sangiorgi | Bisimulation: From the origins to today[END_REF] or language inclusion. A standard bridge between automata theory and logic is provided by the notion of characteristic formula. A characteristic formula is a formula in a temporal logic that completely characterizes the behaviour of an automaton modulo some choosen relation. For the class of timed automata [START_REF] Alur | A theory of timed automata[END_REF], a solution has first been proposed in [START_REF] Laroussinie | From timed automata to logic -and back[END_REF], providing formulae in the logic L ν . Then, these results have been improved in [START_REF] Aceto | Characteristic formulae for timed automata[END_REF], yielding linear constructions.

The subclass of Event-Recording Automata [START_REF] Alur | A determinizable class of timed automata[END_REF] (ERA) is obtained by restricting clocks to be associated with events. This class enjoys good closure properties such as determinization and complementation. It has thus attracted attention to characterize its expressive power in terms of some timed logic [START_REF] Raskin | The logic of event clocks -decidability, complexity and expressiveness[END_REF][START_REF] Souza | A logical characterisation of event clock automata[END_REF], but logics considered there are linear-time. This paper investigates the problem of constructing characteristic formulae for the class of event-recording automata, up to timed similarity and timed bisimilarity, using a branching-time logic devoted to event-based timed specifications. Such a logic, called Event-Recording Logic (ERL) and introduced by Sorea in [START_REF] Sorea | A decidable fixpoint logic for time-outs[END_REF], extends the mu-calculus by allowing the use of event-clocks.

After recalling standard definitions in Section 2, we first study in Section 3 the problem of expressing characteristic formula for timed bisimulation for ERA in the logic ERL. We prove that an existing attempt, which can be found in [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF], is not correct. We also show that this logic can not characterize ERA up to timed bisimilarity, even when we only consider the subclass of ERA without invariants. Then, we consider in Section 4 a new timed logic, called WT µ [START_REF] Nguena-Timo | The logic WTµ[END_REF], to express the characteristic formulae. The definition of this logic is closer from the definition of L ν as it separates quantifications over discrete successors and time successors. We prove that it is indeed strictly more expressive than ERL, and that its model-checking problem over ERA is EXPTIME-complete. Finally, we provide formulae constructions in WT µ for timed (bi)similarity together with complexity issues in Section 5.

Preliminaries

Let Σ be a finite alphabet, Σ * is the set of finite words over Σ. The sets N, Q, Q ≥0 , R and R ≥0 are respectively the sets of natural, rational, non-negative rational, real and non-negative real numbers. Given a real number x, ⌊x⌋ (resp.

x ) denotes its integral part (resp. its fractionnal part). We consider as time domain T the set Q ≥0 or the set R ≥0 . We consider a finite set X of variables, called clocks. A clock valuation over X is a mapping v : X → T that assigns to each clock a time value. The set of all clock valuations over X is denoted T X . Let t ∈ T, the valuation v + t is defined by (v + t)(x) = v(x) + t, ∀x ∈ X . For a clock y ∈ X , we denote by v[y := 0] the valuation such that for each clock x ∈ X , (v[y := 0])(x) = 0 if x = y, and (v[y := 0])(x)= v(x) otherwise. Finally, 0 denotes the valuation mapping every clock on 0.

In the context of event-recording automata, each clock refers to a specific action. Then, we associate clocks with letters of an alphabet. Given an alphabet Σ, we then denote by X Σ the set of clocks {x a | a ∈ Σ}. We may also write T Σ to represent the set of clock valuations T XΣ . Given a set of clocks X Σ , we introduce two sets of clock constraints over X Σ . The most general one, denoted by C(Σ), is defined by the grammar "g ::

= x ∼ c | x -y ∼ c | g ∧ g | tt" where x, y ∈ X Σ , c ∈ Q ≥0 , ∼ ∈ {<, ≤, =, ≥, >}
and tt stands for true. We also use the proper subset C up (Σ) of upper bounds constraints consisting only of conjunctions of constraints of the form x ≺ c with ≺∈ {<, ≤}. We write v |= g when the clock valuation v satisfies the clock constraint g and denote by g the set of clock valuations v such that v |= g holds.

Timed Transition Systems

Timed transition systems describe systems which combine discrete and continuous evolutions. They are used to define the behavior of timed systems such as Timed Automata [START_REF] Alur | A theory of timed automata[END_REF], or Event-Clock Automata [START_REF] Alur | A determinizable class of timed automata[END_REF] (see below).

Definition 1 (Timed Transition System (TTS)). A timed transition system over the alphabet Σ is a transition system S = Q, q 0 , Σ, → , where Q is the set of states, q 0 ∈ Q is the initial state, and the transition relation

→⊆ Q × (Σ ∪ T) × Q consists of continuous transitions q d -→ q ′ (d ∈ T), and discrete transitions q a -→ q ′ (a ∈ Σ).
Moreover, we require the following standard properties for TTS :

-Time-Determinism : if q d -→ q ′ and q d -→ q ′′ with d ∈ T, then q ′ = q ′′ , -0-Delay : q 0 -→ q, -Additivity : if q d -→ q ′ and q ′ d ′ -→ q ′′ with d, d ′ ∈ T, then q d+d ′ ---→ q ′′ , -Continuity : if q d -→ q ′ , then for every d ′ and d ′′ in T such that d = d ′ + d ′′ ,
there exists q ′′ such that q

d ′ -→ q ′′ d ′′ -→ q ′ .
With these properties, a run of S can be defined as a finite sequence of moves

ρ = q 0 d0 -→ q ′ 0 a0 -→ q 1 d1 -→ q ′ 1 a1
-→ q 2 . . . an --→ q n+1 where discrete and continuous transitions alternate. To such a run corresponds the timed word w = (a i , τ i ) 0≤i≤n over Σ where τ i = i j=0 d j is the time at which a i happens, and we say that the timed word w is accepted by S. The language of S, denoted L(S), is defined as the set of timed words that are accepted by S.

Event-Recording Automata

We consider the restriction of Event-Clock Automata to Event-Recording Automata.

Definition 2 (Event-Recording Automata (ERA) [START_REF] Alur | A determinizable class of timed automata[END_REF]). An event-recording automaton over the alphabet Σ is a tuple A = L, ℓ 0 , Σ, T, I where:

-L is a finite set of locations, -ℓ 0 ∈ L is the initial location, -T ⊆ L × C(Σ) × Σ × L is a finite set of transitions, -I : L → C up (Σ)
associates an upper bound constraint with each location.

We say that an ERA is without invariants if the mapping I associates tt to each location. In this case we may remove component I of the definition of A.

Without loss of generality, we assume that the clock constraints of transitions are consistent with invariants. More formally, we have, for any v ∈ T Σ :

∀(ℓ, g, a, ℓ ′ ) ∈ T, v |= g ⇒ (v |= I(ℓ)) ∧ (v[x a := 0] |= I(ℓ ′ ))
The semantics of an event-recording automaton A is defined in the terms of a timed transition system. Intuitively, it manipulates exactly one clock per action, which allows to measure time elapsed since the last occurrence of this action. The formal definition is given by: Definition 3 (Semantics of an ERA). Given an ERA A = L, ℓ 0 , Σ, T, I , its semantics is given by the TTS S A defined by S A = Q, q 0 , Σ, → where

Q = {(ℓ, v) ∈ L × T Σ | v |= I(ℓ)}, q 0 =
(ℓ 0 , 0), and → consists of continous and discrete transitions:

∀(ℓ, v) ∈ Q, Time-elapsing steps: ∀d ∈ T, we have (ℓ, v) d -→ (ℓ, v + d) iff v + d |= I(ℓ), Discrete steps: ∀a ∈ Σ, we have (ℓ, v) a -→ (ℓ ′ , v ′ ) iff there exists a transition t = (ℓ, g, a, ℓ ′ ) ∈ T such that v |= g and v ′ = v[x a := 0].
Finally, we simply denote by L(A) the language of timed words L(S A ).

We say that an ERA is determinisitic whenever, for every location ℓ ∈ L, letter a ∈ Σ and valuation v ∈ T Σ , there exists at most one transition (ℓ, g, a, ℓ ′ ) ∈ T such that v |= g holds.

As it has been introduced for timed automata in [START_REF] Alur | A theory of timed automata[END_REF], a time-abstract bisimulation based on the construction of regions can also be defined for ERA. We briefly recall here this construction, and refer the reader to [START_REF] Alur | A determinizable class of timed automata[END_REF] for more details.

Definition 4 (Clock Region). We consider a constant K ∈ N. A clock region is an equivalence class of the relation ≃ K over clock valuations. For two valuations v, v ′ ∈ T Σ , we have v ≃ K v ′ iff the following conditions hold:

1. ∀x ∈ X Σ , if v(x) ≤ K or v ′ (x) ≤ K, then ⌊v(x)⌋ = ⌊v ′ (x)⌋, 2. ∀x ∈ X Σ s.t. v(x) ≤ K, then v(x) = 0 ⇐⇒ v ′ (x) = 0, 3. ∀x, y ∈ X Σ s.t. |v(x) -v(y)| ≤ K, then v(x) ≤ v(y) ⇐⇒ v ′ (x) ≤ v ′ (y) .
We let R K (Σ) be the set of clock regions for constant K. We recall that the size of R K (Σ) is in 2 O(m. log Km) where m = |Σ| (see [START_REF] Alur | A determinizable class of timed automata[END_REF]). When the constant K is clear from the context, we denote by [v] the clock region that contains v, and by r the set of clock valuations whose clock region is equal to r. To define the region automaton of an ERA A, we can assume that all the constants occuring in its clock constaints are natural numbers (otherwise, all constants need to be multiplied by the least common multiple of the denominators of all rational numbers appearing in clock constraints).

Definition 5 (Region Automaton). Given an ERA A = L, ℓ 0 , Σ, T, I with integral constants. Let K be some positive integer. We define the region automaton of A for constant K, denoted by

R K (A) = R K (A), Σ ∪ {τ }, → , as follows 3 : -R K (A) = {(ℓ, r) ∈ L × R K (Σ) | ∃v ∈ r s.t. v |= I(ℓ)} 3 τ is an action not in Σ intended to represent time elapsing. -(ℓ, r) τ -→ (ℓ, r ′ ) ⇐⇒ ∃δ ∈ T s.t. (ℓ, v) δ -→ (ℓ, v ′ ) in S A , r = [v] and r ′ = [v ′ ] -∀a ∈ Σ, (ℓ, r) a -→ (ℓ, r ′ ) ⇐⇒ ∃(ℓ, v) a -→ (ℓ, v ′ ) in S A s.t. r = [v] and r ′ = [v ′ ]
It is well known that if K is larger than the largest integer constant that appears in the clock constraints of A, then R K (A) is a time-abstract bisimulation of S A .

Event-Recording Logic

Definition 6 (Event-Recording Logic (ERL) [START_REF] Sorea | A decidable fixpoint logic for time-outs[END_REF]). Let Σ be a finite alphabet, Var be a finite set of variables, the formulae of the Event-Recording Logic over Σ and Var are defined by the grammar:

ϕ ::= tt | ff | X | ϕ 1 ∧ ϕ 2 | ϕ 1 ∨ ϕ 2 | [g, a]ϕ | g, a ϕ | µX.ϕ | νX.ϕ
where g ∈ C(Σ), a ∈ Σ and X ∈ Var.

In the timed logic L ν [START_REF] Laroussinie | From timed automata to logic -and back[END_REF], the formulae have their own clocks and the semantics is then defined using a valuation for the clocks of the formula. When defining the semantics of ERL formulae over some alphabet Σ, the clock constraints range over event clocks associated with Σ. Then, the semantics is defined for TTS corresponding to ERA over the same alphabet Σ, and the clock constraints are evaluated over the valuation of the ERA. Moreover, variables of ERL formulae are dealt with using assignment functions. Formally, an assignment function of variables Var over the set Q is a function V : Var → P(Q). The updating notation V[X := Q ′ ] denotes the assignment V ′ that agrees with V on all variables except X, where

V ′ (X) = Q ′ ⊆ Q.
Definition 7 (Semantics of ERL). Let Σ be a finite alphabet, Var be a finite set of variables, A = L, ℓ 0 , Σ, T, I be an ERA4 over Σ and S A = Q, q 0 , Σ, → be its associated TTS. Consider a formula ϕ ∈ ERL over Σ and Var and an assignment function V of Var over Q. The semantics of ϕ for A under V, denoted ϕ A V , is given by the set of states (ℓ, v) ∈ Q for which the formula holds, and is defined inductively as follows:

tt A V := Q ff A V := ∅ X A V := V(X) ϕ 1 ∧ ϕ 2 A V := ϕ 1 A V ∩ ϕ 2 A V ϕ 1 ∨ ϕ 2 A V := ϕ 1 A V ∪ ϕ 2 A V [g, a]ϕ A V := {(ℓ, v) ∈ Q | ∀δ ∈ T, ∀(ℓ, g ′ , a, ℓ ′ ) ∈ T, v + δ |= g ∧ g ′ ⇒ (ℓ ′ , v ′ ) ∈ ϕ A V , where v ′ = (v + δ)[x a := 0]} g, a ϕ A V := {(ℓ, v) ∈ Q | ∃δ ∈ T, ∃(ℓ, g ′ , a, ℓ ′ ) ∈ T s.t. v + δ |= g ∧ g ′ and (ℓ ′ , v ′ ) ∈ ϕ A V , where v ′ = (v + δ)[x a := 0]} µX.ϕ A V := ∩{Q ′ ⊆ Q | ϕ A V[X:=Q ′ ] ⊆ Q ′ } νX.ϕ A V := ∪{Q ′ ⊆ Q | Q ′ ⊆ ϕ A V[X:=Q ′ ] }
Using standard definitions, we say that a variable X is bounded (resp. free) in a formula ϕ whenever it is (resp. it is not) under the scope of a fix-point operator µ or ν. It is easy to verify that if all variables are bounded in a formula ϕ (we say that ϕ is a sentence), then the semantics of ϕ does not depend on the assignment function. In this case, we omit the subscript V, and given an ERA A, and a configuration q of A, for a sentence ϕ, we write A, q |= ϕ whenever we have q ∈ ϕ A . We also use the shortcut A |= ϕ whenever A, q A 0 |= ϕ. Moreover, we say that a bounded variable X is guarded if it is in the scope of an operator • or [•]. According to [START_REF] Sorea | A decidable fixpoint logic for time-outs[END_REF], one can assume that every bounded variable is guarded.

Remark 1 (On greatest fixpoints). To express characteristic formulae, we shall see later that we need greatest fixpoints on systems of inequations. In this case, we will use a slightly different presentation. Given a finit set Var of variables, we will associate to each variable X a formula D(X) over the variables Var. D is then called a declaration, and the semantics associated with this definition is the largest solution of the system of inequations X ⊆ D(X) for any X ∈ Var. It can be proven (see [START_REF] Bekic | Definable operation in general algebras, and the theory of automata and flowcharts[END_REF] or [START_REF] Cleaveland | A linear-time model-checking algorithm for the alternation-free modal mu-calculus[END_REF]) that this presentation is equivalent. To specify the declaration used, we will add it as subscript of the satisfaction relation |=, writing A, q |= D X.

Timed Behavioral Relations

We now recall the standard definitions of timed simulation and timed bisimulation. These definitions are given for TTS and can thus be used for ERA.

Definition 8 (Timed simulation and timed bisimulation). Consider two

TTS S 1 = Q 1 , q 1 0 , Σ, → 1 and S 2 = Q 2 , q 2 0 , Σ, → 2 . A timed simulation be- tween S 1 and S 2 is a relation R ⊆ Q 1 × Q 2 such that whenever q 1 Rq 2 and α ∈ Σ ∪ T, then: -If q 1 α -→ q ′ 1 then there exists q ′ 2 ∈ Q 2 such that q 2 α -→ q ′ 2 and q ′ 1 Rq ′ 2 .
A relation R is a timed bisimulation between S 1 and S 2 iff the relations R and R -1 are timed simulations.

For states q 1 , q 2 , we write q 1 ≺ q 2 (resp. q 1 ∼ q 2 ) if and only if there exists a timed simulation (resp. a timed bisimulation) R with q 1 Rq 2 .

Finally, we say that a TTS S 2 simulates a TTS S 1 (resp. S 1 and S 2 are bisimilar ) whenever there exists a timed simulation (resp. a timed bisimulation) between S 1 and S 2 such that the pair (q 1 0 , q 2 0 ) of their initial states belongs to the relation R, and then we write S 1 ≺ S 2 (resp. S 1 ∼ S 2 ). We naturally extend these notations to ERA.

Definition 9 (Characteristic formulae).

Let A be an ERA. We say that a sentence ϕ ∈ ERL is a characteristic formula for A if and only if, according to the behavioural relation considered, the following equivalence holds:

Simulation: ∀B ∈ ERA, A ≺ B ⇐⇒ B |= ϕ Bisimulation: ∀B ∈ ERA, A ∼ B ⇐⇒ B |= ϕ
The following standard result relates simulation with language inclusion. Proposition 1. Let A 1 and B 2 be two ERA, we have the following implications:

(i) if A 1 ≺ A 2 , then L(A 1 ) ⊆ L(A 2 ), (ii) if A 2 is deterministic and L(A 1 ) ⊆ L(A 2 ), then A 1 ≺ A 2 .
3 On the use of ERL for characterizing bisimulation

As the logic ERL has been introduced to describe behaviours related to events, it is natural to try to write in this logic characteristic formulae for timed bisimulation for ERA. An attempt can be found in Sorea's thesis [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF] for the class of ERA without invariants. We will first show in this section that this attempt is erroneous, by providing two counter-examples to illustrate how the construction fails. Then, we will prove that it is in fact not possible to express timed bisimilarity for ERA (even without invariants) in the logic ERL.

3.1 On the construction proposed in [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF] In [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF], the author addresses the problem of constructing characteristic bisimulation formulae for ERA without invariants using ERL formulae with greatest fixpoints. We recall here the proposed construction and explain why it fails. Before presenting the construction, we introduce some additional notations. Given an ERA without invariants A = L, ℓ 0 , Σ, T , a location ℓ ∈ L and a letter a ∈ Σ, we define:

-the set of a-labelled transitions leaving ℓ:

Out(ℓ, a) = {t = (ℓ, g, a, ℓ ′ ) ∈ T } -the union of clock constraints of a-labelled transitions leaving ℓ:

En(ℓ, a) = {g | ∃(ℓ, g, a, ℓ ′ ) ∈ Out(ℓ, a)} -the set of locations reached by an a from location ℓ:

F(ℓ, a) = {ℓ ′ | ∃(ℓ, g, a, ℓ ′ ) ∈ Out(ℓ, a)}
The formulae defined in [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF] are constructed as follows. One considers a variable Φ A (ℓ) for each location ℓ ∈ L, and then the greatest solution of the system associated with the declaration D defined by:

Φ A (ℓ) D = a∈Σ    (ℓ,g,a,ℓ ′ )∈Out(ℓ,a) g, a Φ A (ℓ ′ ) ∧[En(ℓ, a), a] ℓ ′ ∈F(ℓ,a) Φ A (ℓ ′ ) ∧[¬En(ℓ, a), a]ff    (1)
These definitions should verify the following correctness property: for any ERA B, one has

B |= D Φ A (ℓ 0 ) if and only if A ∼ B.
Note that the construction introduces as clock constraints formulae obtained by disjunctions and negations. They can be rewritten in the syntax of ERL using the property [

g 1 ∨ g 2 , a]ϕ ≡ [g 1 , a]ϕ ∧ [g 2 , a]ϕ.
Before proving that the construction is not correct, we give some intuition on how it fails. To express bisimulation for a finite state automaton A, the standard approach consists in building a formula Φ A (q) for each state q of A, and considering the greatest solution of this system. Roughly, this formula verifies that any behaviour of A can be performed, and conversely that any possible behaviour corresponds to some of A. More formally, the standard formula for state q looks like:

Φ A (q) = a∈Σ       q a -→q ′ ∈A a Φ A (q ′ )    ∧   [a] q a -→q ′ ∈A Φ A (q ′ )       (2)
This is the way characteristic formulae for bisimulation are defined for instance in [START_REF] Laroussinie | From timed automata to logic -and back[END_REF][START_REF] Aceto | Characteristic formulae for timed automata[END_REF]. In the construction of [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF], the first conjunct corresponds to the first part of (2) while the two other conjuncts correspond to the second part of (2). But we can see that both parts are not well encoded. In the first one, notice that the constraint g, a Φ A (ℓ ′ ) implies the existence of at least one time successor in g that corresponds to the edge while all time successors in g should be able to fire this edge. In the second part, it is required that all a-successors occuring in En(ℓ, a) correspond to some a-successor of ℓ. But the a-successors of ℓ may have different clock constraints, and thus should not all be allowed in the whole set En(ℓ, a). We will see in Section 5 that the first point can be solved using the richer logic WT µ , and that the second point can be solved using the region construction.

ℓ ℓ ′ 0 ≤ xa ≤ 1 a A xa = 1 a B Fig. 1.
A counter-example to [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF].

We provide a counter-example exhibiting the first aspect. Consider the two ERA depicted on Figure 1. It is easy to see that A and B are not timed bisimilar. Let us write the formulae for A (Σ = {a}) according to (1):

Φ A (ℓ) = 0 ≤ x a ≤ 1, a [tt, a]ff ∧ [0 ≤ x a ≤ 1, a][tt, a]ff ∧ [x a > 1, a]ff We have B |= D Φ A (ℓ)
, what shows that the construction is not correct. More precisely, this is due to the uncompleteness of the first part of the formula of (1).

Consider the two ERA depicted on Figure 2. It is easy to verify that A and B are not timed bisimilar. However, the formulae for the ERA A according to (1) (with Σ = {a}) are:

Φ A (ℓ 0 ) = 0 ≤ x a ≤ 1, a Φ A (ℓ 1 ) ∧ 1 ≤ x a ≤ 2, a Φ A (ℓ 2 ) ∧[0 ≤ x a ≤ 2](Φ A (ℓ 1 ) ∨ Φ A (ℓ 2 )) ∧ [x a > 2, a]ff Φ A (ℓ 1 ) = x a = 0, a Φ A (ℓ 3 ) ∧ [x a = 0, a]Φ A (ℓ 3 ) ∧ [x a > 0, a]ff Φ A (ℓ 2 ) = [tt, a]ff Φ A (ℓ 3 ) = [tt, a]ff
One can verify that B, q B 0 |= D Φ A (ℓ 0 ) and thus the construction fails. It is worth noticing here that this is due the constraint 2. A second counter-example to [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF].

[0 ≤ x a ≤ 2](Φ A (ℓ 1 ) ∨ Φ A (ℓ 2 )) which is not enough restrictive. ℓ0 ℓ1 ℓ2 ℓ3 0 ≤ xa ≤ 1, a 1 ≤ xa ≤ 2, a xa = 0, a A 0 ≤ xa ≤ 2, a 1 ≤ xa ≤ 2, a xa = 0, a B Fig.

Impossibility Result for ERL

The construction of [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF] was proposed for the class of ERA without invariants.

It would be rather easy to prove that the logic ERL cannot express timed bisimilarity for ERA with invariants as this logic can not quantify over time elapsing independantly of the firing of a discrete transition. We prove here a stronger result by showing that the logic ERL can not express timed bisimulation for the restricted class of ERA without invariants. Following previous discussion, the logic ERL lacks a way to require the existence of a discrete transition for all the time successors satisfying some clock constraint. We will use this remark to prove the following main result:

Theorem 1. The logic ERL can not express timed bisimilarity for ERA, even without invariants.

Proof. We consider the ERA A depicted on Figure 1 and proceed by contradiction. We thus assume that there exists a ERL formula ϕ characterizing A up to timed bisimilarity. To simplify the presentation, we assume that A is defined over the alphabet restricted to letter a, but the result would hold for any alphabet. We denote by d ∈ N >0 the greatest denominator of constants appearing in clock constraints of ϕ. In the sequel, we call granularity of ϕ the value 1 d . In a first step, we simplify the formula ϕ. By Knaster-Tarski theorem, we have the following equalities: (for any ERL formula Φ, any ERA B and any V)

µX.Φ(X) B V = i≥0 Φ i (ff) B V ; νX.Φ(X) B V = i≥0 Φ i (tt) B V
As mentionned before, we can assume that all variables of sentences of ERL are guarded, i.e. are under the scope of the operator • or [•]. A consequence is that when interpreting fixpoints over structures without loops, one can limit above infinite disjunctions and conjunctions up to the maximal length of executions of the structure. For an ERA whose maximal depth 5 is 1 (such as A for instance), we can replace in ϕ the fixpoints operators by the above equations with index i ranging over the set {0, 1, 2}. We denote by Unfold 1 this operation, and by ERA d≤1 the set of ERA whose maximal depth is smaller or equal to 1. Then, we have:

∀B ∈ ERA d≤1 , B |= ϕ ⇐⇒ B |= Unfold 1 (ϕ) (3)
Thus, the outermost operators of the formula Unfold 1 (ϕ) belong to the set {∨, ∧, • , [•]}. We can then transform the formula Unfold 1 (ϕ) in a standard conjunctive normal form and write Unfold 1 (ϕ) = k i=1 mi j=1 Φ i,j where every formula Φ i,j has as outermost operator either • or [•]. Now, as the ERA A is of maximal depth 1 and is naturally timed bisimilar to itself, it satisfies this formula in its initial configuration q A 0 , and thus there exists i ∈ {1, . . . , k} such that A, q A 0 |= Φ i,j for any j ∈ {1, . . . , m i }. To ease the reading, we omit in the sequel the index i. Up to a reordering of the formulae Φ j , we can suppose that there exists an index p such that a formula Φ j has as outermost operator the operator • if and only if j ≤ p.

In this second part, we define an ERA B which is not not bisimilar to A. This ERA B is defined over Σ = {a} and contains exactly two locations, denoted respectively ℓ 1 and ℓ ′ 1 , such that the first one is initial. We denote by q B 0 = (ℓ 1 , 0) the initial configuration of B. In the sequel, we will define a finite set of rational numbers F. We exactly add one edge (ℓ 1 , g f , a, ℓ ′ 1 ) for each f ∈ F, with the constraint g f defined as x a = f . It is easy to verify that A and B are not timed bisimilar as there necessary exists some point in the interval [0, 1] that does no belong to F. We now detail how we build the set F to ensure that B, q B 0 |= ϕ. For each j ∈ {1, . . . , p}, we can write Φ j = g j , a ξ j for some constraint g j and formula ξ j . By construction, we have A, q A 0 |= Φ j , and thus there exists a delay δ ∈ T such that the steps q A 0 δ -→ (ℓ, δ) a -→ (ℓ ′ , 0) exist in A with A, (ℓ ′ , 0) |= ξ j . Note that independantly of the delay after which the a-labelled transition is fired, the configuration reached is the same. As the constraint g j is defined with

granularity 1 d ∈ Q >0 , we can choose δ j ∈ Q ≥0 ∩ [0, 1] such that q A 0 δj -→ (ℓ, δ j ) a -→ (ℓ ′ , 0) with A, (ℓ ′ , 0) |= ξ j . Finally, the finite set of rational values F is defined as F = {δ j | 1 ≤ j ≤ p}.
It remains to prove that the ERA B satisfies the formula ϕ. As the maximal depth of B is 1, and using property (3), it is sufficient to prove that for any j, we have B, q B 0 |= Φ j . First consider formulae Φ j for j > p. In this case the formula is of the form [g j , a]ξ j . Then the property holds because any a-labelled transition firable from q B 0 in B also exists in A, leading to identical configurations (ℓ ′ , 0) and (ℓ ′ 1 , 0), with no actions available in ℓ ′ and ℓ ′ 1 . Second, we consider a formula Φ j with j ≤ p. In this case, the choice of the delay δ j ∈ F ensures that the

transitions q B 0 δj -→ (ℓ 1 , δ j ) a -→ (ℓ ′
1 , 0) exist in B and as A, (ℓ ′ , 0) |= ξ j , we also have B, (ℓ ′ 1 , 0) |= ξ j . 5 The maximal depth of an ERA denotes the length of a longest untimed path.

Finally, we have proven that B |= ϕ holds while A and B are not timed bisimilar, thus yielding a contradiction.

⊓ ⊔

A µ-calculus for Event-Recording Automata

We introduce here a new µ-calculus for ERA, called WT µ [START_REF] Nguena-Timo | The logic WTµ[END_REF]. This stands for Weak Timed µ-calculus, as it can be seen as a timed µ-calculus (as T µ [START_REF] Henzinger | Symbolic model-checking for real-time systems[END_REF] or L ν [START_REF] Laroussinie | From timed automata to logic -and back[END_REF]) devoted to the weak class of timed systems represented by ERA. Its definition differs from ERL in that it separates delay successors and discrete successors, as it is done for instance in the logic L ν . We prove in this section that it is strictly more expressive than the logic ERL and that it preserves the good model-checking properties of ERL. We will show in the next section that it allows to express timed (bi)similarity for ERA.

The Logic WT µ

Definition 10 (Syntax). Let Σ be a finite alphabet and Var be a finite set of variables. A formula ϕ of WT µ is generated using the following grammar:

ϕ ::= tt | ff | X | ϕ ∧ ϕ | ϕ ∨ ϕ | a ϕ | g ϕ | [a]ϕ | [g]ϕ | µX.ϕ | νX.ϕ
where g ∈ C(Σ), a ∈ Σ and X ∈ Var.

As for the logic ERL, we use auxiliary assignment functions, and the notions of free variable, bounded variable, and sentence. Definition 12. Given two sentences ϕ and ϕ ′ in ERL ∪ WT µ , we say that they are equivalent if and only if, for any ERA A, we have ϕ A = ϕ ′ A . We say that a logic L 2 is more expressive than a logic L 1 if for any sentence in L 1 , there exists an equivalent sentence in L 2 .

Then we can state the following property: Proposition 2. Given a sentence ϕ ∈ ERL, we denote by φ the sentence of WT µ obtained by substituting any operator [g, a] (resp. g, a ) by the two operators [g][a] (resp. g a ). Then ϕ and φ are equivalent.

Proof. Proceeding by induction on the length of the formula ϕ, the result directly follows from the definitions.

⊓ ⊔

In the sequel, we consider the ERA A 1 reduced to a single location ℓ with no invariant and no transitions. We denote by Q its set of configurations. Proposition 3. For any ERL sentence ϕ, we have ϕ A1 = ∅ or ϕ A1 = Q.

Proof. According to [START_REF] Sorea | A decidable fixpoint logic for time-outs[END_REF], one can assume that the formulae of ERL are guarded, that is every bounded variable is in the scope of an operator • or [•]. Then, since automaton A 1 contains no transition, we have, for any variable assignment

V, that [g, a]ϕ A1 V = Q = tt A1 V and g, a ϕ A1 V = ∅ = ff A1 V .
As a consequence, any ERL formula in which every variable is bounded (a sentence) is equivalent to a formula obtained by the following grammar:

ϕ ::= tt | ff | ϕ 1 ∧ ϕ 2 | ϕ 1 ∨ ϕ 2 | µX.ϕ | νX.ϕ.
The result follows by induction on the length of the formula. ⊓ ⊔ Theorem 2. The logic WT µ is strictly more expressive than the logic ERL (even for ERA without invariants).

Proof. First, Proposition 2 exactly proves that the logic WT µ is more expressive than the logic ERL.

Second, we have to prove that the converse is false. Consider the ERA A 1 introduced above, and note that A 1 has no invariants. We consider the WT µ sentence ϕ = x a = 1 tt. Then a configuration (ℓ, v) of A 1 satisfies sentence ϕ if and only if v(x a ) ≤ 1. On the other side, Proposition 3 shows that any ERL sentence has a trivial satisfiability set on A. Thus, no ERL sentence is equivalent to ϕ, proving the result. ⊓ ⊔

Model-Checking

We consider the model checking problem of WT µ sentences on ERA models. This problem consists, given a WT µ sentence ϕ and an ERA A, in deciding whether the relation A |= ϕ holds.

Theorem 3. The Model-Checking problem of WT µ on ERA is EXPTIME-complete, even for the fragment of WT µ restricted to greatest fixpoints.

Proof. EXPTIME-membership: This easily follows from the EXPTIME-membership of the model-checking of the logic L + µ,ν over timed automata [START_REF] Aceto | Is your model-checker on time ?[END_REF]. However, to obtain precise complexity results, we sketch here the reasonning.

We first state the following Lemma:

Lemma 1. Let Σ be a finite alphabet. Let A ∈ ERA, ϕ ∈ WT µ be a formula without fixpoint quantifier and K denote the maximal integer constant of A and ϕ. Denote by X 1 , . . . , X n the free variables X 1 , . . . , X n of ϕ and let V be an assignment function over these variables such that for any i, V(X i ) is a union of regions in R K (A). Then, the semantics ϕ A V is also a union of regions of R K (A).

Proof. We proceed by induction on the length of ϕ and consider the type of ϕ:

ϕ = tt or ϕ = ff. The result follows as Q A and ∅ are union of regions.

-ϕ = ϕ 1 ∧ ϕ 2 or ϕ = ϕ 1 ∨ ϕ 2 .
The result follows from the induction property as union of regions are closed under boolean operations. ϕ = X i for some i ∈ {1, . . . , n}. Then ϕ A V = V(X i ) and the result follows from the hypothesis on V.

ϕ = g ϕ ′ or ϕ = [g]ϕ ′ with g ∈ C(Σ). By induction property the semantics ϕ ′ A V is a union of regions. Then, the result follows from the time-abstract bisimulation property of clock regions which implies that the time predecessors of a clock region is a union of clock regions.

ϕ = a ϕ ′ or ϕ = [a]ϕ ′ with a ∈ Σ. By induction property the semantics ϕ ′ A V is a union of regions. Then, the result follows from the time-abstract bisimulation property of regions which implies that the predecessors of a region by a discrete transition is a union of regions. This concludes the proof.

⊓ ⊔

As a corollary, we get:

Lemma 2.
Let Σ be a finite alphabet. Let A ∈ ERA, and ϕ be a sentence in WT µ . Denote by K the maximal integer constant of A and ϕ. Then the semantics of ϕ over A, ϕ A , is a union of regions of R K (A). In other terms, we have:

∀ℓ ∈ L A , ∀v, v ′ ∈ T Σ s.t. v ≃ K v ′ , A, (ℓ, v) |= ϕ ⇐⇒ A, (ℓ, v ′ ) |= ϕ
Proof. As the semantics of formulae of WT µ leads to monotone functions, Knaster-Tarski theorem implies that fixpoint formulae can be evaluated using eventually infinite intersections and unions given by:

µX.ϕ(X) A V = i≥0 ϕ i (ff) A V , νX.ϕ(X) A V = i≥0 ϕ i (tt) A V
As ∅ and Q are both union of regions, Lemma 1 entails that the iterative evaluation of fixpoints leads also to union of regions. As the number of regions if finite, these evaluations terminate, returning also a union of regions.

⊓ ⊔

The proof of Lemma 2 thus yields that the model checking problem is decidable. To obtain results on complexity issues, we reduce the model checking problem to an equivalent model checking problem for standard µ-calculus. Therefore, we define the semantics of WT µ over R K (A). The only operators for which the semantics is non standard are the following:

g ϕ RK (A) V = {(ℓ, r) ∈ R K (A) | ∃r ′ ∈ R K (Σ) s.t. (ℓ, r) τ -→ (ℓ, r ′ ), r ′ ⊆ g and (ℓ, r ′ ) ∈ ϕ RK (A) V } [g]ϕ RK (A) V = {(ℓ, r) ∈ R K (A) | ∀r ′ ∈ R K (Σ) s.t. (ℓ, r) τ -→ (ℓ, r ′ ), if r ′ ⊆ g then (ℓ, r ′ ) ∈ ϕ RK (A) V }
Then, we can prove the correction of this semantics as in [START_REF] Henzinger | Symbolic model-checking for real-time systems[END_REF][START_REF] Laroussinie | From timed automata to logic -and back[END_REF]:

∀v ∈ T Σ , A, (ℓ, v) |= ϕ ⇐⇒ R K (A), (ℓ, [v]) |= ϕ
However, the semantics of WT µ over R K (A) does not exactly match this of standard mu-calculus. This is due to inclusion testing between r ′ and g . To solve this problem, we can for instance introduce atomic propositions corresponding to the clocks constraints g ∈ C(Σ) of the fomula ϕ. A predicate g is satisfied in a region (ℓ, r) if and only if the inclusion r ⊆ g holds. Then, we can write the following equivalences:

g ϕ ≡ τ (g ∧ ϕ); [g]ϕ ≡ [τ ](g → ϕ) ≡ [τ ](¬g ∨ ϕ)
Note that the number of atomic propositions introduced for a formula ϕ ∈ WT µ is linear in the size of this formula. Another approach consists in enlarging the alphabet to include the clock constraints. This approach is described in [START_REF] Nguena-Timo | The logic WTµ[END_REF]. Finally, we obtain the reduction desired to a model checking problem of the standard mu-calculus over the region automaton. This problem, for a mu-calculus formula ϕ and a finite structure S, can be solved in time O((|S|×|ϕ|) n+1 ), where n is the number of alternations of greatest and least fixpoints quantifiers in ϕ [START_REF] Thomas | Languages, automata, and logic[END_REF]. As the size of R K (A) is in |A|×2 O(|Σ|. log K|Σ|) , and n is in O(|ϕ|), we obtain that the model checking problem of WT µ over ERA is in EXPTIME, with a precise time complexity.

EXPTIME-hardness:

We adapt the proof of [START_REF] Aceto | Is your model-checker on time ?[END_REF] to encode the acceptance problem of a word w 0 by a Linear Bounded Alternating Turing Machine (LBATM) M which is EXPTIME-complete [START_REF] Chandra | [END_REF]. One can assume w.l.o.g that the alphabet of M is {a, b}, and let n = |w 0 |. Configurations of M are triples (q, w, i) where i ≤ n denotes the position of the tape head. A transition (q, α, α ′ , δ, q ′ ) of M can be fired from (q, w, i) iff w[i] = α. Then, it writes α ′ instead, and moves left or right according to δ. As M is alternating, Q is partitioned into Q or and Q and . A configuration (q, w, i) with q ∈ Q or (resp. q ∈ Q and ) is winning iff q = q f or there exists an accepting successor configuration (resp. if all its successor configurations are accepting). As we want to build an ERA A while the construction of [START_REF] Aceto | Is your model-checker on time ?[END_REF] is done for timed automata, we make some modifications to control the resets of clocks. Locations of A are pairs (q, i) ∈ Q × N, where i denotes the position of the tape head. The value of cell i of the tape is encoded by the relative values of two clocks, say x ai and x bi . The alphabet of A thus contains Σ = {a i , b i | 1 ≤ i ≤ n}. We add a letter τ not in Σ. A transition (q, α, α ′ , δ, q ′ ) is represented in A by the transitions (q, i) gi,σi ---→ (q ′ , i ′ ), where:

1. g i = x ai < x bi ∧ x τ = 1 if α = a, and g i = x ai > x bi ∧ x τ = 1 otherwise, 2. σ i = x ai if α ′ = a, and σ i = x bi otherwise, 3. i ′ = i + 1 if δ = R and i < n, and i ′ = i -1 if δ = L and i > 1.
To force time elapsing between two transitions corresponding to moves of M, we use letter τ and add transitions (q, i)

xτ =1,τ
-----→ (q, i) for any location (q, i). The initialization of the clocks to represent the word w 0 can be done using a sequence of transitions u i interleaved by transitions labelled by τ . Finally, we use the following WT µ formula, with only greatest fixpoints:

ϕ = [tt][u 1 ][τ ] . . . [tt][u n ][τ ].νX.([accept]ff ∧ [tt][Σ][τ ] tt Σ τ X)
where accept denotes a special letter only firable from the final state of M. Then one can prove that M accepts w 0 iff A |= ϕ. Note that the size of A and ϕ are polynomial in the sizes of M and w 0 .

⊓ ⊔ Remark 2. As in [START_REF] Aceto | Is your model-checker on time ?[END_REF], the hardness proof could be done without diagonal constraints.

Characteristic Formulae Constructions

We describe in this section formulae constructions in the logic WT µ to express timed similarity and timed bisimilarity for ERA with invariants. In the sequel, we consider an ERA A = L A , ℓ A 0 , Σ, T A , I A over the alphabet Σ. Let ℓ ∈ L A and a ∈ Σ, we first introduce an operation, denoted Split(ℓ, a), related to the determinization of ERA. Split(ℓ, a) is a finite set of constraints {g 1 , . . . , g n } ⊆ C(Σ) such that:

(i) it partitions the constraint En(ℓ, a): i g i = En(ℓ, a) and ∀i = j, g i ∩ g j = ∅, (ii) its elements "match" the clock constraints of a-labelled transitions leaving ℓ: ∀i ∈ {1, . . . , n}, ∀(ℓ, g, a, ℓ ′ ) ∈ T A , g i ⊆ g or g i ∩ g = ∅.

We do not investigate here how such an operator can be defined as it is not the purpose of this work. It can for instance be defined using the region construction, and then be optimized using some merging operations on zones. It is worth noticing that in the worst case, the size of Split(ℓ, a) may be |Out(ℓ, a)|×2 O(|Σ| log K|Σ|) , with K the largest integer constant of A (due to the region construction). However, if the ERA A is deterministic, then its size is linear in the size of Out(ℓ, a). Indeed, the determinisim implies that the clock constraints of a-labelled transitions leaving ℓ are disjoint.

Characteristic Formulae for Timed Bisimilation

Definition 13. We define a declaration D ∼A associating a formula to each location ℓ of A, and consider the greatest solution of this system of fixpoint equations.

Φ ∼A (ℓ) D∼A =                                        a∈Σ (ℓ,g,a,ℓ ′ )∈TA [g] a Φ ∼A (ℓ ′ ) (C 1 ) ∧ [I A (ℓ)] Φ ∼A (ℓ) (C 2 ) ∧ a∈Σ g∈Split(ℓ,a) [g][a] (ℓ,g ′ ,a,ℓ ′ )∈TA| g ⊆ g ′ Φ ∼A (ℓ ′ ) (C 3 ) ∧ a∈Σ [¬En(ℓ, a)][a]ff (C 4 ) ∧ [¬I A (ℓ)] ff (C 5 )
Before proving the correctness of this construction, we give some intuition on its definition. Let B be an ERA and analyze how these formulae constrain B. The parts C 1 and C 2 express the simulation constraints (A ≺ B), while the three other constraints express the converse (B ≺ A). More precisely, note that C 1 requires that any discrete transition of A also exists in B: for any transition in A and for all delays after which it is firable, there exists a corresponding transition in B leading to a bisimilar configuration. This combination of a universal quantification over delays with an existential quantification over discrete successors was missing in ERL, as shown in Section 3. In the converse direction, discrete transitions are encoded in C 3 and C 4 . C 4 states that an a transition can only happen in B when it is possible in A. C 3 uses the decomposition Split(ℓ, a) of the clock constraint En(ℓ, a) to express that any a transition in B corresponds to some a transition of A firable from the same valuation. This corrects the corresponding constraint of the construction of [START_REF] Sorea | Verification of Real-Time Systems through Lazy Approximations[END_REF] (see Section 3). Finally, C 2 and C 5 handle the case of delay transitions.

Remark 3 (On the size of formulae Φ ∼A ). Due to the use of the operator Split, these formulae are in the worst case of size |A|×2 O(|Σ| log K|Σ|) , with K the largest integer constant of A, whereas if A is deterministic, then their size is linear in the size of A. We believe that this exponential blow-up is not avoidable. Indeed, for ERA, once a discrete transition labelled by a has been fired, one can not recover the value of clock x a before this firing as it has been reset. Formulae of [START_REF] Aceto | Characteristic formulae for timed automata[END_REF], which have a linear size, compare the clock valuation with the guards after the discrete firing. Moreover, note that this exponential blow-up has no consequences on the theoretical time complexity of timed bisimilarity checking (see Corollary 1), as linear formulae would lead to the same complexity.

The following result states the correctness of the previous construction. Theorem 4. Let A and B be two ERA over Σ and consider ℓ and m two locations of A and B respectively. Then for any valuation v ∈ T Σ , we have :

(ℓ, v) ∼ (m, v) ⇐⇒ B, (m, v) |= D∼A Φ ∼A (ℓ)
In particular, we have: A ∼ B ⇐⇒ B |= D∼A Φ ∼A (ℓ A 0 ) Proof. To prove Theorem 4 we establish successively the two implications:

⇐ If B, (m, v) |= D∼A Φ ∼A (ℓ), then we have (ℓ, v) ∼ (m, v). ⇒ If (ℓ, v) ∼ (m, v), then B, (m, v) |= D∼A Φ ∼A (ℓ) holds.
Let us denote by Q A and Q B the set of configurations of A and B respectively. Proof of ⇐. We consider the relation R

⊆ Q A ×Q B defined as R = {((ℓ, v), (m, v)) | B, (m, v) |= D∼A Φ ∼A (ℓ)}
and show that it is a timed bisimulation. In other terms, we must verify the conditions of Definition 8.

(i) Step in A. Consider σ ∈ Σ ∪ T such that (ℓ, v) σ -→ (ℓ ′ , v ′ ) in A, and show that there exists m ′ ∈ L B such that (m, v) σ -→ (m ′ , v ′ ) in B and (ℓ ′ , v ′ )R(m ′ , v ′ ).
We distinguish two cases according to the nature of σ.

-If σ = a ∈ Σ. Then there exists a transition (ℓ, g, a, ℓ ′ ) ∈ T A corresponding to this firing. In particular, we have v |= g and v ′ = v[x a := 0]. By hypothesis, we have B, (m, v) |= D∼A Φ ∼A (ℓ). In particular the transition of A corresponds to a conjunct in part C 1 of Φ ∼A (ℓ), and we thus have

B, (m, v) |= D∼A [g] a Φ ∼A (ℓ ′ ). As v |= g, this implies the ex- istence of a step (m, v) a -→ (m ′ , v ′′ ) in B, with B, (m, v ′′ ) |= D∼A Φ ∼A (ℓ ′ ). The semantics of ERA implies that v ′′ = v[x a := 0], and then v ′′ = v ′ , what concludes this case. -If σ = δ ∈ T. Then we have (ℓ, v) δ -→ (ℓ, v + δ) in A what implies that v + δ |= I A (ℓ). Part C 2 of Φ ∼A (ℓ)
′ ∈ L A such that (ℓ, v) σ -→ (ℓ ′ , v ′ ) in A and (ℓ ′ , v ′ )R(m ′ , v ′
). Again, we distinguish two cases according to the nature of σ.

-If σ = a ∈ Σ. By hypothesis, we have B, (m, v) |= D∼A Φ ∼A (ℓ). In particular, part C 4 of this formula is satisfied what implies that v |= En(ℓ, a). Then, as Split(ℓ, a) partitions the constraint En(ℓ, a), there exists a unique clock constraint g ∈ Split(ℓ, a) such that v |= g. The corresponding conjunct of part

C 3 implies that B, (m ′ , v ′ ) |= D∼A (ℓ,g ′ ,a,ℓ ′ )∈TA| g ⊆ g ′ Φ ∼A (ℓ ′ ).
The second property of Split(ℓ, a) implies, as g is not empty, that there exists a transition (ℓ, g ′ , a, ℓ ′ ) ∈ T A such that B, (m ′ , v ′ ) |= D∼A Φ ∼A (ℓ ′ ) and with g ⊆ g ′ . As a consequence, we have v |= g ′ and then (ℓ, v)

a -→ (ℓ ′ , v ′′ ) in A, with v ′′ = v[x a := 0] = v ′ ,
what concludes this case.

4. ξ = [¬En(ℓ, a)][a]ff. By contradiction, assume that the property is not satisfied, that is that there exists a delay δ ∈ T such that v + δ ∈ En(ℓ, a) and (m, v + δ) a -→ (m ′ , v ′ ) in B for some configuration (m ′ , v ′ ). By bisimulation, an a-labelled transition is also firable from the configuration (ℓ, v + δ), what contradicts the fact that v + δ ∈ En(ℓ, a). 5. ξ = [¬I A (ℓ)]ff. By contradiction, assume that the property is not satisfied, that is that there exists a delay δ ∈ T such that v + δ |= I A (ℓ) and (m, v)

δ -→ (m, v + δ) in B. By bisimulation, we also have (ℓ, v) δ -→ (ℓ, v + δ) in B what contradicts v + δ |= I A (ℓ).
This concludes the proof of the property (4), and thus the second implication also holds. Proof. Using the previous theorem, this problem reduces to the model checking problem of B against formula Φ ∼A (ℓ A 0 ) under the declaration D ∼A . Note that Φ ∼A contains only greatest fixpoints and thus is alternation-free. As there exists better complexity results for this class (see [START_REF] Cleaveland | A linear-time model-checking algorithm for the alternation-free modal mu-calculus[END_REF]), the proof of Theorem 3 shows that the time complexity of this problem is in O(|R K (B)| × |Φ ∼A |). The result follows from the size of R K (B) and previous remarks on the size of the formulae Φ ∼A . ⊓ ⊔ 5.2 Characteristic Formulae for Timed Simulation Definition 14. We define a declaration D ≻A associating a formula to each location ℓ of A, and consider the greatest solution of this system of fixpoint equations.

Φ ≻A (ℓ)

D≻A =       
a∈Σ (ℓ,g,a,ℓ ′ )∈T

[g] a Φ ≻A (ℓ ′ ) (C ′ 1 )

∧ [I A (ℓ)] Φ ≻A (ℓ) (C ′ 2 )
Note that this construction leads to formulae of size linear in the size of A. The following result states the correctness of the previous construction.

Theorem 5. Let A and B be two ERA over Σ and consider ℓ and m two locations of A and B respectively. Then for any valuation v ∈ T Σ , we have :

(ℓ, v) ≺ (m, v) ⇐⇒ B, (m, v) |= D≻A Φ ≻A (ℓ)
In particular, we have: A ≺ B ⇐⇒ B |= D≻A Φ ≻A (ℓ A 0 ) We omit the proof as it is similar to that of Theorem 4. As for bisimilarity, we obtain an EXPTIME procedure to decide timed similarity: Moreover, this procedure can also be used to decide language inclusion between ERA. More precisely, we have: Corollary 3. Given two ERA A and B, the procedure checking timed simulation leads to an EXPTIME procedure to decide whether L(A) ⊆ L(B) holds or not.

Proof. We first determinize automaton B, resulting in B ′ . Following [START_REF] Alur | A determinizable class of timed automata[END_REF], the number of locations and transitions of B ′ is then exponential in the size of B. Using Proposition 1, language inclusion reduces to A ≺ B ′ , and then to the model checking problem B ′ |= D≻A Φ ≻A (ℓ A 0 ). Using previous analysis, this can be checked in time |R K (B ′ )| × |Φ ≻A |. Finally, we obtain a procedure to decide this language inclusion in time |A| × 2 |B| , which belongs thus to EXPTIME. ⊓ ⊔ Note that the problem of language inclusion is PSPACE-complete [START_REF] Alur | A determinizable class of timed automata[END_REF], thus this procedure is not optimal. However, the known algorithm [START_REF] Alur | A determinizable class of timed automata[END_REF] matching the lower bound consists in guessing a path in the region automaton. A zone-based version of this procedure may thus be an interesting alternative in practice.

Conclusion

In this paper, we focused on the construction of characteristic formulae for ERA up to timed (bi)similarity. After having shown that the problem could not be solved in the logic ERL, we have introduced the new logic WT µ , and have proven that it is strictly more expressive than ERL and that its model checking problem over ERA is EXPTIME-complete. We have finally provided characteristic formulae constructions in WT µ for the whole class of ERA with invariants.

Compared to existing results of [START_REF] Aceto | Characteristic formulae for timed automata[END_REF] for timed automata which can also be applied to ERA using natural translations, we obtain procedures in the same class of complexity (EXPTIME), but our time complexity are more precise. For instance, for a fixed alphabet Σ and if constants are encoded in unary, then timed (bi)simulation can be checked in polynomial time! Moreover, our algorithm for model checking WT µ against ERA should also be more efficient than going through L ν and timed automata as it involves only one copy of the event-clocks. Finally, we obtain a non-optimal procedure for inclusion checking between ERA, which we believe could lead to good results in practice.

As future work, we plan to study how the good decidability results of the satisfiability problem for ERL transfer to WT µ . We also envisage to adapt the implementation of the procedures of [START_REF] Aceto | Characteristic formulae for timed automata[END_REF] done in the tool CMC [START_REF] Laroussinie | CMC: A tool for compositional model-checking of real-time systems[END_REF] to this framework for ERA.

  then implies the existence of the transition (m, v) δ -→ (m, v + δ) in B, such that B, (m, v + δ) |= D∼A Φ ∼A (ℓ), as desired. This shows that the relation R is a timed simulation between A and B. (ii) Step in B. Conversely, we show that the relation R -1 is a timed simulation between B and A. As above, let us consider σ ∈ Σ ∪ T such that (m, v) σ -→ (m ′ , v ′ ) in B, and show that there exists ℓ

Corollary 1 .

 1 This concludes the proof of Theorem 4.⊓ ⊔ One can decide timed bisimilarity of two ERA A and B over Σ in time |A| × |B| × 2 O(|Σ| log K|Σ|) (K denotes the largest constant of A and B).

Corollary 2 .

 2 One can decide timed similarity of two ERA A and B over Σ in time |A| × |B| × 2 O(|Σ| log K|Σ|) (K denotes the largest constant of A and B).

Note that we extend the definition of[START_REF] Sorea | A decidable fixpoint logic for time-outs[END_REF] to ERA with invariants.

 Definition 11 (Semantics). For a given ERA A = L, ℓ 0 , Σ, T, I with associated TTS S A = Q, q 0 , Σ, → , a given formula ϕ ∈ WT µ , and an assignment function V : Var → P(Q), the set of states satisfying the fomula, denoted ϕ A V , is inductively defined as follows:

The cases of atomic and boolean formulae are standard.

Expressivity

As expected, the logic WT µ increases the expressive power of ERL. Before proving this main result, we introduce some definitions.

-

This concludes the proof that R -1 is also a timed simulation between B and A, and thus R is a timed bisimulation as desired. This concludes the proof of the first implication.

Proof of ⇒. We assume that the property B, (m, v) |= D∼A Φ ∼A (ℓ) holds and want to show that the two configurations (ℓ, v) and (m, v) are timed bisimilar. Recall that the formulae Φ ∼A (ℓ) are defined as the greatest solution of a system of inequations. Using the notion of coinduction [START_REF] Sangiorgi | Bisimulation: From the origins to today[END_REF], any solution of these inequations also satisfies these formulae. We consider the assignment function V over the variabes Φ ∼A (ℓ) defined by

It is then sufficient to prove the following inclusions:

The proof proceeds by considering each conjunct ξ of D ∼A (Φ ∼A (ℓ)).

1. ξ = [g] a Φ ∼A (ℓ ′ ) for some transition (ℓ, g, a, ℓ ′ ) ∈ T A . We consider two cases whether this transition can be fired from the configuration (ℓ, v) or not. If it is not the case, that is ∀δ ∈ T, v +δ |= g, then we trivially have B, (m, v) |= ξ.

Otherwise, there exists a delay δ ∈ T such that v + δ |= g. Then, we have

. By bisimulation property and by time determinism, we have that (ℓ, v + δ) ∼ (m, v + δ) and then that there exists a configuration (m

Semantics of ERA implies that v ′ = v ′′ and thus the result follows since, by definition of V, we have (m

By bisimulation property and time determinism, we then have (ℓ, v + δ) ∼ (m, v + δ). This concludes this case. 3. ξ = [g][a] (ℓ,g ′ ,a,ℓ ′ )∈TA| g ⊆ g ′ Φ ∼A (ℓ ′ ), for some clock constraint g ∈ Split(ℓ, a).

Consider, if some exists, a delay δ ∈ T such that v + δ |= g and (m, v)

Then, we must show that the following holds: B, (m ′ , v ′ ) |= D∼A (ℓ,g ′ ,a,ℓ ′ )∈TA| g ⊆ g ′ Φ ∼A (ℓ ′ ). First, we have by bisimulation and time-determi-nism that (ℓ, v) δ -→ (ℓ, v + δ) exists in A and that (ℓ, v + δ) ∼ (m, v + δ) holds. Bisimulation then imples that there exists a transition (ℓ, v + δ) a -→ (ℓ ′ , v ′′ ) in B such that (ℓ ′ , v ′′ ) ∼ (m ′ , v ′ ). This implies that there exists a transition (ℓ, g ′ , a, ℓ ′ ) in T A such that v + δ |= g ′ . By the second property of Split(ℓ, a), this implies that g ⊆ g ′ , and thus this transition belongs to the disjunction of ξ. In particular, we thus have B, (m ′ , v ′ ) |= D∼A Φ ∼A (ℓ ′ ), as required.