Florence Chaubert-Pereira

Yann Guédon
email: guedon@cirad.fr

Christian Lavergne
email: lavergne@math.univ-montp2.fr

Catherine Trottier
email: trottier@math.univ-montp2.fr

Markov and semi-Markov switching linear mixed models 1 Markov and semi-Markov switching linear mixed models used to identify forest tree growth components

Keywords: Individual random effect, Markov switching model, MCEM algorithm, Plant structure analysis, Semi-Markov switching model

Tree growth is assumed to be mainly the result of three components: (i) an endogenous component assumed to be structured as a succession of roughly stationary phases separated by marked change points that are asynchronous between individuals, (ii) a time-varying environmental component assumed to take the form of synchronous fluctuations between individuals, (iii) an individual component corresponding mainly to the local environment of each tree. In order to identify and characterize these three components, we propose to use semi-Markov switching linear mixed models, i.e. models that combine linear mixed models in a semi-markovian manner. The underlying semi-Markov chain represents the succession of growth phases and their lengths (endogenous component) while the linear mixed models attached to each state of the underlying semi-Markov chain represent -in the corresponding growth phase-both the influence of time-varying climatic covariates (environmental component) as fixed effects, and inter-individual heterogeneity (individual component) as random effects. In this paper we address the estimation of Markov and semi-Markov switching linear mixed models in a general framework. We propose a MCEM-like algorithm whose iterations decompose into three steps: (i) sampling of state sequences given random effects, (ii) prediction of random effects given state sequences, (iii) maximization. The proposed statistical modeling approach is illustrated by the analysis of successive annual shoots along Corsican pine trunks influenced by climatic covariates.

Introduction

The analysis of plant structure at macroscopic scales is of major importance in forestry and different fields of agronomy; see [START_REF] Durand | Analysis of the plant architecture via tree-structured statistical models: The hidden Markov tree models[END_REF] and [START_REF] Guédon | Analyzing growth components in trees[END_REF] for illustrations. Tree development can be reconstructed at a given observation date from morphological markers corresponding to past events. Observed growth, as given for instance by the length of successive annual shoots along a tree trunk, is assumed to be mainly the result of three components: an endogenous component, an environmental component and an individual component. The endogenous component is assumed to be structured as a succession of roughly stationary phases that are asynchronous between individuals [START_REF] Guédon | Analyzing growth components in trees[END_REF] while the environmental component is assumed to take the form of fluctuations that are synchronous between individuals. This environmental component is thus assumed to be a "population" component as opposed to the individual component. The environmental factors that modulate plant development are mainly of climatic origin such as rainfall or temperature. The individual component may cover effects of diverse origins but always includes an individual genetic effect in the case of a tree population characterized by a genetic diversity. Other effects correspond to the local environment of each individual, e.g. competition between trees for light or nutrient resources.

A Gaussian hidden Markov chain (or Markov switching Gaussian model) [START_REF] Ephraim | Hidden Markov processes[END_REF][START_REF] Cappé | Inference in hidden Markov models[END_REF] can be defined as a pair of stochastic processes {S t , Y t } where the output process {Y t } is related to the state process {S t }, which is a finite state Markov chain, by the Gaussian distribution Y t |S t = s t ∼ N (µ st , Γ 2 st ). Markov switching linear models extend the class of Gaussian hidden Markov chains by incorporating the influence of covariates as fixed effects in the output process; see [START_REF] Fruhwirth-Schnatter | Finite Mixture and Markov Switching models[END_REF] for an overview of Markov switching models with different applications. In the literature, hidden

Markov models with random effects in the output process have been used in a limited way. [START_REF] Altman | Mixed hidden Markov models: An extension of the hidden Markov model to the longitudinal data setting[END_REF] introduced Markov switching generalized linear mixed models and applied these models to longitudinal data; see also [START_REF] Rijmen | Qualitative longitudinal analysis of symptoms in patients with primary and metastatic brain tumours[END_REF]. Both [START_REF] Altman | Mixed hidden Markov models: An extension of the hidden Markov model to the longitudinal data setting[END_REF] and [START_REF] Rijmen | Qualitative longitudinal analysis of symptoms in patients with primary and metastatic brain tumours[END_REF] assumed that the individual random effect and its variance are independent of the non-observable states.

A major drawback with hidden Markov models is the inflexibility in describing the time spent in a given state which is assumed to be geometrically distributed. It is unlikely that such a type of implicit state occupancy distribution is an appropriate model for tree growth phases. [START_REF] Guédon | Analyzing growth components in trees[END_REF] applied hidden semi-Markov chains with simple observation distributions to forest tree growth data. Hidden semi-Markov chains [START_REF] Guédon | Estimating hidden semi-Markov chains from discrete sequences[END_REF] generalize hidden Markov chains with the distinctive property of explicitly modeling the sojourn time in each state. Here we introduce semi-Markov switching linear mixed models that generalize both hidden semi-Markov chains and Markov switching linear mixed models. In our application context, the underlying semi-Markov chain represents the succession of growth phases and their lengths while the linear mixed models attached to each state of the semi-Markov chain represent -in the corresponding growth phase-both the effect of time-varying climatic covariates as fixed effects and inter-individual heterogeneity as random effects. In this paper, we will consider two types of random effect: (i) an individual random effect (common across states), and (ii) individual state-specific random effects. In both cases, we assume that the random effect variance depends on the non-observable state. The objective is both to characterize the tree population and analyse the behavior of each individual within the population.

Since both the states of the underlying Markov (respectively semi-Markov) chain and the random effects are non-observable, Markov (respectively semi-Markov) switching linear mixed models involve two latent structures and remain difficult to estimate. [START_REF] Altman | Mixed hidden Markov models: An extension of the hidden Markov model to the longitudinal data setting[END_REF] proposed a deterministic and a stochastic approximation method for estimating Markov switching generalized linear mixed models. The deterministic approximation approach combines numerical integration by Gaussian quadrature and quasi-Newton methods, and relies on the fact that the hidden Markov model likelihood can be written as a product of matrices.

Since the hidden semi-Markov model likelihood cannot be written as a product of matrices, this deterministic approximation method cannot be transposed to the semi-Markovian case.

Moreover, the deterministic approximation approach can only be applied in the case of a few random effects. The stochastic approximation method is a Monte Carlo EM (MCEM)

algorithm [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithm[END_REF] where the M-step involves quasi-Newton routines. Altman underlined some limitations of the two proposed methods such as sensitivity to starting values, slowness to converge and a strong computation burden. Since conditional independence assumptions within a Markov switching generalized linear mixed model can be represented by a directed acyclic graph, [START_REF] Rijmen | Qualitative longitudinal analysis of symptoms in patients with primary and metastatic brain tumours[END_REF] proposed to implement the E-step of the EM algorithm by a junction tree algorithm. The M-step involves numerical integration by Gaussian quadrature and Fisher scoring methods. Since conditional independence assumptions within a hidden semi-Markov model cannot be efficiently represented by a directed acyclic graph, this method cannot be transposed to the semi-Markovian case. Moreover, the approaches proposed by [START_REF] Altman | Mixed hidden Markov models: An extension of the hidden Markov model to the longitudinal data setting[END_REF] and [START_REF] Rijmen | Qualitative longitudinal analysis of symptoms in patients with primary and metastatic brain tumours[END_REF] cannot be transposed to our context where the random effects may be attached to the states. [START_REF] Kim | Segmental hidden Markov models with random effects for waveform modeling[END_REF] proposed an estimation method for a "left-right" semi-Markov switching linear mixed model with individual state-specific random effects. Thus, the states are ordered and each state can be visited at most once. Their proposed method, which is basically an application of the EM algorithm based on a forward-backward algorithm for the E-step, relies heavily on the two specific model assumptions (state visited at most once and individual state-specific random effects). Its complexity is cubic in the sequence length (because of the computation of the marginal observation distributions for each possible state segment location).

For stochastic models involving a latent structure such as Markov or semi-Markov switching models, Markov chain Monte Carlo (MCMC) algorithms can be difficult to implement and suffer from problems such as poor mixing and the difficulty of diagnosing convergence; see [START_REF] Fearnhead | Computational methods for complex stochastic systems: a review of some alternatives to MCMC[END_REF] This paper is organized as follows. Markov and semi-Markov switching linear mixed models are formally defined in Section 2. The maximum likelihood estimation of both Markov and semi-Markov switching linear mixed models using the proposed MCEM-like algorithm is presented in Section 3. The semi-Markov switching linear mixed model is illustrated in Section 4 by the analysis of successive annual shoots along Corsican pine trunks. Section 5 consists of concluding remarks.

Model definition

Let {S t } be a Markovian model with finite state space {1, . . . , J}. A J-state Markov chain is defined by the following parameters:

• initial probabilities π j = P (S 1 = j), j = 1, . . . , J with j π j = 1,

• transition probabilities p ij = P (S t = j|S t-1 = i), i, j = 1, . . . , J with j p ij = 1.

A J-state semi-Markov chain is defined by the following parameters:

• initial probabilities π j = P (S 1 = j), j = 1, . . . , J with j π j = 1,

• transition probabilities -nonabsorbing state i: for each j = i, p ij = P (S t = j|S t = i, S t-1 = i), with j =i p ij = 1 and p ii = 0, -absorbing state i: p ii = P (S t = i|S t-1 = i) = 1 and for each j = i, p ij = 0.

• An explicit occupancy (or sojourn time) distribution is attached to each nonabsorbing state:

d j (u) = P (S t+u+1 = j, S t+u-v = j, v = 0, . . . , u -2|S t+1 = j, S t = j), u = 1, 2, . . . , U j
where U j denotes the upper bound to the time spent in state j and

U j u=1 d j (u) = 1.
We define as possible parametric state occupancy distributions binomial distributions, Poisson distributions and negative binomial distributions with an additional shift parameter h (h 1) which defines the minimum sojourn time in a given state; see Web Appendix A for definitions of these parametric state occupancy distributions.

Since t = 1 is assumed to correspond to a state entering, the following relation is verified

P (S t+1 = j, S t-v = j, v = 0, . . . , t -1) = d j (t) π j .
Let Y a,t be the observation and S a,t be the non-observable state for individual a (a = 1, . . . , N ), at time t (t = 1, . . . , T a ). Let Y Ta a,1 = (Y a,1 , Y a,2 , . . . , Y a,Ta ) denote the T a -dimensional vector of observations for individual a, and Y the T -dimensional vector of all observations with T = N a=1 T a . The vectors of non-observable states, S Ta a,1 and S, are defined analogously.

A Markov (respectively semi-Markov) switching linear mixed model can be viewed as a pair of stochastic processes {S a,t , Y a,t } where the output process {Y a,t } is related to the state process {S a,t }, which is a finite state Markov (respectively semi-Markov) chain, by a linear mixed model. We introduce two nested families of Markov (respectively semi-Markov) switching linear mixed models which differ in the assumptions made concerning inter-individual heterogeneity in the output process:

• Individual random effect:

Given state S a,t = s a,t , Y a,t = X a,t β sa,t + τ sa,t ξ a + ǫ a,t , (1)
ξ a ∼ N (0, 1), ǫ a,t | S a,t = s a,t ∼ N (0, σ 2 sa,t ).
The individual status (compared to the average individual) within the population is common to all the states. The random effect ξ a is thus common to all the states but the random effect standard deviation τ sa,t depends on the state.

• Individual state-specific random effect:

Given state S a,t = s a,t , Y a,t = X a,t β sa,t + τ sa,t ξ a,sa,t + ǫ a,t , (2)
ξ a,sa,t ∼ N (0, 1), ǫ a,t | S a,t = s a,t ∼ N (0, σ 2 sa,t ).
The individual status is different in each state. The random effect ξ a,sa,t (with the attached standard deviation τ sa,t ) depends thus on the state.

In these definitions, X a,t is the Q-dimensional row vector of covariates for individual a at time t. Given the state S a,t = s a,t , β sa,t is the Q-dimensional fixed effect parameter vector and σ 2 sa,t is the residual variance. The individuals are assumed to be independent.

In our context, the random effects ξ a or ξ a,sa,t are standardized (i.e. follow the standard Gaussian distribution). This formalization is in particular appropriate for the individual random effect model where the random effect variance changes with state; see [START_REF] Foulley | Heterogeneous variances in Gaussian linear mixed models[END_REF]. In the individual state-specific random effect models, the random effects for an individual a are assumed to be mutually independent (cov(ξ a,j , ξ a,j ′ ) = 0; j = j ′ ).

Including random effects in the output process removes the assumption that the successive observations for an individual are conditionally independent given the non-observable states.

In the individual state-specific random effect models, observations in different states for an individual a are assumed to be conditionally independent given states (for

t = t ′ , cov(Y a,t , Y a,t ′ |S Ta a,1 = s Ta a,1 ) = 0 if s a,t = s a,t ′ and cov(Y a,t , Y a,t ′ |S Ta a,1 = s Ta a,1 ) = τ 2 sa,t if s a,t = s a,t ′ ).
In the individual random effect models, observations in different states for an individual a are not assumed to be conditionally independent given states (cov(Y a,t , Y a,t ′ |S Ta a,1 = s Ta a,1 ) = τ sa,t τ s a,t ′ for t = t ′ ). In state j, the introduction of random effects makes it possible to decompose the total variance Γ 2 j = τ 2 j + σ 2 j into two parts: variance due to inter-individual heterogeneity τ 2 j and residual variance σ 2 j .

Maximum likelihood estimation with a Monte Carlo EM-like algorithm

The maximum likelihood estimation is presented in the case of Markov switching linear mixed models with individual state-specific random effects. The extension to semi-Markov switching linear mixed models and the transposition to individual random effect models are straightforward. Markov switching linear mixed model parameters can be divided into two categories: parameters π = (π j ; j = 1, . . . , J) and P = ( p ij ; i, j = 1, . . . , J) of the underlying Markov chain and parameters β = (β j ; j = 1, . . . , J), τ = (τ j ; j = 1, . . . , J)

and σ 2 = (σ 2 j ; j = 1, . . . , J) of the J linear mixed models. In the following, we denote by θ = (π, P , β, τ, σ 2 ) the set of parameters to be estimated.

Let ξ J a,1 = (ξ a,j ; j = 1, . . . , J) be the J-dimensional random effect vector for individual a.

The likelihood function of the observed data is given by

L(y; θ) = N a=1 s Ta a,1 f (s Ta a,1 , ξ J a,1 , y Ta a,1 ; θ) dξ J a,1 = N a=1 s Ta a,1 f (s Ta a,1 ; θ)f (ξ J a,1 ; θ)f (y Ta a,1 | s Ta a,1 , ξ J a,1 ; θ) dξ J a,1 , (3)
where s Ta a,1 means "sum over every possible state sequences of length T a for individual a".

Since both the states of the underlying Markov chain and the random effects are nonobservable, the EM algorithm [START_REF] Mclachlan | The EM algorithm and extensions[END_REF] at first sight appears to be a natural candidate for estimating Markov switching linear mixed models. Let us consider the complete-data log-likelihood where both the observed data y, the random effects ξ and the states s of the underlying Markov chain are observed:

log f (s, ξ, y; θ) = N a=1 log f (s Ta a,1 , ξ J a,1 , y Ta a,1 ; θ) = N a=1 log f (s Ta a,1 ; θ) + log f (ξ J a,1 ; θ) + log f (y Ta a,1 | s Ta a,1 , ξ J a,1 ; θ) = N a=1 J j=1 I s a,1 = j log π j + N a=1 Ta t=2 J i,j=1 I s a,t = j, s a,t-1 = i log p ij + N a=1 J j=1
log φ(ξ a,j ; 0, 1)

+ N a=1 Ta t=1 J j=1 I s a,t = j log φ(y a,t ; X a,t β j + τ j ξ a,j , σ 2 j ), (4)
where φ(y; µ, σ 2 ) is the Gaussian density with mean µ and variance σ 2 , and I() is the indicator function.

The E-step of the EM algorithm requires calculating the conditional expectation of log f (s, ξ, y; θ) given the observed data y and the current value of θ. But the EM algorithm for hidden Markov models cannot be applied because the successive observations for an individual are not conditionally independent given the non-observable states; see Section 2.

The EM algorithm for a finite mixture of linear mixed models [START_REF] Celeux | Mixture of linear mixed models for clustering gene expression profiles from repeated microarray experiments[END_REF] cannot be adapted because the distributions of ξ a,j | Y Ta a,1 = y Ta a,1 and ξ a,j | S a,t = j, Y Ta a,1 = y Ta a,1 cannot be analytically derived for each individual a at time t. Thus, the Monte Carlo EM algorithm [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithm[END_REF], where the quantities computed in the deterministic E-step are approximated using Monte Carlo methods, provides an alternative approach.

For the presentation of the estimation algorithm, we adopted the framework of restorationmaximization (RM) algorithms proposed by [START_REF] Qian | Estimation of parameters in hidden Markov models[END_REF]. The MCEM algorithm proposed by [START_REF] Altman | Mixed hidden Markov models: An extension of the hidden Markov model to the longitudinal data setting[END_REF] can be interpreted as a RM algorithm with two restoration steps for the two latent structures, an unconditional stochastic one for the random effects and a conditional deterministic one for the state sequences (the unconditional/conditional qualifier refers to the other latent structure). We cannot adopt a similar approach with an unconditional restoration step since, in our definition of Markov switching linear mixed models (see Section 2), the random effects may be attached to the states. Hence, following [START_REF] Shi | Latent variable models with mixed continuous and polytomous data[END_REF], we chose rather to perform two conditional restoration steps, one for the state sequences given the random effects (and the observed data) and one for the random effects given the state sequences (and the observed data).

The proposed RM algorithm takes thus the following form:

1. Choose starting values θ (0) and ξ (0) for k = 0.

2. a. Conditional restoration step (R-step) for state sequences:

For each individual a, sample M k state sequences s Ta a,1 (m) from the conditional distribution P (S Ta a,1 = s Ta a,1 | ξ J a,1 (m), Y Ta a,1 = y Ta a,1 ; θ (k) ) using:

• a direct application of the forward-backward algorithm proposed by [START_REF] Chib | Calculating posterior distributions and modal estimates in Markov mixture models[END_REF] in the Markov switching model case; see Web Appendix B;

• a direct application of the forward-backward algorithm proposed by Guédon For each individual a, compute for each sampled state sequence, m = 1, . . . , M k , the best posterior prediction

ξ J a,1 (m) = E ξ J a,1 | S Ta a,1 = s Ta a,1 (m), Y Ta a,1 = y Ta a,1 ; θ (k) .

Prediction of random effects given state sequences

The predicted vector for the random effects ξ J a,1 (m) attached to the mth state sequence sampled for individual a is:

ξ J a,1 (m) = E ξ J a,1 | S Ta a,1 = s Ta a,1 (m), Y Ta a,1 = y Ta a,1 ; θ (k) = ΩU ′ a (m) U a (m)Ω 2 U ′ a (m) + diag(U a (m)σ 2 ) -1 y Ta a,1 - J j=1 I a,j (m)X a β j , (5)
where

• s Ta a,1 (m) is the mth state sequence sampled for individual a,

• Ω = diag(τ j ; j = 1, . . . , J) is the J × J random standard deviation matrix,

• U a (m) is the T a × J design matrix associated with state sequence s Ta a,1 (m), composed of 1 and 0, where u a,t (m) = I(s a,t (m) = 1), . . . , I(s a,t (m) = J) is the tth row with row sums, j U a (m)(t, j) = 1 and the sum of all elements, t j U a (m)(t, j) = T a ,

• σ 2 = (σ 2 1 , . . . , σ 2 J ) ′ is the J-dimensional residual variance vector,

• diag(U a (m)σ 2 ) is the T a ×T a diagonal matrix with (u a,t (m)σ 2 ; t = 1, . . . , T a ) on its diagonal,

• I a,j (m) = diag I(s a,t (m) = j), t = 1, . . . , T a is a T a × T a diagonal matrix,
• X a is the T a × Q matrix of covariates.

Maximization-step

The M-step is presented in the case of Markov switching linear mixed models. The Mstep for semi-Markov switching linear mixed models is presented in Web Appendix D. In the proposed MCEM-like algorithm, the conditional expectation of the complete-data loglikelihood (Equation 4) given the observed data is approximated at iteration k by

E log f (S, ξ, Y ; θ) | Y = y; θ (k) = N a=1 E log f (S Ta a,1 , ξ J a,1 , Y Ta a,1 ; θ) | Y Ta a,1 = y Ta a,1 ; θ (k) ≈ 1 M k N a=1 M k m=1 log f s Ta a,1 (m), ξ J a,1 (m), y Ta a,1 ; θ (k) ≈ 1 M k N a=1 M k m=1 J j=1 I s a,1 (m) = j log π (k) j + 1 M k N a=1 M k m=1 Ta t=2 J i,j=1 I s a,t (m) = j, s a,t-1 (m) = i log p (k) ij + 1 M k N a=1 M k m=1 J j=1
log φ(ξ a,j (m); 0, 1)

+ 1 M k N a=1 M k m=1 Ta t=1 J j=1 I s a,t (m) = j log φ(y a,t ; X a,t β (k) j + τ (k) j ξ a,j (m), σ 2(k) j ). ( 6
)
At iteration k, the new values for the parameters of the Markov switching linear mixed model are obtained by maximizing the different terms of Equation 6, each term depending on a given subset of θ.

For the parameters of the underlying Markov chain, we obtain:

• initial probabilities π (k+1) j = a m I s a,1 (m) = j N M k ,
• transition probabilities

p (k+1) ij = a m Ta t=2 I s a,t (m) = j, s a,t-1 (m) = i a m Ta t=2 I s a,t-1 (m) = i .
For the parameters of the J linear mixed models, we obtain:

• fixed effect parameters

β (k+1) j = a m X ′ a I a,j (m)X a -1 a m X ′ a I a,j (m) y Ta a,1 -τ (k) j ξ a,j (m) , (7)
• random effect standard deviations

τ (k+1) j = a m t I s a,t (m) = j ξ a,j (m) y a,t -X a,t β (k) j a m t I s a,t (m) = j ξ 2 a,j (m) , (8)
• residual variances

σ 2(k+1) j = a m y Ta a,1 -X a β (k) j -τ (k) j ξ a,j (m) ′ I a,j (m) y Ta a,1 -X a β (k) j -τ (k) j ξ a,j (m)
a m tr(I a,j (m))

.

These reestimation formulas are thus similar to standard maximum likelihood estimators for the different parameters.

For the computation of standard errors, [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithm[END_REF] proposed a Monte Carlo version of Louis method where the expectation is replaced by a Monte Carlo procedure when using a MCEM algorithm. The Wei and Tanner method can directly be applied in the proposed MCEM-like algorithm.

MCEM-like algorithm for individual random effect models

The application to individual random effect models (linear mixed model ( 1)) is straightforward. Since the individual random effects are incorporated in the output process, the main difference concerns the conditional R-step of random effect prediction given a state sequence (Equation 5). Using the notations introduced in Section 3.1, the predicted random effect ξ a for each individual a is given by

ξ a (m) = E ξ a | S Ta a,1 = s Ta a,1 (m), Y Ta a,1 = y Ta a,1 ; θ (k) = τ ′ U ′ a (m) U a (m)τ τ ′ U ′ a (m) + diag(U a (m)σ 2 ) -1 y Ta a,1 - J j=1 I a,j (m)X a β j ,
where τ = (τ 1 • • • τ J ) ′ is the J-dimensional random effect standard deviation vector. In the forward-backward algorithm (see Web Appendix B for Markov switching linear mixed models and Web Appendix C for semi-Markov switching linear mixed models) and in the M-step (Equations 7, 8 and 9), the random effects ξ J a,1 are replaced by ξ a .

Initialisation of the algorithm

Various simulations were conducted using different starting values. The more distant the starting values from true values, the worse the parameter estimates. We recommend choosing as starting values the parameters estimated by the EM algorithm for a simple Markov or semi-Markov switching linear model (i.e. without random effects).

Convergence of the algorithm

Under the assumption of convergence of random effect predictions, we chose to monitor the convergence of the proposed MCEM-like algorithm by the difference between successive

iterations log P (Y = y | ξ (k+1) ; θ (k+1) ) -log P (Y = y | ξ (k) ; θ (k) ). ( 10
)
The quantity log P (Y = y | ξ (k) ; θ (k) ) is obtained directly as a byproduct of the forward recursion.

Sample size

The conditional R-steps rely on the restoration of several pairs (s Ta a,1 , ξ J a,1 ) for each individual a. As discussed by [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithm[END_REF], it is inefficient to start with a large number of sampled state sequences M k . They recommended increasing M k as the current approximation moves closer to the true maximizer.

Application to Corsican pine growth

The use of semi-Markov switching linear mixed models is illustrated here by an analysis of forest tree growth. The data set comprised four sub-samples of Corsican pines planted in a forest stand in the "Centre" region (France): 31 6-year-old trees, 29 12-year-old trees (first year not measured), 30 18-year-old trees (first year not measured) and 13 23-year-old trees (three first years not measured). Tree trunks were described by annual shoot from the top to the base where the length (in cm) was recorded for each annual shoot (Figure 1).

The annual shoot is defined as the segment of stem established within a year. The trees were chosen in order to cover the entire range of behaviors and were not subject to any silvicultural interventions.

[Figure 1 about here.]

We applied the practical approach discussed in [START_REF] Guédon | Analyzing growth components in trees[END_REF] for determining the number of growth phases i.e. the number of states of the underlying non-ergodic semi-Markov chain. A "left-right" three-state semi-Markov switching linear mixed model composed of two successive transient states followed by a final absorbing state was estimated. Since the last year of measurement is arbitrary with regard to tree development, the length of the last growth phase is assumed to be systematically right-censored and cannot be modeled. In temperate regions, rainfall can have a one-year-delayed effect (on the number of leaves) or an immediate effect (on shoot elongation) depending on whether it occurs during organogenesis or elongation. We chose to use an intercept and the centered cumulative rainfall (in mm)

during a period covering one organogenesis period and one elongation period as fixed effects

for each linear mixed model. The linear mixed model attached to state j is:

Given state S a,t = j, y a,t = β j1 + β j2 X t + τ j ξ a,j + ǫ a,t , ξ a,j ∼ N (0, 1), ǫ a,t | S a,t = j ∼ N (0, σ 2 j ),
where y a,t is the length of the annual shoot for individual a at time t, β j1 is the intercept, X t is the centered cumulative rainfall at time t (E(X t ) = 0), β j2 is the cumulative rainfall parameter and ξ a,j is the random effect for individual a in state j. Because of the centering of the climatic covariate, the intercept β j1 is directly interpretable as the average length of successive annual shoots in state j.

[Figure 2 about here.]

The MCEM-like algorithm was initialized with parameters π, P , d, β and σ 2 estimated without taking the random effects into account (hence, ξ = 0). Once the random effects had converged, the convergence of the algorithm was monitored by the difference between two consecutive iterations of the observed data log-likelihood given the random effects (Equation 10). The plot of the values against the iteration rank showed that the estimation algorithm converged rapidly by, say, about 70 iterations with M k = k state sequences sampled for each tree at the kth iteration (Figure 2). It should be noted that various strategies were tested for increasing the number of sampled state sequences: linearly [START_REF] Mcculloch | Maximum likelihood variance components estimation for binary data[END_REF], in stages [START_REF] Mcculloch | Maximum likelihood algorithms for generalized linear mixed models[END_REF] and geometrically [START_REF] Caffo | Ascent-based Monte Carlo expectationmaximization[END_REF]. For the Corsican pine data set, the parameter estimates were not sensitive to the selected strategy.

The MCEM-like algorithm was initialized with a "left-right" model such that π j > 0 for each state j, p ij = 0 for j i and p ij > 0 for j > i for each transient state i. The fact that states 1 and 2 are the only possible initial states (with π 1 = 0.95 and π 2 = 0.05 at convergence) and that state 2 cannot be skipped (i.e. p 13 = 0 at convergence) is the result of the iterative estimation procedure; see Figure 3. The estimated transition probability matrix is thus degenerate i.e. for each transient state i, p i i+1 = 1 and p ij = 0 for j = i + 1. This deterministic succession of states supports the assumption of a succession of growth phases.

It should be noted that the estimated state occupancy distributions for states 1 and 2 are far from geometric distributions in terms of shape and relative dispersion; see Figure 3. This is an a posteriori justification of the semi-Markovian modeling of the growth phases.

The state occupancy distributions estimated for a Gaussian hidden semi-Markov chain (GHSMC) (i.e. without taking into account the effect of climatic covariates and interindividual heterogeneity) and for a semi-Markov switching linear mixed model (SMS-LMM) were compared. When the effect of climatic covariates and inter-individual heterogeneity were taken into account, this rendered the growth phases more synchronous between individuals; see the estimated state occupancy distributions, in particular their standard deviations, for the two models in Table 1 and Figure 3.

[Figure 3 about here.]

At convergence of the MCEM-like algorithm, the median predicted random effects were computed for each individual based on the random effects predicted for each state in each sampled state sequence. The most probable state sequence given the median predicted random effects was computed for each observed sequence using a Viterbi-like algorithm [START_REF] Guédon | Estimating hidden semi-Markov chains from discrete sequences[END_REF]. This restored state sequence can be viewed as the optimal segmentation of the corresponding observed sequence into sub-sequences, each corresponding to a given state.

The optimal segmentations of the observed sequences were used, in particular, to compute the mean centered cumulative rainfall and the average cumulative rainfall effect in each state (see below).

The marginal observation distribution of the linear mixed model attached to state j is the Gaussian distribution N (µ j , Γ 2 j ) with µ j = β j1 + β j2 E j (X) and Γ 2 j = τ 2 j + σ 2 j where E j (X)

is the mean centered cumulative rainfall in state j. The marginal observation distribution represents the length of the annual shoots in state j. The marginal observation distributions for the different states are well separated (little overlap between marginal observation distributions corresponding to two successive states); compare the mean difference µ j+1µ j between consecutive states with the standard deviations Γ j and Γ j+1 in Table 1. The fixed part of the three linear mixed models (i.e. β j1 + β j2 X t for each state j) for 18-year-old and 23-year-old trees is shown in Figure 4. This confirms that the states are well separated with little overlap and correspond to a growth increase.

[Table 1 about here.]

The average cumulative rainfall effect (i.e. the average amplitude of the climatic fluctuations) was computed as β j2 × mad j (X) for each state j where mad j (X) is the mean absolute deviation of the centered cumulative rainfall in state j. The effect of cumulative rainfall was weak in the first state (of slowest growth) while it was stronger in the last two states (a little less in the second state than in the third state); see Table 1. The proportion of interindividual heterogeneity, defined by the ratio between the random variance τ 2 j and the total variance Γ 2 j in state j, was greater in early plant life (first two states with approximately 55%) and decreased slightly in the last state (approximately 47%).

[Figure 4 about here.]

The correlation coefficient between the predicted random effect in state 1 and the predicted random effect in state 2 was 0.28 while the correlation coefficient between the predicted random effect in state 2 and the predicted random effect in state 3 was 0.61. Hence, the behavior of an individual is quite strongly related between the last two states but only loosely related between the first two states. A 95% prediction interval [START_REF] Hulting | Some Bayesian and non-Bayesian procedures for the analysis of comparative experiments and for small area estimation: Computational aspects, frequentist properties, and relationships[END_REF] was computed to check whether the influence of the predicted random effect for each state is significant:

[-t 0.975 (N -1) σ ξ j √ N ; t 0.975 (N -1) σ ξ j √ N ]
where σ ξ j is the empirical standard deviation of the predicted random effects for state j.

Of the 103 Corsican pines, 50 had a significant random effect in each state and of these 50 individuals, 29 showed a common behavior over all growth phases (i.e. growing either more rapidly or slowly than the "average tree" in all growth phases) but 21 showed different behaviors between growth phases. Hence, the more general assumption of individual statespecific random effect models (linear mixed model ( 2)) compared to individual random effect models (linear mixed model ( 1)) is more representative of Corsican pine behavior.

Complementary biological results concerning Corsican pine and sessile oak growth can be found in [START_REF] Chaubert-Pereira | Identifying ontogenetic, environmental and individual components of forest tree growth[END_REF].

Concluding remarks

In the proposed MCEM-like algorithm, the conditional restoration step for state sequences given random effects relies on simulations while the conditional restoration step for random effects given state sequences is deterministic. In this latter case, an alternative solution would be to sample random effects applying a Metropolis-Hastings algorithm; see [START_REF] Mcculloch | Maximum likelihood algorithms for generalized linear mixed models[END_REF].

The estimation algorithms proposed in this paper can be directly transposed to other families of hidden Markov models such as, for instance, hidden Markov tree models; see [START_REF] Durand | Analysis of the plant architecture via tree-structured statistical models: The hidden Markov tree models[END_REF] and references therein. Another interesting direction for further research would be to develop the statistical methodology for semi-Markov switching generalized linear mixed models to take into account categorical or count response variables (for instance, apex death/life, non-flowering/flowering character, number of growth units per annual shoot in the plant architecture context). Since the conditional expectation of random effects given state sequences cannot be derived analytically, the proposed MCEM-like algorithm for semi-Markov switching linear mixed models cannot be transposed to semi-Markov switching generalized linear mixed models, and other conditional restoration steps, for instance based on a Metropolis-Hastings algorithm, have to be derived for the random effects. Cumulative rainfall parameter β j2 cm.mm -1 2.7 10 -3 16.5 10 -3 30.9 10 -3 (s.e.)

(1.2 10 -3 ) (2 10 -3 ) (2.4 10 6.99, 3.24 25.88, 9.48 50.32, 12.09 Table 1 Comparison of the estimated Gaussian hidden semi-Markov chain (GHSMC) parameters (i.e. where the influence of covariates and the inter-individual heterogeneity are not taken into account) with the estimated semi-Markov switching linear mixed model (SMS-LMM) parameters (state occupancy distributions and marginal observation distributions). The regression parameters, the cumulative rainfall effect and the variability decomposition are given for each observation linear mixed model. Standard errors are given in brackets.

(

  2007) in the semi-Markov switching model case; see Web Appendix C. b. Conditional restoration step for random effects; see Section 3.1 :

  Ta a,1 (m), ξ J a,1 (m), y Ta a,1 ; θ(k) .d. Sample size increase step:Set k = k + 1 and increase the number of sampled state sequences. If M k+1 > M k , sample with replacement M k+1 random effects among the M k predicted random effect.3. If convergence is achieved, then declare θ (k+1) to be maximum likelihood estimator;otherwise, return to step 2.Since Markov (respectively semi-Markov) chain parameters and linear mixed model parameters form disjoint sets and influence the complete-data log-likelihood separately (see Equation 4 in the Markov switching model case and Web Equation 15 in Web Appendix D in the semi-Markov switching model case), Markov (respectively semi-Markov) chain parameters can be updated when the state sequences S are sampled, and linear mixed model parameters can be updated when the random effects ξ are predicted. It makes sense to re-estimate the parameters immediately before performing the conditional R-step for the other latent structure in order to speed up the convergence. The forward-backward algorithm for sampling state sequences given the random effects can be decomposed into two passes, a forward recursion which is similar to the forward recursion of the usual forward-backward algorithm and a backward pass for sampling state sequences; see Web Appendix B for Markov switching linear mixed models and Web Appendix C for semi-Markov switching linear mixed models.

Figure 1 .

 1 Figure 1. Length of successive annual shoots along Corsican pine trunks: (a) 6-year-old trees, (b) 12-year-old trees, (c) 18-year-old trees, (d) 23-year-old trees.

Figure 2 .Figure 3 .

 23 Figure 2. Difference of observed data log-likelihood given random effects between successive iterations (Equation10) for the estimated semi-Markov switching linear mixed model with individual state-specific random effects.

Figure 4 .

 4 Figure 4. Fixed part of the three observation linear mixed models (i.e. β j1 + β j2 X t for each state j; •: state 1, * : state 2, △: state 3) represented by point lines and observed annual shoot lengths (points): (a) 18-year-old trees, (b) 23-year-old trees.

Acknowledgements

The authors are grateful to Yves Caraglio and Céline Meredieu for providing the Corsican pine data and to Météo France for freely providing the meteorological data. The authors thank the Associate Editor and two referees for helpful comments that led to improve the presentation of this paper.


Odile Valiron
email: valiron<odile.valiron@ujf-grenoble.fr

Isabelle Arnal

Didier Nicolas Caudron

Didier Job

GDP-Tubulin Incorporation into Growing Microtubules Modulates Polymer Stability

Keywords: Animals, Biochemistry, methods, Centrosome, metabolism, Cryoelectron Microscopy, methods, Dimerization, Filtration, Guanosine Diphosphate, chemistry, Guanosine Triphosphate, chemistry, Humans, Hydrolysis, Microtubules, chemistry, metabolism, Nucleotides, chemistry, Polymers, chemistry, Tubulin, chemistry, Tubulin Modulators, chemistry

HAL is

Author Keywords Cytoskeleton ; Microtubules ; Nucleotide ; Protein ; Self-assembly ; Tubulin ; GTP Hydrolysis ; GTP-Tubulin ; Microtubule Dynamics In mammalian cells, microtubules are centrally involved in many vital processes such as cell morphogenesis and motility. Microtubule arrays display substantial variations in their dynamic behaviour, depending on the cell cycle or on the cell type, and this dynamic character is crucial to microtubular functions. The building blocks of microtubules are -tubulin heterodimers. Tubulin subunits associate laterally αβ and longitudinally into growing microtubules ( , ) in the form of either tubulin dimers or oligomers (

). Microtubules shorten 1 2

3 -5

through tubulin oligomer loss ( ). It is currently assumed that the GDP-tubulin subunits (GDP-tub) which build up the microtubule wall 6 originate from the incorporation of GTP-tubulin complexes (GTP-tub) followed by GTP hydrolysis in the polymer wall ( ). 7

Microtubule growth displays spontaneous transitions between growing and shrinking states, known as dynamic instability ( ). 8

Microtubule length variations are currently viewed as principally governed by the behaviour of the microtubule ends ( ). It has long been 9

assumed that GTP hydrolysis at the extremity of microtubules determined tubulin addition and loss ( ). Recent studies show that, for a 10 given nucleotide bound state of the tubulin dimers, microtubule length fluctuations are also dependent on structural events occurring at polymer ends. For example, microtubule tip-binding proteins such as EB1 can regulate dynamics and tip structure of microtubules assembled from purified tubulin (

) and structural differences between microtubule tips have been shown in studies on 11 kinetochore-microtubule interactions ( ). This reveals that microtubule ends can also experiment structural plasticity ( ). Although 12 " " 13 the switch-like behaviour of ends is essential for dynamic instability, microtubules could also exhibit structural plasticity along their length ( , ). Such plasticity could imply that the tubulin subunits constituting the microtubule wall exist in several structural and/or 13 14 biochemical states which may influence microtubule dynamic properties ( , ). 13 14 Although the dominant view is currently that microtubule growth proceeds exclusively from the incorporation of GTP-tub, there are scattered reports of direct GDP-tub incorporation into growing polymers ( , ). However, other studies indicate that GDP-tub does not 15 16 significantly participate to elongation ( ). Additionally, recent structural studies have revealed differences between GTP-tub and 17 -19 GDP-tub intra-dimer and inter-dimer interactions. This led to the suggestion that GDP-tub could not be directly incorporated in microtubules under any conditions ( ). Here, we have re-examined the possibility of a direct GDP-tub incorporation into growing 7 microtubules, using a minimal tubulin assembly system composed of nucleotide bound tubulin dimers, in the absence of excess free nucleotide ( , ). Within the framework of our study, such a minimal system had substantial advantage over usual systems in which 15 20 -22 tubulin assembles in the presence of excess of free GTP. The proportion of GTP-tub and GDP-tub added in solution could be controlled at will without the complication of excess free GTP competing with GDP for the tubulin nucleotide binding site; when GTP-tub assembles above the critical concentration in bulk assembly tests, microtubules undergo a phase of assembly followed by a phase spontaneous disassembly ( , ), which allow monitoring of both assembly and disassembly dynamics; additionally, in the present study, a 15 20 -22 similar minimal assembly system proved to be usable for study of individual microtubule dynamics at tubulin concentrations below the critical concentration, using short microtubule seeds or centrosomes to nucleate tubulin assembly. Using such minimal assembly systems, we show that substantial amounts of GDP-tub can be incorporated in growing microtubules during both microtubule bulk assembly and seed or centrosome nucleated microtubule assembly. Microtubules assembled from GTP-tub and GDP-tub mixtures (GTP-tub/GDP-tub mix) display altered dynamics. Our results suggest that the GDP-tub constituting the microtubule wall may be in different structural states according to their initial nucleotide bound state, with resulting variations in intrinsic microtubule disassembly properties. This may reveal a novel form of microtubule structural plasticity.

Experimental procedures

Tubulin preparation

Tubulin was purified from fresh bovine brain as described previously ( ). To prepare GTP-tub or H -GTP-tub, pure tubulin was 23

[ 3 ]
incubated in PEM buffer (100 mM Pipes (pH 6.7), 1 mM EGTA, 1 mM MgCl ), for 10 minutes at 4 C, in the presence of either 1 mM 2

°GTP

(for GTP-tub) or 0.5 mM GTP supplemented with 100 Ci/ M H -GTP (for H -GTP-tub). Free nucleotides were removed using

μ μ [ 3 ] [ 3 ]
Biogel P30 chromatography.

GDP-tub or H -GDP-tub was obtained as the cold disassembly product of microtubules initially polymerized from pure tubulin (100

[ 3 ]
M) during 20 min at 35 C in PEM buffer in the presence of 5 mM MgCl , 30 (v/v) glycerol, and either 1 mM GTP or 0. after dilution in PEM buffer at 4 C of pellets of microtubules assembled from respectively either GTP alone or GTP and H -GTP

°[3 ]
mixtures.

Tyr-tub was prepared according to Paturle ( ). Briefly, tubulin was equilibrated in MEM buffer (100 mM Mes, 1 mM EGTA, 1 mM 24 MgCl ) through Biogel P30 chromatography. Then tubulin (44 M final concentration) was incubated in MEM buffer containing ATP (5

Microtubule assembly conditions

GTP-tub, GTP-tub/GDP-tub mix or radioactive tubulin mixtures (either H -GTP-tub alone or H -GTP-tub/GDP-tub mix or

[ 3 ] [ 3 ]
GTP-tub/ H -GDP-tub mix) were aliquoted (20 l) at 4 C in tubes. Tubulin assembly was initiated by immersing tubes in a water bath at

[ 3 ] μ °35
C. At selected time points, reactions were stopped. Stop procedures and further processing of the samples were adapted to the °parameter to be measured, as described below (

) and in ( Filter assay Supplemental Results Analysis of microtubule nucleation and mean and ). All measurements were done in triplicate. length Microtubule sedimentation assay

Filter assay

To estimate polymeric radioactive tubulin concentration at selected time points, we used a previously developed filter assay ( ), with 25 minor modifications. Briefly, either H -GTP-tub alone or H -GTP-tub/GDP-tub mix or GTP-tub/ H -GDP-tub mix were prepared at 4

[ 3 ] [ 3 ] [ 3 ] °C,
aliquoted and assembled as above. Tubulin assembly was stopped by adding to samples 1 ml of 100 mM Mes (pH 6.7) containing 1 mM EGTA, 1 mM MgCl2, 0.75 glutaraldehyde and 50 (w/v) sucrose ( ). GF/F glass fiber filters (one filter per time point aliquot) % %

25

were placed in a vacuum filtration device, and washed with 4 ml of PEM buffer containing 25 glycerol (buffer PEM-G). Microtubule % suspensions were applied to filters under negative pressure. At this step, H -labelled cross-linked microtubules were trapped on filters

[ 3 ]
and the bulk of unassembled subunits passed through filters. Then filters were washed 3 times with 4 ml of PEM-G buffer. To extract H

[ 3 ]
-radioactivity from filters, they were incubated under shaking during 30 min with 2 ml ethanol in vials (1 vial per filter). Then liquid scintillation cocktail was added (10 ml per vial) and radioactivity was counted. Polymerized H -tubulin concentrations were estimated

[ 3 ]
from the radioactivity measured on filters, and from the specific activity of H -tubulin obtained during tubulin preparation. The H

[ 3 ] [ 3 ]
-tubulin specific activity (cpm/ mole of tubulin) was calculated as the ratio of the number of cpm contained in an aliquot of the starting μ tubulin solution over the amount of tubulin contained in the same aliquot.

Preparation of EGS cross-linked microtubule seeds (EGS-seeds)

Covalently cross-linked microtubule seeds were prepared using ethylene glycol-bis-succinimidylsuccinate (EGS) ( ). Tubulin (100 26 μ M) in a total volume of 200 l was assembled for 15 min at 35 C in 80 mM Pipes (pH6.7), 1 mM EGTA, 50 (v/v) glycerol, 5 mM MgCl μ °% 2 and 1 mM GTP. Microtubules were crosslinked by incubation for 15 min at 35 C after addition of 3.4 mM final concentration of EGS. To °quench the EGS in excess, the mixture was diluted into 1.8 ml of a buffer containing 80 mM Pipes (pH 6.7), 1 mM EGTA, 50 sucrose, % 10 mM glutamate and 1 mM MgCl , and incubated for 1 hr at room temperature. The solution was then centrifuged at 200,000g for 1 hr at 2 room temperature. The pellet containing cross-linked microtubules (EGS-seeds) was resuspended in 80 l of PEM buffer.

μ

Microtubule immunostaining

In microtubule self-assembly conditions (over-critical GTP-tub concentration), after stopping the assembly reaction microtubules were diluted in PEM buffer supplemented with 10 glycerol. They were then centrifuged on coverslips at 77,000g during 30 min at room % temperature and fixed with methanol for 6 min at 20 C. Coverslips were processed for indirect immunofluorescence analysis as -°previously described ( ) using primary anti Detyr-tub antibody ( ) and primary anti Tyr-tub antibody (clone YL1/2; 29). At sub-critical 27 28

GTP-tub concentration, microtubules were nucleated on centrosomes as described ( ) and immunolabelled as above. 30

Video-Microscopy and Data Analysis

Video-microscopy and analysis were performed as previously described ( ). Briefly, samples were prepared in perfusion chambers. 31

Purified centrosomes were first perfused into the chamber on ice. Samples (60 l) containing either GTP-tub or GTP-tub/GDP-tub mix μ were then perfused in the chamber and microtubule assembly was observed at 37 C under an Olympus BX51 microscope equipped with °Differential Interference Contrast prisms and a camera (Sony, XC-ST70/CE). Images were recorded every 2 sec and microtubule dynamics measurements and data analysis were performed using NIH-Image and Kaleidagraph softwares. For growth and shrinkage rates, standard deviations were calculated as standard error of the mean, assuming a normal data distribution. For catastrophe frequencies, standard deviations were calculated as catastrophe frequency/( n), where n is the number of events counted, assuming a Poisson distribution ( ).

√

32

RESULTS

GDP-tub incorporates in self-assembled microtubules and modulates their assembly-disassembly properties

We tested the incorporation of GTP-tub or GDP-tub into growing microtubules in a chemically simple system ( , ) in which 15 20 -22 microtubules were assembled from solutions containing GTP-tub/GDP-tub mix, in the absence of free nucleotide. Starting mixtures contained either H -GTP-tub and unlabelled GDP-tub or unlabelled GTP-tub and H -GDP-tub. We verified (

[ 3 ] [ 3 ]
Supplemental Results, Fig.

) that our tubulin preparations were devoid of detectable amounts of nucleotide-diphosphokinase (NDP) activity, which could induce S1 conversion of free GDP-tub to free GTP-tub during our experiments ( ). The incorporation of either GTP-tub or of GDP-tub in 33 microtubules assembled from GTP-tub/GDP-tub mix could then be selectively monitored by counting the H -nucleotide radioactivity

[ 3 ]
associated to microtubules trapped on filters (see Experimental Procedures, 22).

In the absence of GDP-tub, GTP-tub assembly follows a bell shaped curve, with a phase of microtubule assembly followed by a phase of microtubule disassembly ( ), as described before ( ). In our experiments, in agreement with previous reports ( , ), GDP-tub Fig. 1A 22 7 19 alone was unable to polymerize (not shown). Addition of GDP-tub to GTP-tub did not detectably modify the maximum of GTP-tub incorporation into microtubules ( ( )). Interestingly, when H -GDP-tub was mixed with GTP-tub in starting suspension, H Fig. 1A

D - ■ [ 3 ] [ 3 ]
-GDP-tub was incorporated in assembling microtubules ( ( Fig. 1B D -)). GDP-tub incorporation occurred in amounts proportional to the initial GDP-tub/GTP-tub ratio (

). As a result of the Fig. 1E incorporation of GDP-tub, the total tubulin assembly level increased at increasing initial GDP-tub concentrations ( ( Fig. 1A D -)). The assembly phase was prolonged whereas, in a quantitative analysis, initial microtubule nucleation and elongation seemed little affected (

). The disassembly phase was conspicuously prolonged at increasing GDP-tub concentrations, to Supplemental Results, Fig. S2 such an extent that instead of exhibiting a characteristic symmetry bell shaped aspect (

), the tubulin assembly-disassembly curve Fig. 1A became right skewed (

). Compared to the control, the bulk microtubule half disassembly time increased 5 fold at a 1.5/1 initial Fig. 1D GDP-tub/GTP-tub ratio (

). Fig. 1F In a series of control experiments, tubulin assembly was monitored using either turbidity measurements (Supplemental Results, Fig.

) or microtubules sedimentation assays (

). In turbidity assays, GDP-tub alone was unable to polymerize S3 Supplemental Results, Fig. S4 whereas addition of excess free GTP produced the expected sustained microtubule assembly ( ). The Supplemental Results, Fig. S1A assembly of GTP-tub/GDP-tub mix yielded assembly plots similar to those observed in the same conditions using filter assays (Fig. 1A D -), which was also the case of assembly plots derived from microtubule sedimentation assays ( ). Supplemental Results, Fig. S4 These results provide compelling evidence that GDP-tub can be directly co-incorporated into growing microtubules, together with GTP-tub. Additionally, our results suggest that GDP-tub incorporation in microtubules could impair microtubule disassembly.

GDP-tub co-assembly occurs along the whole polymer length but does not induce detectable changes in the microtubule lattice organization

To directly visualize GDP-tub incorporation into microtubules we used GDP-tub or GTP-tub made of different tubulin tyrosination variants. Previous work has shown that tubulin tyrosination variants such as tyrosinated tubulin (Tyr-tub) and detyrosinated tubulin (Detyr-tub) have distinct immunoreactivity properties ( ) but have indistinguishable in vitro assembly properties ( ). To test the 34 24 reliability of tubulin variants as reporters of the tubulin bound nucleotide state, GDP-Detyr-tub was mixed with GTP-Tyr-tub or, in symmetric experiments GTP-Detyr-tub was mixed with GDP-Tyr-tub. Tubulin mixtures were assembled at 35 C and sedimented °microtubules were assayed for their tubulin composition. The proportion of Detyr-tub or Tyr-tub incorporated in microtubules was a function of the initial Detyr-tub or Tyr-tub bound nucleotide state and was in good agreement with the ratio of GTP-tub/GDP-tub incorporation estimated from radioactive measurements ( And ). Thus, over the time course of our Fig. 2A Supplemental Results

experiments, the association of each tubulin variant with its bound nucleotide was stable enough to allow qualitative visualization of GDP-tub incorporation in growing microtubules. Microtubules assembled from Tyr-GTP-tub/Detyr-GDP-tub mix were then sedimented and double stained with Tyr-tub and Detyr-tub antibodies. Immunofluorescence analysis showed that polymers were uniformly and homogeneously stained by both antibodies ( ). These results indicate that GDP-tub is incorporated in growing microtubules during Fig. 2B the whole assembly process.

We then used cryo-electron microscopy to examine the structure of microtubules assembled either from GTP-tub alone or from GTP-tub/GDP-tub mix (

). Based on the moire pattern observed on microtubule images, we determined the microtubule Fig. 2C protofilament number and the frequency of lattice defects as previously described ( ,

). We also assessed the aspect of microtubule 35 36 ends in polymerizing and depolymerizing conditions and we analyzed the disassembly products of polymers exposed to high calcium concentrations ( , ). We found no significant difference between GTP-tub or GTP-tub/GDP-tub polymers (not shown). Thus, 37 38 incorporation of GDP-tub into microtubules did not induce any obvious changes in the polymer structure.

GDP-tub incorporated into microtubules nucleated from seed or centrosomes

The experiments shown above were performed at high initial tubulin concentrations, compatible with efficient spontaneous microtubule nucleation. We investigated whether GDP-tub could also be incorporated in growing microtubules at sub-critical tubulin concentrations, in seed-or centrosome-nucleated tubulin assembly conditions. In our experiments, critical concentration was 30 M. The μ incorporation of GDP-tub in microtubules elongating from EGS-seeded microtubules was monitored using filter assays as described above (

).

shows that microtubules assembled from tubulin solutions containing GDP-tub and GTP-tub in a 0.5/1 Fig. 1B D -Fig. 3(A B) proportion incorporated nearly 1 molecule of GDP-tub for 2 molecules of GTP-tub. We then tested whether GDP-tub incorporation occurred in microtubules nucleated on centrosomes by using immunofluorescence microscopy and tubulin variants as markers of GDP-tub and GTP-tub incorporation. Again, results showed apparently homogeneous qualitative incorporation of GDP-tub during assembly (Fig.

). 3C

These results indicate robust and direct incorporation of GDP-tub into seed or centrosome nucleated microtubules, at sub-critical tubulin concentrations.

Dynamic behavior of centrosome nucleated microtubules assembled in the presence of GDP-tub

Microtubules nucleated on centrosomes were visualized by standard video-DIC microscopy for direct tests of GDP-tub effects on individual polymer dynamics. Centrosomes were pre-adsorbed on the surface of a perfusion chamber. Samples containing GTP-tub (25 μ M) without or with increasing GDP-tub concentrations (from 5 to 25 M) were then injected into the chamber and observed under the μ microscope at 37 C. In these conditions, the average microtubule growth rate was unaffected by addition of up to 3.8 M of GDP-tub. In °μ the 3.8 M 12.5 M GDP-tub concentration range, the average growth rate decreased linearly with the GDP-tub concentration ( μ -μ Fig. 4A,B ). GDP-tub concentrations above 12.5 M inhibited microtubule nucleation on centrosomes.

μ

The catastrophe frequency increased when microtubules were assembled in the presence of GDP-tub concentration from 7.6 M and μ above ( ). It has been shown that, at sub-critical GTP-tub concentration, microtubule catastrophe rates increase with decreased Fig. 4A elongation rates ( , ). Here, the elongation rate was lower in the presence of GDP-tub compared to control. To know whether the 32 39 decrease in elongation rate accounted for increased catastrophe rate, we compared catastrophe frequencies versus growth rate in the absence ( , black symbols) and in the presence ( , empty symbols) of added GDP-tub. Plot analysis showed similar Fig. 4C Fig. 4C catastrophe rates at equal elongation rates, regardless of the presence or absence of added GDP-tub at the onset of assembly. These results indicate that, in our experimental conditions, GDP-tub addition increased the catastrophe frequency due to a decrease in the polymer growth rate, in the absence of detectable modification of intrinsic microtubule catastrophe properties.

Interestingly, the microtubule shrinkage rate decreased linearly with increasing GDP-tub concentration above 3.8 M, being circa 50

μ % inhibited at a 12.5 M GDP-tub concentration ( ). Microtubule shrinkage has been previously demonstrated to be a zero order μ Fig. 4D reaction, depending on intrinsic structural properties of microtubules, not on the composition of the soluble tubulin pool ( ). Accordingly, shrinkage rates were similar when chambers containing microtubules grown from centrosomes were perfused with PEM buffer alone or in the presence of added GDP-tub complexes ( ). The shrinkage rate was also in the same range during the spontaneous catastrophes Table 1 observed in control samples, containing only GTP-tub in the soluble pool (

). We also checked that, in the range of GTP Table 1 concentrations used in our study, the shrinkage rate of microtubules assembled at various GTP-tub concentrations was not significantly correlated with the microtubule growth rates, as previously shown ( ). 32

Collectively, our results indicate that direct GDP-tub incorporation into the microtubule wall occurs in both self-assembled microtubule and seed nucleated polymers. GDP-tub incorporation apparently induces a decrease of the microtubule growth rate and the microtubule shrinkage rate. The decrease of the shrinkage rate most probably results from a modification of the intrinsic stability properties of microtubules.

DISCUSSION

In this study, we show substantial GDP-tub incorporation into polymerizing microtubules with resulting impaired microtubule dynamic parameters. The main novelty of our study is the use of a simple system containing only GTP-tub and GDP-tub without any free nucleotide, allowing direct measurements and visualization of the incorporation of both GTP-tub and GDP-tub.

The possibility of a direct incorporation of GDP-tub into polymerizing microtubules has been subject of controversy. Based on turbidimetry measurements, previous studies have suggested GDP-tub incorporation in the microtubule wall ( , ), which has been 15 16 questioned in subsequent work (

). Additionally, direct GDP-tub incorporation seemed precluded according to recent structural 17 -19 studies indicating that the shape of GDP-tub does not fit the microtubule lattice ( ). Our data show a different picture, demonstrating that 7

robust and substantial direct GDP-tub incorporation into growing microtubules can be achieved. However, GDP-tub incorporation in microtubules occurs only if GTP-tub is also present at a concentration sufficient to support microtubule assembly, indicating obligatory co-incorporation of GTP-tub with GDP-tub. It has been suggested that microtubule assembly can involve the incorporation of tubulin oligomers ( ). Maybe, GDP-tub can be incorporated by hitch hiking polymerizing GTP-tub oligomers, although other models are 3

" "
possible. Such a co-incorporation of GDP-tub with GTP-tub fits with recent studies indicating that GTP-tub and GDP-tub may be similarly bent and that subunits straighten only after their incorporation into microtubules ( , ). 40 41

Our data indicate that neither initial microtubule nucleation nor initial microtubule elongation were sizeably affected by GDP-tubulin incorporation during microtubule assembly in our bulk microtubule assembly conditions. According to previous work ( ), at the high 22 tubulin concentrations used in bulk assembly conditions, the rate of microtubule elongation is limited during most of the assembly phase by intrinsic structural factors such as the speed of tube closure ( ). In this view our data could indicate that tube closure is not impaired 22 by co-incorporation of GDP-tub together with GTP-tub. In apparent contrast, in our study, the growth rate of centrosome nucleated microtubules decreased at increasing GDP-tub concentrations. However, at sub-critical tubulin concentrations, the availability of free tubulin dimers becomes rate limiting for microtubule growth ( ). Our data would then indicate that GDP-tubulin behaves as a 32 competitive inhibitor of GTP-tub, when the tubulin concentration becomes rate limiting.

In the present study, microtubules assembled from GDP-tub/GTP-tub mix displayed impaired disassembly behavior, involving a decrease of the individual polymer shrinkage rate. Previous works have established that hydrolysis of the tubulin bound nucleotide within the microtubule wall is required for subsequent microtubule disassembly ( , ). Thus, microtubule disassembly is dramatically 42 43 impaired when microtubules were assembled in the presence of a slowly hydrolysable analog of GTP, GMPCPP ( ). Obviously, in our 42 experiments, the incorporation of GDP-tub in the polymer wall is not followed by hydrolysis of the bound nucleotide. GDP may thus function as a natural non hydrolysable analog of GTP, with resulting impairment of the disassembly properties of microtubules assembled in the presence of GDP-tub.

Whether microtubules are assembled from GTP-tub alone or from GTP-tub/GDP-tub mix, they are ultimately composed of GDP-tub.

Yet, microtubules have different dynamic properties according to the composition of the starting tubulin solution. Such dynamic differences uncoupled to the bound nucleotide state in the microtubule wall provide a striking illustration of the recently proposed concept of microtubule structural plasticity ( ). 13

Microtubules assembled with GMPCPP show detectable modifications in their lattice organization ( ) and in the structure of their 44 oligomeric breakdown products ( ). We have not detected such modifications in polymers assembled in the presence of GDP-tub. 38

Maybe, in our conditions, lattice modifications are blurred by the mixed incorporation of GTP-and GDP-tub, compared to a homogeneous incorporation of GMPCPP-tubulin complexes in previous work. Alternatively, structural alterations in polymers assembled with GDP-tub may be truly undetectable at the level of electron-microscopy resolution.

The fundamental discovery that GTP hydrolysis in the microtubule wall is required for microtubule disassembly ( ), not for microtubule regulations. However, exchanging GTP for a slowly hydrolysable analog is the only known way to modulate the bound nucleotide hydrolysis in the microtubule wall during tubulin assembly and the existence of a naturally occurring non hydrolysable analog of GTP may look as a remote possibility. Our data indicate that GDP may represent such an analog and a regulation of microtubule disassembly through GDP-tub incorporation in cellular microtubules is an attractive possibility. In the absence of a structural signature allowing visual identification of polymers assembled with GDP-tub, tests of such a possibility may rely on the characterization and manipulation of putative regulatory systems which could modulate the GDP-tub/GTP-tub ratio at the end of cellular microtubules.

Footnotes:

* We thank Dr. Annie Andrieux and Dr. Christian Delphin for critical reading of the manuscript. This work was supported by grants from La Ligue Nationale Contre le Cancer (O.V., D.J.) and ANR-07-PCVI-0012 (I.A.).

The abbreviations used are

Detyr-tub : detyrosinated tubulin GDP-tub : GDP-tubulin complex

2 μ

 2 mM), DTT (5 mM), MgCl (25 mM), tyrosine (0.5 mM), KCl (100 mM) and tubulin tyrosine ligase (50 l/ml), at 30 C for 30 min. Sample 2 μ °was cooled down on ice and centrifuged at 200,000g at 4 C for 10 min to remove tubulin aggregates. The supernatant containing Tyr-tub °was subjected to gel filtration in Biogel P30 equilibrated in PEM buffer. To prepare Detyr-tub, purified tubulin (50 M) was incubated with carboxypeptidase A (2 g/ml) during 15 min at 30 C. Reaction was μ μ °stopped with 20 mM DTT addition. Aggregates were removed by centrifugation (200,000g, 10 min, 4 C). DTT was removed by gel °filtration through Biogel P30 chromatography.

Fig. 1 3 ] 2 -Fig. 2 3 H

 13223 Fig. 1 GDP-tub incorporation in spontaneous nucleated microtubules. (A D) GTP-tub (75 M) was assembled in the absence (A) or in the presence μ of GDP-tub (37.5 M (B), 75 M (C) or 120 M (D)). Either GTP-tub or GDP-tub was H -labelled and H incorporation in microtubules μ μ μ [ 3 ] [ 3 ] was determined to measure GTP-tub ( ) or GDP-tub ( ■ ) incorporation. Total polymeric tubulin ( ) was calculated by summation of the GDP-tub and GTP-tub curves. (E) Plot of (polymeric GDP-tub concentration/polymeric GTP-tub concentration) ratios at maximum assembly time versus (initial GDP-tub concentration/initial GTP-tub concentration) ratio. (F) Microtubule half-disassembly times. Plot of the microtubule half-disassembly time t versus (initial GDP-tub concentration/initial GTP-tub concentration) 1/2

Fig. 4

 4 Fig. 4 Microtubules grown from centrosomes: effects of GDP-tub on microtubular dynamics. (A) Dynamic parameters of microtubules polymerized from centrosomes and GTP-tub (25 M) in the presence of increasing GDP-tub concentrations. Standard deviations are represented in μ parentheses. n, number of events. (B) Microtubule growth rate versus GDP-tub concentration. (C) Catastrophe frequency versus growth rate of microtubules assembled with increasing concentrations (18, 22, 25 and 28 M) of GTP-tub (black squares) or with 25 M of GTP-tub in μ μ the presence of increasing GDP-tub concentrations (empty squares). (D) Microtubule shrinkage rate versus GDP-tub concentration.

  Drechsel DN , Hyman AA , Cobb MH , Kirschner MW . 1992 ; Mol Biol Cell . 3 : 1141 -1154 40 . Rice LM , Montabana EA , Agard DA . 2008 ; PNAS . 105 : 5378 -5383 41 . Buey RM , Diaz JF , Andreu JM . 2006 ; Biochemistry . 45 : 5933 -5938 42 . Hyman AA , Salser S , Drechsel DN , Unwin N , Mitchison TJ . 1992 ; Mol Biol Cell . 3 : 1155 -1167 43 . Caplow M , Ruhlen RL , Shanks J . 1994 ; J Cell Biol . 127 : 779 -788 44 . Hyman AA , Chr tien

	GMPCPP : guanylyl-( , )-methylene-diphosphonate α β
	GTP-tub : GTP-tubulin complex
	NDP kinase : nucleosidediphosphate kinase
	Tyr-tub : tyrosinated tubulin
	video-DIC : video-differential interference contrast
				é	D . 2004 ; Current Biol . 14 : 2086 -2095
	32 . Walker RA , O Brien '	ET , Pryer NK , Soboeiro MF , Voter WA , Erickson HP , Salmon ED . 1988 ; J Cell Biol . 107 : 1437 -1448
	33 . Nickerson JA , Wells WW . 1978 ; Biochem Biophys Res Comm . 85 : 820 -826
	34 . Gundersen GG , Kalnoski MH , Bulinski JC . 1984 ; Cell . 38 : 779 -789
	35 . Wade RH , Chr tien é	D , Job D . 1990 ; J Mol Biol . 212 : 775 -786
	36 . Chr tien é	D , Metoz F , Verde F , Karsenti E , Wade R . 1992 ; J Cell Biol . 117 : 1031 -1040
	37 . Tran PT , Joshi P , Salmon ED . 1997 ; J Struct Biol . 118 : 107 -118
	38 . M ller-Reichert ü	T , Chr tien é	D , Severin F , Hyman AA . 1998 ; Proc Natl Acad Sci USA . 95 : 3661 -3666
	39 .

é D , Arnal I , Wade RH . 1995 ; J Cell Biol . 128 : 117 -125

Table 1

 1 Microtubule shrinkage in different conditions. Microtubule were assembled from GTP-tub (65 M) at room temperature. Microtubule shrinkage rates were measured either during spontaneous μ catastrophes (Control) or after chamber perfusion with PEM buffer alone or PEM buffer containing 65 M GDP-tub as indicated. Standard deviations are represented in parentheses. n, number of events.

	μ

μ Control GTP-tub (65 M) μ

10

10

10

10

microtubule assembly as previously assumed, has been a substantial surprise for scientists in the microtubule field. This discovery indicated a strong link between GTP hydrolysis and microtubule stability, a remarkable microtubule feature that cells could use for

10

10

10

10

10


G Simonin
email: simonin@lirmm.fr

R Giroudeau

J C König
email: konig@lirmm.fr

Complexity and approximation for scheduling problem for coupled-tasks in presence of compatibility tasks

Keywords: coupled-tasks, compatibility graph, complexity, approximation

Introduction

In this paper, we present the problem of data acquisition according to compatibility constraints in a submarine torpedo. The torpedo is used in order to make cartography, topology studies, temperature measures and many other tasks in the water. The aim of this torpedo is to collect and process a set of data as soon as possible on a mono processor. In this way, it possesses few sensors, a mono processor and one type of tasks which must be scheduled: acquisition tasks. The acquisition tasks A = {A 1 , . . . , A n } can be considered as coupled-tasks introduced by [START_REF] Shapiro | Scheduling coupled-tasks[END_REF], indeed the torpedo sensors emit a wave which propagates in the water in order to collect the data. Each acquisition task A i has two sub-tasks, the first one a i sends an echo, the second one b i receives it. For a better reading, we will denote the processing time of each sub-task a i and b i . Between the sub-tasks, there is an incompressible idle time L i which represents the spread of the echo in the water.

At last, there exists compatibility constraints between acquisition tasks, due to the fact that some acquisition tasks cannot be processed at the same time as another task. In order to represent this constraint, a compatibility graph G c = (A, E c ) is introduced, where A is the set of coupled-tasks and E c represents the edges which link two coupled-tasks which can be executed simultaneously. In other words, at least one sub-task of a task A i may be executed during the idle time of another task A j (see example in Figure 1).

A 1 A 2 A 3 a 1 b 1 a 2 b 2 a 3 b 3 L 1 L 2 L 3 Compatibility graph Fig. 1. Example of compatibility constraints with ai = bi = 1, L1 = 3, L2 = L3 = 2
The aim of this problem is to produce a shortest schedule (i.e. to minimize the completion time of the last processed task) denoted by C max , in presence of compatibility contraints between acquisition tasks. In the scheduling theory, a problem is categorized by its machine environment, job characteristic and objective function. So using the notation scheme α|β|γ proposed by [START_REF] Graham | Optimization and approximation in deterministic sequencing and scheduling theory: a survey[END_REF], this problem will be defined as

1|(a i , L i , b i ), G c |C max .
Our work consists in measuring the impact of the compatibility graph on the complexity and approximation of scheduling problems with coupled-tasks on a mono processor. This paper is focusing on the limit between the polynomiality and the N P-completeness of our problem, when the compatibility constraint is introduced with different topologies.

The complexity of the scheduling problem, with coupled-tasks and a complete compatibility graph1 , has been investigated first by [START_REF] Orman | On the complexity of coupled-task scheduling[END_REF], [START_REF] Błażewicz | Scheduling of coupled tasks with unit processing times[END_REF][START_REF] Ahr | An exact algorithm for scheduling identical coupled tasks[END_REF]. Nevertheless, in this article we study a different problem in which coupled-tasks (or acquisition tasks) must respect a compatibility graph. By comparing the results of [START_REF] Orman | On the complexity of coupled-task scheduling[END_REF] and those obtained by relaxing the constraint of compatibility, we can measure the impact of compatibility constraint on this kind of problem. In such context, several results (Simonin et al. 2009) have been recently obtained according to coupled-tasks parameters. In the following, we consider the specific problem

Π = 1|(a i = L i = b i ), G c |C max
where all the coupled-tasks have a different length2 . Then, we study the complexity of Π according to the topology of the compatibility graph.

2 Polynomial-time approximation algorithm for Π

For general compatibility graph

We will develop a polynomial-time approximation algorithm for the coupled-tasks scheduling problem in presence of compatibility constraints where the duration of all the tasks is different. Therefore, for all A i and A j , i = j, we have

a i = L i = b i = a j = L j = b j .
With these hypotheses, the tasks are sorted according to their length L i (we have

L 1 < L 2 < . . . < L n ).
Theorem 1. A simple algorithm, which consists in scheduling the tasks consecutively, admits a ratio ρ equal to3 2 . Proof. The length of the schedule obtained by our algorithm is i 3L i , whereas the optimal length is greater or equal to i 2L i . Thus we obtain a ratio ρ ≤ 3 2 . The bound is trivial, and in the following we will propose a better ratio for a specific compatibility graph.

For a complete oriented bipartite graph

This section is devoted to the complexity and approximation in the case where the compatibility graph is a complete oriented 3 bipartite graph G c = (X, Y, E). Let X = {X 1 , . . . , X k } where x i is the length of task X i , and Y = {Y 1 , . . . , Y m } where y j is the length of task Y i . Without lost of generality, we have x 1 > . . . > x k ≥ 3y 1 > . . . > 3y m .

Computational complexity

Theorem 2. The problem

Π ′ : 1|a i = L i = b i , G c = complete oriented bipartite|C max is N P-complete.
Proof. In the following we consider a complete oriented bipartite graph G b c = (X, Y, E), with |X| = 2. The proof is established by a reduction from even partition problem [START_REF] Garey | Computers and Intractability, a Guide to the Theory of N P-Completeness[END_REF]:

Instance: A finite set K of n elements {e 1 , . . . , e n }, a bound B ∈ IN + and a size s(e i ) ∈ IN for each e i ∈ K such that each s(e i ) is even and such that ei∈K s(e i ) = 2B. Remark that the problem remains N P-complete even if ∀i = j, s(e i ) = s(e j ).

Question: Can K be partitioned into two disjoint sets K 1 , K 2 of K such that B = ei∈K1 s(e i ) = ei∈K2 s(e i ) ∈ IN? Let I * be an instance from even partition, we construct an instance I for our scheduling problem where C max = 18B +3 as follows:

We consider (n + 2) tasks denoted by, X 1 , X 2 ∈ X of length x 1 , x 2 , and Y 1 , . . . , Y n ∈ Y of length y 1 , . . . , y n such that:

-X 1 = (3B + 1, 3B + 1, 3B + 1) and X 2 = (3B, 3B, 3B)

-Y i = (a i = s(e i ), L i = s(e i ), b i = s(e i )), ∀i ∈ {1, . . . , n}
We add edges between (X i , Y j ), ∀i ∈ {1, 2}, ∀j ∈ {1, . . . , n}. Obviously, this transformation can be computed in polynomial time.

• Assume that K can be partitioned into two disjoint sets K 1 , K 2 of K such that B = ei∈K1 s(e i ) = ei∈K2 s(e i ) ∈ IN. The task X 1 (resp. X 2 ) is executed at t = 0 (resp. at t = 9B + 3). The task Y i , corresponding to the element e i ∈ K 1 (resp. e i ∈ K 2 ) of length s(e i ), is scheduled as soon as possible in the idle slot of the task X 1 (resp. X 2 ).

• Reciprocally, we suppose that the length of the schedule is C max = 18B + 3. The two tasks X 1 and X 2 cannot be executed into each other. It is clear that the sum of durations of these two tasks is 18B + 3. All the other tasks must be executed in idle slots of length 3B or 3B + 1, which implies the existence of a partition.

Polynomial-time Approximation Algorithm

We propose the following algorithm: initially the X-tasks, |X| = k ≥ 1, and the Y -tasks, |Y | = m ≥ 1, are sorted in non-increasing order with respect to their duration, separately. The tasks are executed in the input order according to the following rule: X-tasks are executed consecutively. A Y -task is assigned to the first used X-task where there is enough available space to include it; if no X-tasks can contain the Y -tasks, they are assigned after the last executed task in the schedule.

This algorithm is similar to the classical Best-Fit algorithm. The sets X and Y are sorted in O(nlog(n)), where |X| + |Y | = n. Definition 1. The set of Y -tasks, which are not executed in the X-tasks, is called remainder R.

Theorem 3. The previous algorithm admits a relative performance bounded by 7/6 for the problem

Π ′ : 1|(a i = L i = b i = x i ), G c = complete oriented bipartite|C max .
Proof. We suppose that there exists a remainder R (otherwise there would be an optimal solution). We consider

C 1 = m i=1 3y i and C 2 = k i=1 x i .
Lemma 1. The size of the slot (i.e. the inactivity time of a couple-task) x i , ∀X i ∈ X, is at least greater than xi 2 for an execution of the Y -tasks if there exists a remainder. Proof. By contradiction, let X i be the first task such that the used time in its slot is less than xi 2 , and let Y j be the first task of the remainder. We have y j > xi 6 , otherwise the task Y j may be executed into the task X i . In such context, there is no task from the Y -tasks which are executed into X i when Y j is considered, because the Y -tasks are sorted in nonincresang order, so y j > xi 3 (in contradiction with the fact that the graph G c is an oriented bipartite complete graph, and that any Y -task can be scheduled in the X-tasks).

Thus, we have

0 ≤ R ≤ C 1 -C2
2 , and one upper bound is:

C h max ≤ Remainder + sequential duration of the X-tasks = C 1 -C2 2 + 3C 2 = 5C2
2 + C 1 . Two cases must be considered:

-If C 1 ≤ C 2 then C opt max ≥ 3C 2 .
Indeed, in the best case the Y -tasks can be executed in the idle slots of the X-tasks. Thus, the ratio is ρ =

C h max C opt max ≤ 5C 2 2 +C1 3C2 ≤ 7 6 . -If C 1 > C 2 then C opt max ≥ 2C 2 + C 1 .
Indeed, in the best case all the inactivity times of the X-tasks are used. The ratio is ρ =

C h max C opt max ≤ 5C 2 2 +C1 2C2+C1 ≤ 5C 2 2 +y 2C2+y ≤ 7 6 , ∀y ≤ C 1 .
The worst case analysis of the relative performance is not rigorous enough. We conjecture that the ratio is 8 7 . In such context, we may exhibit an instance for which the bound is tight. We consider the four tasks

X 1 ∈ X, Y 1 ∈ Y , Y 2 ∈ Y et Y 3 ∈ Y with following char- acteristics: X 1 = (M, M, M ), Y 1 = ( M 6 + 2ǫ, M 6 + 2ǫ, M 6 + 2ǫ), Y 2 = ( M 6 + ǫ, M 6 + ǫ, M 6 + ǫ), Y 3 = ( M 6 -ǫ, M 6 -ǫ, M 6 -ǫ) with M ∈ IN * .
Using the previous algorithm, we obtain a schedule of length C h max = 4M. An optimal solution consists in executing the tasks Y 2 and Y 3 within the idle time of the task X 1 . Therefore, we obtain C opt max = 7M 2 . The ratio is equal to ρ = 8 7 .

Conclusion

We investigate a particular coupled-tasks scheduling problem 1|a i = L i = b i , G c |C max in presence of a compatibility graph with regard to the complexity and appproximation. We design a 3 2 -approximation algorithm for the general case. We also establish the N Pcompleteness for the specific case where there is a bipartite compatibility graph. In such context, we propose a 7 6 -approximation algorithm and the bound is tight. For further research, it would be interesting to propose a polynomial-time approximation algorithm with a non-trivial ratio for non complete oriented bipartite graph.

Note that the lack of compatibility graph is equivalent to a fully connected graph. In this way, all the tasks may be compatible with each other.

The problem with complete compatibily graph is N P-complete (see[START_REF] Orman | On the complexity of coupled-task scheduling[END_REF], so the relaxed version is obviously N P-complete).

the graph is oriented because of the Y -tasks can be scheduled in the X-tasks.