Vlad Stefan Barbu
email: barbu@univ-rouen.fr

Nicolas Vergne
email: nicolas.vergne@univ-rouen.fr

Reliability and survival analysis for drifting Markov models: modeling and estimation

Keywords: Drifting Markov chains, Multi-state systems, Reliability theory, Survival analysis, Estimation, Asymptotic properties Mathematics Subject Classification (2000) 60J10, 60K15, 90B25

HAL is

Note that a hypothesis used in many mathematical models built as modeling tools for real applications is the homogeneity with respect to time.

On the one hand, it is clear that in many applications this homogeneity is inappropriate, at least for two reasons: first, there is often a non-homogeneity (or heterogeneity) intrinsic to the evolution of the systems; second, we can have interactions between the system and the environment, that can generate a non-homogeneity in the functioning of the system.

On the other hand, from a practical point of view, considering general non-homogeneous processes could be inappropriate, as the estimation procedures could be not tractable in practice. A possible solution is to consider a non-homogeneity that is "smooth", controlled, of a known shape. An example of this type in a Markov framework consists in the so called drifting Markov chains, introduced in [START_REF] Vergne | Drifting Markov models with polynomial drift and applications to DNA sequences[END_REF]. For these processes, the Markov transition matrix is a linear (polynomial) function/mixture of two (several) Markov transition matrices. In other words, instead of fitting a transition matrix on a whole sequence of observations (homogeneous Markov model) or different transition matrices on different homogeneous parts of the observation (like we do in hidden Markov models), we allow the transition matrix to vary (to drift) from the beginning to the end of the observed phenomenon. Consequently, at each position, we obtain a different transition matrix, but with a controlled, known shape. These models are thus constrained non-homogeneous Markov models. Thus we obtain the desired "smooth" non-homogeneity. In this paper, we focus on polynomial drift ans associated reliability measures and estimation techniques. To the best of our knowledge, the use of such models in reliability/survival studies is a completely new approach.

Note that the type of models we propose is appropriate for modeling phenomena that evolve somehow smooth along time. On the one hand, it is worth to stress that this evolution is not only in one way, for instance aging evolution. The polynomial drift allows the models to take several possible shapes. On the other hand, if one needs piecewise homogeneity with abrupt transitions between the homogeneous zone, one possible setting is furnished by hidden Markov or hidden semi-Markov processes (see, e.g., [26], [START_REF] Cappé | Inference in hidden Markov models[END_REF], [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications-their use in reliability and DNA analysis[END_REF] for general references on hidden Markov and semi-Markov models).

A remark needs to be done here. It is obvious that the family of stochastic models we propose in this article is a particular case of a non-homogeneous Markov chain. Consequently, using a general non-homogeneous Markov chain could be seen as a better modeling approach (and really could be in some situations). Nonetheless, note that the number of parameters of a (fully) non-homogeneous Markov chain could be huge as compared to the number of parameters of our models: two (for linear drift) or several (i.e., few) Markov transition matrices (for polynomial drift). Moreover, in order to be able to estimate them all, one needs an extremely huge number of copies of same type of trajectories (at every time k, we need an important or at least reasonable number of observations of transitions from any state to any another state).

The work developed along this paper belongs to MSS reliability methods, that have received a substantial amount of attention in the last decades. For theoretical advances and significant applications in MSS reliability theory in recent years, the reader is referred to [START_REF] Levitin | The universal generating function in reliability analysis and optimization[END_REF], [START_REF] Lisnianski | Multi-state system reliability analysis and optimization for engineers and industrial managers[END_REF], [START_REF] Lisnianski | Multi-state system reliability: assessment, optimization and applications[END_REF] and [START_REF] Natvig | Multistate systems reliability. Theory with applications[END_REF].

Typical tools for investigating the evolution and performance of such systems are the Markov and semi-Markov processes. For general references on Markov processes and applications, as well as associated reliability topics, in continuous or discrete time, with corresponding estimation topics, one can see, e. g., [START_REF] Sadek | Asymptotic properties for maximum likelihood estimators for reliability and failure rates of Markov chains[END_REF], [START_REF] Sadek | Nonparametric estimation of reliability and survival function for continuous-time finite Markov processes[END_REF], [START_REF] Chiquet | Estimating stochastic dynamical systems driven by a continuous-time jump Markov process[END_REF], [START_REF] Vasiliadis | On the distributions of the state sizes of closed continuous time homogeneous Markov systems[END_REF], [START_REF] Vassiliou | G-inhomogeneous Markov systems of high order[END_REF], [START_REF] D'amico | Discrete time homogeneous Markov processes for the study of the basic risk processes[END_REF], [START_REF] Ruiz-Castro | Markov counting and reward processes for analysing the performance of a complex system subject to random inspections[END_REF]. Semi-Markov processes generalize renewal processes as well as Markov jump processes by allowing general distributions for sojourn times; this is the reason why the semi-Markov processes are more adapted for applications in general and for reliability studies in particular. One can see [18] for an introduction to continuous-time semi-Markov processes and [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications-their use in reliability and DNA analysis[END_REF] for a presentation of the corresponding discrete-time framework.

For general references on continuous-time semi Markov systems, associated applications and reliability topics, one can see for example [START_REF] Limnios | Empirical estimators of reliability and related functions for semi-Markov systems[END_REF], [START_REF] Limnios | Nonparametric estimation of some important indicators in reliability for semi-Markov processes[END_REF], [29], [START_REF] Alvarez | Estimation in stationary Markov renewal processes, with application to earthquake forecasting in Turkey[END_REF], [27], [START_REF] Votsi | Estimation of the expected number of earthquake occurrences based on semi-Markov models[END_REF], [START_REF] Barbu | Semi-Markov modelling for multi-state systems[END_REF], [START_REF] Hou | On the existence and uniqueness of solution of MRE and applications[END_REF], [START_REF] Koutras | Optimization of the dependability and performance measures of a generic model for multistate deteriorating systems under maintenance[END_REF], [START_REF] Silvestrov | Reward algorithms for semi-Markov processes[END_REF]. Discrete-time semi-Markov modeling and estimation, and reliability computation and estimation can be found in [START_REF] Barbu | Discrete time semi-Markov processes for reliability and survival analysis-a nonparametric estimation approach[END_REF], [START_REF] Barbu | Nonparametric estimation for failure rate functions of discrete time semi-Markov processes[END_REF], [START_REF] Barbu | Reliability of semi-Markov systems in discrete time: modeling and estimation[END_REF], [START_REF] Chryssaphinou | Multi-state reliability systems under discrete time semi-Markovian hypothesis[END_REF], [START_REF] D'amico | Performance analysis of second order semi-Markov chains: an application to wind energy production[END_REF]. A unified semi-Markov framework is described in [START_REF] Limnios | Reliability measures of semi-Markov systems with general state space[END_REF]. This paper is structured as follows. We first introduce several types of drifting Markov models in the next section and present basic properties. We also obtain explicit expressions for reliability/survival indicators of drifting Markov models, like reliability, availability, maintainability and failure rates. Then, under different statistical settings, in Section 3 we estimate the parameters of the models, obtain plug-in estimators of the associated reliability/survival indicators and investigate the consistency of the estimators. The quality of the proposed estimators is also illustrated in Section 4, where numerical experiments are carried out.

Drifting Markov models

In this section we introduce several types of drifting Markov models, namely the linear drifting Markov chain of order 1 (Definition 1) and the polynomial drifting Markov chain of order k (Definition 2), and present basic properties. Then we obtain explicit expressions for associated reliability/survival indicators, like reliability, availability, maintainability and failure rates; note that two distinct definitions of failure rates are considered.

Definitions and notation

Let us consider a random system with finite state space E = {1, . . . , s}, s < ∞, and let (Ω, A, P) be a probability space.

All along this paper we assume that the time evolution of the system is governed by discrete-time stochastic processes defined on (Ω, A, P) with values in E. Let Π 0 = (Π 0 (u, v)) u,v∈E and Π 1 = (Π 1 (u, v)) u,v∈E be two Markov transition matrices (of order 1) over a state space E.

Definition 1 (linear drifting Markov chain of order 1 and of length n) A sequence X 0 , X 1 , . . . , Xn with state space E = {1, 2, . . . , s} is said to be a linear drifting Markov chain (of order 1) of length n between the Markov transition matrices Π 0 and Π 1 if the distribution of X t , t = 1, . . . , n, is defined by

P(X t = v | X t-1 = u, X t-2 , . . .) = Π t n (u, v), u, v ∈ E, (1)
where

Π t n (u, v) = 1 - t n Π 0 (u, v) + t n Π 1 (u, v), u, v ∈ E. (2)
Let us denote by α = (α(1), . . . , α(s)) the initial distribution of the chain, that is the distribution of X 0 , α(u) = P(X 0 = u) for any state u ∈ E.

The following result, although straightforward, will be useful in the sequel.

Lemma 1 For a linear drifting Markov chain and

k 1 ≤ k 2 , k 1 , k 2 ∈ N, we have P(X k2 = j | X k1-1 = i) =   k2 l=k1 Π l n   (i, j).
Proof.

The proof is straightforward and we omit it. In a most general case, it is given in Lemma 2

2

Note that one could consider a k order linear drifting Markov chain, k ∈ N * , by allowing a dependence of order k in Equation ( 1) and letting Π 0 et Π 1 be the transition probabilities of Markov chains of order k. Note also that the linear drifting can be generalized to polynomial drifting of a certain degree d. See [START_REF] Vergne | Drifting Markov models with polynomial drift and applications to DNA sequences[END_REF] for more details on these points.

Let us generalize linear drifting Markov model of order 1 to polynomial drifting Markov model of order k and degree d. Let

Π i d = (Π i d (u 1 , . . . , u k , v))
P(X t = v | X t-1 = u k , X t-2 = u k-1 , . . .) = Π t n (u 1 , . . . , u k , v), u 1 , . . . , u k , v ∈ E, (3)
where

Π t n (u 1 , . . . , u k , v) = d i=0 A i (t)Π i d (u 1 , . . . , u k , v), u 1 , . . . , u k , v ∈ E, (4)
Indeed v∈E Π t n (u 1 , . . . , u k , v) = d i=0 A i (t) = A(t)
where A is a polynomial of degree d equal to one in (d + 1) points; then

A is a constant polynomial equal to one.

Example 1 (Polynomial drifting Markov model of degree 2) For a polynomial drifting Markov model of degree 2, we have

Π t n (u, v) = A 0 (t)Π 0 (u, v) + A 1 (t)Π 1 2 (u, v) + A 2 (t)Π 1 (u, v),
with A 0 , A 1 , A 2 polynomials of degree 2 and Π 0 , Π 1 2 , Π 1 parameter matrices.

We want to find A 0 , A 1 , A 2 such that Π t n be a stochastic matrix and :

-for t = 0, Π t n = Π 0 , i.e. A 0 (0) = 1, A 1 (0) = 0, A 2 (0) = 0; -for t = n/2, Π t n = Π 1 2 , i.e. A 0 (n/2) = 0, A 1 (n/2) = 1, A 2 (n/2) = 0; -for t = n, Π t n = Π 1 , i.e. A 0 (n) = 0, A 1 (n) = 0, A 2 (n) = 1.
Solving these three systems we obtain

Π t n (u, v) = 2 t 2 n 2 -3 t n + 1 Π 0 (u, v) + -4 t 2 n 2 + 4 t n Π 1 2 (u, v) + 2 t 2 n 2 - t n Π 1 (u, v).
Obviously, when dealing with a k order Markov chain, we can always transform it into a Markov chain of order 1, as it is stated in the following remark.

Remark 3 (From a Markov chain of order k to a Markov chain of order 1) As it is well known, it is always possible to transform a Markov chain of order k into a Markov chain of order 1. We only have to change the state space: for an order k Markov chain over a state space E with transition π, we define an order 1 Markov chain over a state space E = E k with transition π such as for u ∈ E, v, v ∈ E k-1 and w ∈ E, π (uv; vw) = π(uv; w) and π (uv

; v w) = 0 if v = v .
The following result is the analogous of Lemma 1 for a polynomial drifting Markov chain of order k instead of a linear drifting Markov chain of order 1.

Lemma 2 For a polynomial drifting Markov chain of order k and

k 1 , k 2 ∈ N, such that k 2 -k 1 ≥ k -1, we have P(X k2 = j | X k1-1 = i k , . . . , X k1-k = i 1 ) = v1,...,v k-1 ∈E   k2 t=k1 Π t n   ((i 1 , . . . , i k ); (v 1 , . . . , v k-1 , j))
with Π t n the corresponding transition matrix of the Markov chain of order 1 in the transformed alphabet E (see Remark 3).

Proof.

We immediately have

P(X k2 = j | X k1-1 = i k , . . . , X k1-k = i 1 ) = i k 1 ,i k 1 +1 ,...,i k 2 -1 P(X k2 = j | X k2-1 = i k2-1 , . . . , X k2-k = i k2-k ) ×P(X k2-1 = i k2-1 | X k2-2 = i k2-2 , . . . , X k2-k-1 = i k2-k-1 ) × . . . × P(X k1 = i k1 | X k1-1 = i k , . . . , X k1-k = i 1 ) = i k 1 ,i k 1 +1 ,...,i k 2 -1 Π k 2 n (i k2-k , . . . , i k2-1 , j)Π k 2 -1 n (i k2-k-1 , . . . , i k2-2 , i k2-1 ) . . . Π k 1 n (i 1 , . . . , i k , i k1 ) = v1,...,v k-1 ∈E   k2 t=k1 Π t n   ((i 1 , . . . , i k ); (v 1 , . . . , v k-1 , j))
where in the last line, we take into account all the possibilities for the k states before state j on the position k 2 .

2

We would like to make some comments on the interest in practice and applicability of the drifting Markov models introduced in Definitions 1 and 2.

Remark 4 (practical interest of drifting Markov models) Note that the family of models introduced in this paper is fairly general, with a wide potential for real life applications. First, note that the classical Markov models of arbitrary order are particular cases of drifting Markov models. Second, in particular applications where the time evolution of the system is modelled by a Markov chain, it could make sense to assume that the evolution is different at the beginning of the lifetime, as compared to the end of the lifetime of the system. In such a situation, a linear drifting should be appropriate.

Another typical situation is when the time evolution of a system could be split in three parts, corresponding to early failure/ stable failure/wear-out failure periods, respectively. In this case, a polynomial drifting Markov model of degree 2 could be appropriate. Surely, this decomposition in different epochs of the lifetime of a system can be considered not only for 2 or 3 epochs, but also for a different number, corresponding to different real situations depending on the application.

As mentioned in the Introduction, the drifting Markov models are particular non-homogeneous Markov models, with a controlled non-homogeneity. This is in fact one core point of the practical interest of drifting Markov models: it comes from the fact that we can adapt a parsimonious non-homogeneous model to a real application, with a reasonable number of parameters (as compared to a fully non-homogeneous model).

Reliability and survival measures

In order to undertake a reliability analysis of a repairable system modeled by a drifting Markov chain, let us assume that the state space of the system is partitioned into working and failure states, E = U ∪ D, with U ∩ D = ∅, where U = {1, . . . , s 1 } are the working states and D = {s 1 + 1, . . . , s} the failure states of the system. Note that the states of U can be seen as representing different operating modes or performance levels of the system (with state "1" associated with the nominal performance of the system), while the states of D can represent different failure modes of the systems (with state "s" associated with the total failure).

According to this partition of the state space we partition any matrix of vector we are working with. For instance, we have

U D Π 0 =   Π U U 0 Π DU 0 Π U D 0 Π DD 0   U D , U D α = α U α D .
For a system starting to function at time l = 0 (that is α D = 0) the reliability at time l ∈ N is the probability that the system has functioned without failure in the period [0, l].

Denoting by T D hitting time of subset D, called the lifetime of the system, i.e.,

T D := inf{t ∈ N; X t ∈ D}, with inf ∅ := ∞, the reliability at time l ∈ N is R(l) := P(T D > l) = P(X t ∈ U, t = 0, . . . , l).
The following results will present the expressions of reliability measures for the "simple" case of linear drifting Markov chains of order 1 (Propositions 1, 3, 5 and 7) and for the general case of polynomial drifting Markov chains of order k (Propositions 2, 4, 6).

Proposition 1 For a linear drifting Markov chain of order 1 (X t ) 0≤t≤n , the reliability at time l, l ∈ N, is given by

R(l) = α U l t=1 1 - t n Π U U 0 + t n Π U U 1 1 U , where 1 U = (1, • • • , 1 s1 ) . (5)
Proof.

We have

R(l) = P(X t ∈ U, t = 0, . . . , l) = i,j∈U P(X l = j, X l-1 ∈ U, . . . , X 1 ∈ U | X 0 = i)P(X 0 = i) = i,j∈U α(i) l t=1 Π U U t n (i, j),
where

Π U U t n := 1 -t n Π U U 0 + t n Π U U 1
and we have used Lemma 1, with

Π U U t n , Π U U 1 and Π U U 0 instead of Π t n
, Π 1 and Π 0 respectively.

2 Proposition 2 For a polynomial drifting Markov chain of order k (X t ) 0≤t≤n , the reliability at time l, l ∈ N, is given by

R(l) = α U j∈U   v1,...,v k-1 ∈U l t=k Π U U t n   (v 1 , . . . , v k-1 , j). (6)
Proof.

The proof is analogous to the proof of Proposition 1, using the passage from a Markov chain of order k to a Markov chain of order 1 (cf. Remark 3). We have

R(l) = i0,...,i k-1 ,j∈U P(X l = j, X l-1 ∈ U, . . . , X k ∈ U | X k-1 = i k-1 , . . . , X 0 = i 0 )P(X k-1 = i k-1 , . . . , X 0 = i 0 ) = i0,...,i k-1 ,j∈U α(i 0 , . . . , i k-1 )   v1,...,v k-1 ∈U l t=k Π U U t n ((i 0 , . . . , i k-1 ); (v 1 , . . . , v k-1 , j))   ,
where

Π U U t n := 1 -t n Π U U 0 + t n Π U U 1 and we have used Lemma 2, with Π U U t n , Π U U 1 and Π U U 0 instead of Π t n
, Π 1 and Π 0 respectively (that is to say that we considered the restriction of the initial k order Markov chain to the subspace U , and, consequently, the associate restriction of the extended Markov chain of order 1).

Let us now consider the pointwise (or instantaneous) availability; for any l ∈ N, the availability at time l it is defined as the probability that the system is in a working state at time l, independently of the fact that the system worked or not during the time interval [0; l), A(l) := P(X l ∈ U ).

Proposition 3 For a linear drifting Markov chain of order 1 (X t ) 0≤t≤n , the pointwise (or instantaneous) availability at time l, l ∈ N, is given by

A(l) = α l t=1 1 - t n Π 0 + t n Π 1 1 E,U , where 1 E,U = (1, • • • , 1 s1 , 0, • • • , 0 s-s1
) .

Proof.

As in the proof of the expression of the reliability, we immediately have

A(l) = i∈E,j∈U P(X l = j | X 0 = i)P(X 0 = i) = i∈E,j∈U α(i) l t=1 Π t n (i, j).
2 Proposition 4 For a polynomial drifting Markov chain of order k (X t ) 0≤t≤n , the pointwise (or instantaneous) availability at time l, l ∈ N, is given by

A(l) = α j∈U   v1,...,v k-1 ∈E l t=k Π t n   (v 1 , . . . , v k-1 , j). (8)
Proof.

The proof is similar to the proof of Proposition 2.

2

A remark needs to be done here; for a similar remark in a discrete-time semi-Markov framework, see [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications-their use in reliability and DNA analysis[END_REF].

Remark 5 When computing reliability measures, one can consider either the case of ergodic (i.e., repairable) systems, or the case of non-ergodic (i.e., non-repairable) systems. For non-repairable systems, the set D is absorbing; consequently, in this case the availability and reliability coincide and are given by Formula (6). Obviously, for repairable systems availability and reliability do not coincide, and one has to consider the different formulas presented above.

To conclude, although in this paper we consider ergodic (repairable) systems, explicit forms of the reliability indicators can be immediately obtained even in the non-ergodic (non-repairable) case.

Let us assume that we deal with a repairable system for which the failure occurs at time l = 0 (in other words, we take α U = 0). Its maintainability at time k ∈ N is the probability that the system is repaired up to time l, given that it has failed at time l = 0. Let us denote by T U the first hitting time of subset U, called the duration of repair or repair time, that is,

T U = inf{t ∈ N; X t ∈ U }. Thus, the maintainability at time l ∈ N is defined by M (l) := P(T U ≤ l) = 1 -P(T U > l) = 1 -P(X t ∈ D, t = 0, . . . , l).
Proposition 5 For a linear drifting Markov chain of order 1 (X t ) 0≤t≤n , the maintainability at time l, l ∈ N, is given by

M (l) = 1 -α D l t=1 1 - t n Π DD 0 + t n Π DD 1 1 D , where 1 D = (1, • • • , 1 s-s1
) .

Proof.

The proof is similar to the one of Proposition 1, exchanging the places of U and D.

2 Proposition 6 For a polynomial drifting Markov chain of order k (X t ) 0≤t≤n , the maintainability at time l, l ∈ N, is given by

M (l) = 1 -α D j∈D   v1,...,v k-1 ∈D l t=k Π DD t n   (v 1 , . . . , v k-1 , j). ( 10
)
Proof.

The proof is similar to the proof of Proposition 2, exchanging the places of U and D.

2

In the sequel we will be interested in the failure rate of the system. To this purpose, we will consider two different definitions of the failure rate: the first one is the classical failure rate introduced in [START_REF] Barlow | Properties of probability distributions with monotone hazard rate[END_REF] and usually considered for continuous time systems; this one will be called the BMP-failure rate and will be denoted by λ(l), l ∈ N. The second one, proposed in [START_REF] Roy | Classification of discrete lives[END_REF], is a failure rate adapted to the work in discrete time; we will call it the RG-failure rate and denote it by r(l), l ∈ N.

The interested reader could see also [START_REF] Xie | Redefining failure rate function for discrete distributions[END_REF] for detailed justification of the use of the RG-failure rate in discrete time, instead of the classical one. The use of this failure rate for discrete time Markov chains can be found in [START_REF] Sadek | Asymptotic properties for maximum likelihood estimators for reliability and failure rates of Markov chains[END_REF], while for discrete time semi-Markov chains can be found in [START_REF] Barbu | Nonparametric estimation for failure rate functions of discrete time semi-Markov processes[END_REF] and [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications-their use in reliability and DNA analysis[END_REF].

Definition 3 (failure rates; cf. [START_REF] Barlow | Properties of probability distributions with monotone hazard rate[END_REF], [START_REF] Roy | Classification of discrete lives[END_REF], [START_REF] Xie | Redefining failure rate function for discrete distributions[END_REF]) Let us consider a discrete time system starting to work at time l = 0 (that is α D = 0).

1. The BMP-failure rate at time l ∈ N is the conditional probability that the failure of the system occurs at time l, given that the system has worked until time l -1, i.e., λ(l)

:= P(T D = l | T D ≥ l) =    1 -R(l) R(l-1) , if R(l -1) = 0, 0, otherwise. (11)
2. The RG-failure rate at time l ∈ N is defined by

r(l) :=    -ln R(l) R(l-1) , if l ≥ 1, -ln R(0), if l = 0, (12)
for R(l) = 0. If R(l) = 0, we set r(l) := 0.

Note that the two failure rates are related by r(l) = -ln(1 -λ(l)), l ∈ N.

Taking into account the expression of the reliability of a drifting Markov chain given in Proposition 1, we immediately obtain the expression of the failure rates for a linear drifting Markov chain of order 1. The general case can be immediately written using Proposition 2.

Proposition 7 For a linear drifting Markov chain of order 1 (X t ) 0≤t≤n , the BMP-failure rate λ(l) and the RG-failure rate r(l) are given by:

λ(l) =    1 - α U l t=1 ((1-t n )Π U U 0 + t n Π U U 1 ) 1 U α U l-1 t=1 ((1-t n )Π U U 0 + t n Π U U 1 ) 1 U , if R(l -1) = 0, 0, otherwise,
where

1 U = (1, • • • , 1 s1 ) . ( 13
)
r(l) =    -ln α U l t=1 ((1-t n )Π U U 0 + t n Π U U 1 ) 1 U α U l-1 t=1 ((1-t n )Π U U 0 + t n Π U U 1 ) 1 U , if l ≥ 1, -ln R(0), if l = 0, (14)
for R(l) = 0. If R(l) = 0, we set r(l) := 0.

The purpose of this section is twofold: first, continuing the direction developed in [START_REF] Vergne | Drifting Markov models with polynomial drift and applications to DNA sequences[END_REF], we will consider different types of data for which the estimators of the characteristics of a drifting Markov chain will be derived. Second, we will estimate the associated reliability indicators.

For estimation purposes we first consider two types of data: the first one consists in one sample path of a linear drifting Markov chain; the second one consists in several sample paths of the same fixed length. Then we consider various types of data that can occur in practice and we adapt the estimation procedures that we have previously obtained in order to get estimators of the drifting Markov models in these cases. These different types of data arise from three ingredients: an observed sample path varying between two Markov transition matrices could be complete or incomplete; we can have at our disposal either one or several sample paths; several observed sample paths could be of the same length or of different length. The possible types of data that we take into account are the following:

A. One can observe a sample path, that will be denoted by H(m, n) := (X 0 , X 1 , . . . , Xm), where m denotes the length of the sample path and n the length of the drifting Markov chain. Two cases can be considered:

A1 : m = n (a complete sample path);

A2 : m ≤ n (an incomplete sample path).

B. One can also observe H i.i.d. sample paths, H i (m i , n i ), i = 1, . . . , H. Four cases can be considered here: B4 : m i ≤ n i for all i = 1, . . . , H (incomplete sample paths of drifting Markov chains of different lengths).

B1 : m i = n i = n for all i = 1, . . . ,
Notation 1 We introduce the following notations:

P 1 (H, m, n) = m t=1 H h=1 1 {X h t-1 =u} t n 2 , P 2 (H, m, n) = m t=1 H h=1 1 {X h t-1 =u,X h t =v} 1 - t n , P 3 (H, m, n) = m t=1 H h=1 1 {X h t-1 =u} 1 - t n t n , P 4 (H, m, n) = m t=1 H h=1 1 {X h t-1 =u,X h t =v} t n , P 5 (H, m, n) = m t=1 H h=1 1 {X h t-1 =u} 1 - t n 2
, with m = (m 1 , . . . , m H ) and n = (n 1 , . . . , n H ).

Estimation of drifting Markov models: the usual types of data

From [START_REF] Vergne | Drifting Markov models with polynomial drift and applications to DNA sequences[END_REF], we have the following result.

Proposition 8 For (X 0 , X 1 , . . . , Xn) a sample path of a linear drifting Markov chain between two Markov transition matrices (of order 1) Π 0 and Π 1 , for any states u, v ∈ E, the least square estimators of Π 0 (u, v) and Π 1 (u, v) are given by:

Π 0;n (u, v) = P 1 (1, n, n)P 2 (1, n, n) -P 3 (1, n, n)P 4 (1, n, n) P 5 (1, n, n)P 1 (1, n, n) -P 3 (1, n, n) 2 (15)
Π 1;n (u, v) = P 5 (1, n, n)P 4 (1, n, n) -P 3 (1, n, n)P 2 (1, n, n) P 5 (1, n, n)P 1 (1, n, n) -P 3 (1, n, n) 2 . ( 16
)
Consequently, from Equation ( 2), for any u, v ∈ E and t = 1, . . . , n, the estimator of Π t n (u, v) is given by

Π t n (u, v) = 1 - t n Π 0;n (u, v) + t n Π 1;n (u, v), (17)
with Π 0;n (u, v) and Π 1;n (u, v) given in ( 15) and ( 16), respectively.

Second, let us now consider H i.i.d. sample paths of a linear drifting Markov chain (X k ) 0≤k≤n , H i (m, n), i = 1, . . . , H.

That is, we assume that the sample paths have the same length n and that the linear drift within the H sample paths takes place between the same two Markov transition matrices (of order 1) Π 0 and Π 1 (framework B1).

Proposition 9 Under this setting, starting with H sample paths of a linear drifting Markov chain between two Markov transition matrices (of order 1) Π 0 and Π 1 , for any states u, v ∈ E, the estimators of Π 0 (u, v) and Π 1 (u, v) are given by:

Π 0;(n,H) (u, v) = P 1 (H, n, n)P 2 (H, n, n) -P 3 (H, n, n)P 4 (H, n, n) P 5 (H, n, n)P 1 (H, n, n) -P 3 (H, n, n) 2 , ( 18
)
Π 1;(n,H) (u, v) = P 5 (H, n, n)P 4 (H, n, n) -P 3 (H, n, n)P 2 (H, n, n) P 5 (H, n, n)P 1 (H, n, n) -P 3 (H, n, n) 2 , ( 19
)
with n = (n 1 , . . . , n H ) = (n, . . . , n).

Proof.

As in [START_REF] Vergne | Drifting Markov models with polynomial drift and applications to DNA sequences[END_REF], we minimize a quadratic form of the different parameters which is the sum of "prediction errors". For each sequence h, at each position t in the sequence, knowing X t-1 = u, we want Π t n (u, v) to be the nearest possible to 1 if

X t = v or the nearest possible to 0 if X t = v.
We minimize the sum of squares of the errors

H h=1 1 {X h t-1 =u} Π t n (u, v) - H h=1 1 {X h t-1 =u;X h t =v} .
We choose a quadratic distance and then we minimize the following function:

n t=1 u∈A k v∈A H h=1 1 {X h t-1 =u} Π t n (u, v) -1 {X h t =v} 2 .
Then, we obtain (18) and ( 19) (similarly as ( 15) and ( 16) in Proposition 8).

2

Consequently, from Equation (2), for any u, v ∈ E and t = 1, . . . , n, the estimator of Π t n (u, v) is given by

Π t n ;H (u, v) = 1 - t n Π 0;(n,H) (u, v) + t n Π 1;(n,H) (u, v), (20)
with Π 0;(n,H) (u, v) and Π 1;(n,H) (u, v) given in ( 18) and ( 19), respectively.

Theorem 1 For a linear drifting Markov chain (X k ) 0≤k≤n between two Markov transition matrices (of order 1) Π 0 and Π 1 , for any states u, v ∈ E and time t, t = 1, . . . , n, the estimators of Π 0 (u, v), Π 1 (u, v) and Π t n (u, v) given in Equations ( 18), [START_REF] Limnios | Empirical estimators of reliability and related functions for semi-Markov systems[END_REF] and [START_REF] Limnios | Nonparametric estimation of some important indicators in reliability for semi-Markov processes[END_REF], respectively, are strongly consistent, i.e.,

Π 0;(n,H) (u, v) a.s. ----→ H→∞ Π 0 (u, v), Π 1;(n,H) (u, v) a.s. ----→ H→∞ Π 1 (u, v), Π t n ;H (u, v) a.s. ----→ H→∞ Π t n (u, v).
Proof.

First, for any states u, v ∈ E and time t = 1, . . . , n, note that from the Strong law of large numbers and from the definition of a drifting Markov chain, we immediately have

1 H H h=1 1 {X h t-1 =u} a.s. ----→ H→∞ P(X t-1 = u), 1 H H h=1 1 {X h t-1 =u,X h t =v} a.s. ----→ H→∞ P(X t-1 = u)Π t n (u, v). ( 21
)
Second, we obtain the following limits

Π 0;(n,H) (u, v) a.s. ----→ H→∞ (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u)Π t 2 n (u, v) t 1 n t 1 n 1 - t 2 n - t 2 n 1 - t 1 n (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u) 1 - t 1 n t 2 n t 2 n - t 1 n = (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u)Π 0 (u, v) 1 - t 2 n t 1 n t 1 n - t 2 n (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u) 1 - t 1 n t 2 n t 2 n - t 1 n + (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u)Π 1 (u, v) t 2 n t 1 n t 1 n - t 2 n (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u) 1 - t 1 n t 2 n t 2 n - t 1 n Thus Π 0;(n,H) (u, v) a.s.
----→ H→∞ Π 0 .

Π 1;(n,H) (u, v) a.s. ----→ H→∞ (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u)Π t 2 n (u, v) 1 - t 1 n t 2 n - t 1 n (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u) 1 - t 1 n t 2 n t 2 n - t 1 n = (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u)Π 0 (u, v) 1 - t 2 n 1 - t 1 n t 2 n - t 1 n (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u) 1 - t 1 n t 2 n t 2 n - t 1 n = (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u)Π 1 (u, v) t 2 n 1 - t 1 n t 2 n - t 1 n (t1,t2)∈ 1,n 2 ,t1 =t2 P(X t1-1 = u)P(X t2-1 = u) 1 - t 1 n t 2 n t 2 n - t 1 n . Thus Π 1;(n,H) (u, v) a.s.
----→ H→∞ Π 1 .

2

In the case of H i.i.d. sample paths of a linear drifting Markov chain (X k ) 0≤k≤n , it is immediate to show that the MLE of the initial distribution of the chain is given by

α H (i) = N (H) i,0 H , (22)
where

N (H) i,0 := H l=1
         n t=k 1uA 0 (t)A 0 (t) • • • n t=k 1uA 0 (t)A d (t) . . . . . . n t=k 1uA d (t)A 0 (t) • • • n t=k 1uA d (t)A d (t)               Π 0;n (u 1 , . . . , u k , v)
. . .

Π 1;n (u 1 , . . . , u k , v)      =          n t=k A 0 (t)1uv . . . n t=k A d (t)1uv          . ( 23
)
Consequently, from Equation ( 4), for any u 1 , . . . , u k , v ∈ E and t = 1, . . . , n, the estimator of Π t n (u1, . . . , u k , v) is given by

Π t n (u 1 , . . . , u k , v) = d i=0 A i (t) Π i d ;n (u 1 , . . . , u k , v), u 1 , . . . , u k , v ∈ E, (24)
with Π i d ;n (u, v) given by solving the system [START_REF] Lisnianski | Multi-state system reliability analysis and optimization for engineers and industrial managers[END_REF]. Second, as in the framework of Proposition 9, let us now consider H i.i.d. sample paths of a polynomial drifting Markov chain of order k and degree d, (X k ) 0≤k≤n . That is, we assume that the sample paths have the same length n and that the polynomial drift within the H sample paths takes place between the same d Markov transition matrices (of

order k) Π i d
. The proof is similar to proof of Proposition 9.

Proposition 11 Let us note 1 h u = 1 {X h t-k ...X h t-1 =u1...u k } and 1 h uv = 1 {X h t-k ...X h t-1 =u1...u k ,X h t =v} .
Under this setting, starting with H sample paths of a polynomial drifting Markov chain of order k and degree d, for any states u 1 , . . . , u k , v ∈ E, the estimators of Π i d (u 1 , . . . , u k , v) are obtained by solving the linear system

          n t=k H h=1 1 h A 0 (t)A 0 (t) • • • n t=k H h=1 1 h u A 0 (t)A d (t) . . . . . . n t=k H h=1 1 h u A d (t)A 0 (t) • • • n t=k H h=1 1 h u A d (t)A d (t)                Π 0;(n,H) (u 1 , . . . , u k , v) . . . Π 1;(n,H) (u 1 , . . . , u k , v)      =           n t=k A 0 (t) H h=1 1 h uv . . . n t=k A d (t) H h=1 1 h uv           . (25)
Proposition 12 For a polynomial drift Markov chain of order k, degree d and length n within the H sample paths, the estimators

Π i d ;(n,H) (u, v) verify v∈E Π i d ;(n,H) (u, v) = 1.
Proof.

From the system of Proposition 11, we provide the following system by summing over v :

          n t=k H h=1 1 h u A 0 (t)A 0 (t) • • • n t=k H h=1 1 h u A 0 (t)A d (t) . . . . . . n t=k H h=1 1 h u A d (t)A 0 (t) • • • n t=k H h=1 1 h u A d (t)A d (t)                   v∈A Π 0;(n,H) (u, v) . . . v∈A Π 1;(n,H) (u, v)         =           n t=k H h=1 A 0 (t)1 h u . . . n t=k H h=1 A d (t)1 h u           .
The solution of this system is (1, . . . , 1) .

Estimation for other types of data

As we said at the beginning of this section, different types of data can be considered, that we have denoted by A1, A2, B1, B2 and B3. Note that the estimators of Π 0 and Π 1 in the framework A1 were given in Proposition 8, while in the framework B1 were given in Proposition 9 and their convergence was the object of Theorem 1. In the sequel we will provide the estimators of Π 0 and Π 1 in the other frameworks, namely A2, B2, B3 and B4. The proofs of these results are similar to those of Propositions 8 and 9. For the sake of simplicity, we confine here to a linear drifting Markov model of order 1. The polynomial and k order Markov case can be considered in a similar way.

Proposition 13 Under the statistical framework A2, for (X 0 , X 1 , . . . , Xm), m ≤ n, an incomplete sample path of a linear drifting Markov chain (of order 1) of length n between Π 0 and Π 1 , for any states u, v ∈ E, the estimators of Π 0 (u, v) and Π 1 (u, v) are given by:

Π 0;n (u, v) = P 1 (1, m, n)P 2 (1, m, n) -P 3 (1, m, n)P 4 (1, m, n) P 5 (1, m, n)P 1 (1, m, n) -P 3 (1, m, n) 2 , (26)
Π 1;n (u, v) = P 5 (1, m, n)P 4 (1, m, n) -P 3 (1, m, n)P 2 (1, m, n) P 5 (1, m, n)P 1 (1, m, n) -P 3 (1, m, n) 2 , (27)
with m ≤ n.

Proposition 14 Under the statistical framework B2, for (X h

0 , X h 1 , . . . , X h m h ), m h ≤ n, h = 1, . . . , H, H i.i.d. incomplete
sample paths of a linear drifting Markov chain (of order 1) of length n between Π 0 and Π 1 , for any states u, v ∈ E, the estimators of Π 0 (u, v) and Π 1 (u, v) are given by:

Π 0;(n,H) (u, v) = P 1 (H, m, n)P 2 (H, m, n) -P 3 (H, m, n)P 4 (H, m, n) P 5 (H, m, n)P 1 (H, m, n) -P 3 (H, m, n) 2 , (28)
Π 1;(n,H) (u, v) = P 5 (H, m, n)P 4 (H, m, n) -P 3 (H, m, n)P 2 (H, m, n) P 5 (H, m, n)P 1 (H, m, n) -P 3 (H, m, n) 2 , ( 29
)
with n = (n 1 , . . . , n H ) = (n, . . . , n) and m = (m 1 , . . . , m H ) where m h ≤ n for h = 1, . . . , H.

Proposition 15 Under the statistical framework B3, for (X h 0 , X h 1 , . . . , X h n h ), h = 1, . . . , H, H i.i.d. complete sample paths of H linear drifting Markov chains (of order 1), each of them of length n i , i = 1, . . . , H, between the same Π 0 and Π 1 , for any states u, v ∈ E, the estimators of Π 0 (u, v) and Π 1 (u, v) are given by:

Π 0;(n,H) (u, v) = P 1 (H, n, n)P 2 (H, n, n) -P 3 (H, n, n)P 4 (H, n, n) P 5 (H, n, n)P 1 (H, n, n) -P 3 (H, n, n) 2 , (30)
Π 1;(n,H) (u, v) = P 5 (H, n, n)P 4 (H, n, n) -P 3 (H, n, n)P 2 (H, n, n) P 5 (H, n, n)P 1 (H, n, n) -P 3 (H, n, n) 2 , (31)
with n = (n 1 , . . . , n H ).

Proposition 16 Under the statistical framework B4, for (X h 0 , X h 1 , . . . , X h m h ), m h ≤ n h , h = 1, . . . , H, H i.i.d. incomplete sample paths of H linear drifting Markov chains (of order 1), each of them of length n i , i = 1, . . . , H, between the same Π 0 and Π 1 , for any states u, v ∈ E, the estimators of Π 0 (u, v) and Π 1 (u, v) are given by:

Π 0;(n,H) (u, v) = P 1 (H, m, n)P 2 (H, m, n) -P 3 (H, m, n)P 4 (H, m, n) P 5 (H, m, n)P 1 (H, m, n) -P 3 (H, m, n) 2 , (32)
Π 1;(n,H) (u, v) = P 5 (H, m, n)P 4 (H, m, n) -P 3 (H, m, n)P 2 (H, m, n) P 5 (H, m, n)P 1 (H, m, n) -P 3 (H, m, n) 2 , ( 33
)
with m = (m 1 , . . . , m H ) and n = (n 1 , . . . , n H ).

As we did before in Equations ( 17) and [START_REF] Limnios | Nonparametric estimation of some important indicators in reliability for semi-Markov processes[END_REF]

Estimation of reliability and survival measures

Using the expression of the reliability indicators of a drifting Markov chain obtained in Subsection 2.2 and the estimators obtained in Subsection 3.2 for the characteristics of a drifting Markov chain, one immediately obtain the associated plug-in estimators. We will write in the sequel these estimators for a linear drifting Markov chain of order 1, under the statistical framework B1. We will also give some comments in Remark 6 about the associated estimation under the other frameworks.

indicators, one has to consider (fix) a certain value of the length n of the drifting model of interest (the one that appears in formulas like Formula (5) for the reliability).

Numerical example

In this section we carry out numerical experiments in order to illustrate the quality of the proposed estimators and to numerically validate our procedures. Thus, in Section 4.1 we investigate the quality of the estimators of Π 0 and Π 1 .

Then we present the estimators of reliability indicators in Section 4.2. We conclude by retrieving results on reliability of discrete-time Markov systems obtained in [START_REF] Sadek | Asymptotic properties for maximum likelihood estimators for reliability and failure rates of Markov chains[END_REF] as a particular case of our model.

Let us consider a system with state space E = {1, 2, 3, 4}, partitioned into the up-states U = {1, 2, 3} and down-states D = {4}. In order to undertake numerical illustration of the results presented in this paper, we consider four different linear drifting models, M (1) , M (2) , M (3) and M (4) , on the state space E, each of them varying between Π

i , i = 0, 1, k = 1, . . . , 4 :

d(Π (k) i , Π (k)
i;(n,H) ) = u,v∈E Π (k) i (u, v) -Π (k) i;(n,H) (u, v) 2 . (39)
The obtained distances are given in Table 1. We observe a convergence with the number of sequences H for all the models, fact that is in accordance with Theorem 1. Note that even for small values of H we have excellent results for the convergence.

Model H Π 0 Π 1 5

Reliability measures

In order to observe the convergence of the estimators of reliability measures as the number of sequences H increases (Theorem 2), we simulate H sequences (H = 5, 10 and 100) of length n = 1000, according to the model M (1) .

In Figure 1 we have represented the reliability function as well as the corresponding estimators obtained for these values of H. Notice that the estimator approaches the true value as the number of sequences increases, which is in accordance with Theorem 2.

Similarly, in Figure 2 we have represented the availability function as well as the corresponding estimators, while in Figure 3 we have represented the classical failure rate function and associated estimators. Notice that we represented the entire length of the trajectory and also a zoom up to time 30.

We can notice that the estimator approaches the true value as the number of sequences increases, in accordance with Theorem 2. Note also that in both entire figures we observe the drift between the two Markov models.

Specification to Markov case

We would like also to considered the particular case where Π 0 and Π 1 are the same, equal to Π given before. In other words, this is the classical discrete time Markov system for which reliability measures and associated estimation have been investigated in [START_REF] Sadek | Asymptotic properties for maximum likelihood estimators for reliability and failure rates of Markov chains[END_REF].

In Figure 4 we have represented Reliability, Availability and Failure Rate together with the corresponding estimators obtained for a number of trajectories H = 5, 10 and 100. Here, we represented only up to time 50.

We can notice that the obtain results are in agreement with those presented in [START_REF] Sadek | Asymptotic properties for maximum likelihood estimators for reliability and failure rates of Markov chains[END_REF]. The main novelty of the present article comes from proposing in reliability theory a particular case of non-homogeneous Markov model that is obtain through a gradual variation between two or more Markov models. The type of variation we considered ranges from the simplest linear one to a general polynomial.

Several directions for further research can be investigated.

A natural ongoing research direction is developing analogous notions of semi-Markov drifting models and carrying out applications of this kind of models, for instance in reliability theory.

Also another natural research direction is to consider continuous-time Markov drifting models. It is worth mentioning here the importance of having both developments of discrete-time and continuous-time drifting Markov models. On the one hand, note that continuous-time models are those that are used in most of the applications. On the other hand, note that there is also an intrinsic interest in considering discrete-time models (along with continuous-time versions of these models). Indeed, there are applications (in reliability and survival among others) where the time is discrete by nature; we think about the cycles of an industrial application, the various applications with counting processes, the DNA applications, etc. As we have different theories of discrete-time and continuous-time Markov processes, in our opinion it is of interest to have different developments of discrete-time and continuous-time drifting Markov models. Moreover, the techniques used in discrete and continuous time are different.

Another direction is considering an infinite countable state space for the drifting Markov chains ; generally, in reliability/survival analysis studies (in MSS studies in general) the state space is assumed to be finite. Nonetheless, there could be situations where one could be interested in considering systems with infinite number of states. Obviously, although in practice we always deal with finite quantities (states, number of performance rates, etc.), it could happen that this number of states is not necessarily known in advance. Thus a general and unifying framework could be to consider an infinite countable state space for drifting Markov chains and to develop associated probabilistic/statistical theory.

Nonetheless, several problems may arise when working with infinite matrices: the matrix multiplication is not always defined (for obvious reasons of series convergence), it is not always associative, and the existence of left and right inverse is not an obvious matter in the general case. The reader interested in infinite matrices can see the references and discussions in [START_REF] Barbu | Some algebraic methods in semi-Markov processes[END_REF], that deals with countable state space semi-Markov chains.

Finally, a more advanced research direction could be to consider that the drifting system has continuous-state space.

For continuous-state systems, one can refer to [START_REF] Lisnianski | Estimation of boundary points for continuum-state system reliability measures[END_REF] and the references within. This case could be important for applications where the scale of the measures is continuous. As in practice continuous-state systems are tractable with difficulty, people generally do discretizations, by considering associated discrete-state systems instead of continuous-state systems. One can see [START_REF] Lisnianski | A multi-state Markov model for a short-term reliability analysis of a power generating unit[END_REF] for such an application for the reliability of a coal power generating unit.

2

 2

4 . 4 . 1

 441 We consider the following numerical examples for the couples of Markov transition matrices: Convergence of the estimators In order to illustrate the convergence of the estimators of Π 0 and Π 1 as the number of sequences increases, in according to the results obtained in Theorem 1, we simulate H sequences (H = 5, 10 and 100) of length n = 100 for each of the four different models considered, M (i) , i = 1, . . . , 4. We compute the following distance d between each matrix of parameters Π (k) i and their estimators Π

Fig. 2 Fig. 3

 23 Fig. 1 Reliability estimation

  u1,...,u k ,v∈E be d Markov transition matrices (of order k) over a state space E.Definition 2 (polynomial drifting Markov chain of order k and of length n) A sequence X 0 , X 1 , . . . , Xn with state space E = {1, 2, . . . , s} is said to be a polynomial drifting Markov chain of order k and of length n if the distribution of X t , t = 1, . . . , n, is defined by

  with A i polynomials of degree d such as, for any i, j ∈ {0, 1, . . . , d}, A i Polynomials A i are Lagrange polynomials : they are chosen to have Π t

			nj d	= 1 {i=j} .
	Remark 1 Note that the coherence between notations implies the choice of the notation Π i	. Indeed, for t = ni/d, we
				d
	have Π t n	d = Π i	.
	Remark 2

n a stochastic matrix for each 0 ≤ t ≤ n.

1

  {X l 0 =i} represents the number of sample paths starting in state i among the H sample paths. 10 Let us note 1u = 1 {X t-k ...Xt-1=u1...u k } and 1uv = 1 {X t-k ...Xt-1=u1...u k ,Xt=v} . For (X 0 , X 1 , . . . , Xn) a sample path of a polynomial drifting Markov chain of order k and degree d, for any states u 1 , . . . , u k , v ∈ E,

	Proposition the estimators of Π i d	(u1, . . . u k , v) are given by solving the following linear system

  for the frameworks A1 and B1, it is clear that in any of these new statistical frameworks, from Equation (2), for any u, v ∈ E one can immediately obtain an estimator of Π t = 1, . . . , n h , h = 1, . . . H, for the frameworks B3 and B4.

		n	(u, v), t = 1, . . . , n, for
	the frameworks A2 and B2, respectively of Π t n h	(u, v), t

Table 1

 1 Convergence of the estimators: different values of the distance[START_REF] Xie | Redefining failure rate function for discrete distributions[END_REF]

			0.0001249	0.0000496
	M (1)	10	0.0001058	0.00003858
		100	0.000004377	0.000002932
		5	0.009501	0.003593
	M (2)	10	0.00304	0.002224
		100	0.0001349	0.0002502
		5	0.005092	0.002636
	M (3)	10	0.001491	0.0008556
		100	0.0000884	0.0001761
		5	0.001448	0.000993
	M (4)	10	0.0006423	0.0004372
		100	0.00006185	0.00009181

Acknowledgments This research work was partially supported by the projects MMATFAS -Mathematical Methods for System Reliability and Survival Analysis (2014-2015) and MOUSTIC-Random Models and Statistical, Informatics and Combinatorics Tools (2016-2019), within the Large Scale Research Networks from the Region of Normandy, France. The authors would like to thank their colleagues Dr. Anatoly Lisnianski and Dr. Gregory Levitin from the Israel Electric Corporation Ltd. for pointing out important aspects and references on multi-state systems and also to Prof. Ilia Frenkel and his colleagues from Sami Shamoon College of Engineering, Israel, for organizing the series of conferences International Symposium on Stochastic Models in Reliability Engineering, Life Science and Operations Management, to which, indirectly, the research developed in this article owes much.

Availability Estimator for H=5 Estimator for H=10 Estimator for H=100 True value Availability Estimator for H=5 Estimator for H=10 Estimator for H=100 True value Availability Estimator for H=5 Estimator for H=10 Estimator for H=100 True value

Proposition 17 Under the statistical framework B1, for (X h 0 , X h 1 , . . . , X h n ), h = 1, . . . , H, H i.i.d. complete sample paths of a linear drifting Markov chain (of order 1) of length n between Π 0 and Π 1 , the estimators of the reliability, availability, maintainability, BMP-failure rate, RG-failure rate at time l, l ∈ N, are given by:

for R (n,H) (l) = 0 (if R (n,H) (l) = 0, we set r (n,H) (l) := 0), where, for any states u, v ∈ E, the estimators α

were obtained in Equations ( 22), ( 18) and ( 19),

) .

Using the strong consistency of the estimators ----→ H→∞ R(l); A (n,H) (l) a.s.

----→ H→∞ A(l); M (n,H) (l) a.s.

----→ H→∞ M (l); λ (n,H) (l) a.s.

----→ H→∞ λ(l); r (n,H) (l) a.s.

----→ H→∞ r(l).

Remark 6

In Proposition 17 we considered the estimators of the reliability indicators under the statistical framework B1, that is starting with several complete sample paths of a linear drifting Markov chain (of order 1) of the same length between Π 0 and Π 1 . We would like to give some comments about the other statistical frameworks.

1. Note that the initial distribution α can be estimated only if we have at our disposal several sample paths of the drifting Markov chain (the statistical frameworks B1-B4). In the other cases, one can consider that the initial distribution is known (given), as it is the case when we know that the system starts in a perfect working state. Another possibility is to consider that the initial distribution is the stationary distribution of the Markov transition matrix Π 0 .

2. On the one hand, if we are under the statistical frameworks A1, A2 or B2, we can obtain the estimators of the reliability indicators as we did in Proposition 17 for the statistical framework B1, by considering the estimators of Π 0 and Π 1 given in Propositions 8, 13 and 14, respectively.

3. On the other hand, if we are under the statistical frameworks B3 and B4, obviously that we can estimate the transition matrices Π 0 and Π 1 using Propositions 15 and 16, respectively. Note nonetheless that these two frameworks assume that we start from data issued from different linear drifting Markov chains (consequently, of different theoretical lengths) between the same Π 0 and Π 1 . Thus, in order to obtain the corresponding plug-in estimators of the reliability


Jae-Ho Jung
email: jhjung@gwnu.ac.kr

Ji Hye Moon

Kyung-Min Park

Sanghee Kim

John R Dolan

Eun Jin Yang

Novel insights into the genetic diversity of Parafavella based on mitochondrial CO1 sequences

Keywords: cryptic species, marine plankton, polymorphism, Tintinnidea, 2, JUNG ET AL

  L'archive ouverte pluridisciplinaire

INTRODUCTION

There is a great deal of interest in assessing boreal flora and fauna given the rapid changes that have already occurred and are anticipated to continue in high latitude ecosystems. In a recent report on arctic marine biodiversity, the chapter on plankton [START_REF] Lovejoy | State of the Arctic marine biodiversity: Key findings and advice for monitoring. conservation of Arctic flora and fauna[END_REF] specifically points out the scarcity and need for species data on phytoplankton and heterotrophic protists of the microzooplankton. Here we report on some tintinnid species of the microzooplankton commonly found in waters of the Arctic Ocean as well as the North Pacific and North Atlantic Oceans, forms that have long considered being 'troublesome'. These are forms of the genus Parafavella. The trouble investigators have had is distinguishing morphotypes described as distinct species when intermediate forms appear to exist. This trouble has existed for some time now. Likely the first skeptic was [START_REF] Ostenfeld | Marine plankton from the East-Greenland Sea (W. of 6° W. long. and N. of 73°30' N. lat.) collected during the Danmark Expedition 1906-1908: II Protozoa[END_REF] reporting on the existence of morphologically intermediate forms between those described as 'species' of Parafavella earlier by [START_REF] Brandt | Die Tintinnen[END_REF]. Later Schultz andWulf (1927, 1929) in two large studies argued that there was but one species of Parafavella based on the variety of loricas found in the Barents Sea. [START_REF] Burkovsky | Variability of Parafavella denticulata in the White Sea[END_REF] postulated that lorica morphology of Parafavella changed seasonally in the White Sea. [START_REF] Gold | Tintinnida of the New York Bight: Loricae of Parafavella gigantea, P. parundentata, and Ptychocylis obtusa[END_REF] documented a large variability in empty Parafavella lorica. The last concerted attempts to cast doubt on lorica-morphology based species boundaries in these common boreal tintinnids were studies by Davis concerning Parafavella [START_REF] Davis | Variations of the lorica in the genus Parafavella (Protozoa: Tintinnida) in northern Norway waters[END_REF].

His remarks echoed those of earlier researchers; Davis stated "... I experienced many vexations in attempting the classification of specimens." [START_REF] Davis | Variations of the lorica in the genus Parafavella (Protozoa: Tintinnida) in northern Norway waters[END_REF](Davis, , p. 1823)).

Adding to the troubles and confusion of those attempting to put names on forms is the fact that in the classic monograph on tintinnid ciliates, [START_REF] Kofoid | A conspectus of the marine and freshwater Ciliata belonging to the suborder Tintinnoinea, with descriptions of new species principally from the AGASSIZ expedition to the eastern tropical Pacific, 1904-1905[END_REF] many of the illustrations are not drawn to the scale of the source illustrations (e.g., Figure S1 showing some Parafavella spp.). Thus investigators may assign different names to the same morphotype depending on apparent relative sizes (with no indication of variability) as shown in [START_REF] Kofoid | A conspectus of the marine and freshwater Ciliata belonging to the suborder Tintinnoinea, with descriptions of new species principally from the AGASSIZ expedition to the eastern tropical Pacific, 1904-1905[END_REF] or the dimensions (and notes on variability) given in the original descriptions.

In recent years nuclear ribosomal genes have been applied to distinguish ciliate species and infer their phylogeny as DNA barcode [START_REF] Pawlowski | CBOL protist working group: Barcoding eukaryotic richness beyond the animal, plant, and fungal kingdoms[END_REF][START_REF] Stoeck | The D1-D2 region of the large subunit ribosomal DNA as barcode for ciliates[END_REF]. Ciliates are suspected to contain a huge species diversity considering estimates of the numbers of as yet undescribed species [START_REF] Foissner | Diversity and geographic distribution of ciliates (Protista: Ciliophora)[END_REF].

3, JUNG ET AL.

An existing challenge is to separate closely related species efficiently using molecular methods. It has become apparent that 18S rDNA, one of most common genes for ciliate taxonomy/phylogeny, appears to have an insufficient resolution for species identification. For instance, the 18S rDNA is completely identical between some species of Tetrahymena [START_REF] Lynn | Species of Tetrahymena identical by small subunit rRNA gene sequences are discriminated by mitochondrial cytochrome c oxidase I gene sequences[END_REF] and shows more than 99.9% genetic similarity between some congeners in Pleurotricha [START_REF] Park | Morphology, 18, JUNG ET AL. morphogenesis, and molecular phylogeny of Pleurotricha oligocirrata nov. spec. (Ciliophora: Spirotrichea: Stylonychinae)[END_REF].

Because of the problematic resolution, there have been efforts to establish new genetic marker as a barcode for ciliates [START_REF] Barth | Intraspecific genetic variation in Paramecium revealed by mitochondrial cytochrome c oxidase I sequences[END_REF][START_REF] Stoeck | The D1-D2 region of the large subunit ribosomal DNA as barcode for ciliates[END_REF]. For tintinnids, Santoferrara et al. (2013) assessed the nuclear genes as a potential barcoding tool, and she and colleagues (2016) suggested a guideline for data collection and evaluation. In addition, [START_REF] Santoferrara | Discrimination of closely related species in tintinnid ciliates: New insights on crypticity and polymorphism in the genus Helicostomella[END_REF] reported the crypticity and polymorphism in the genus Helicostomella based on the ribosomal genes.

As mentioned by [START_REF] Santoferrara | Updating biodiversity studies in loricate protists: The case of the tintinnids (Alveolata, Ciliophora, Spirotrichea)[END_REF], the mitochondrial cytochrome c oxidase subunit I (CO1) gene had not yet been amplified for the tintinnids. After the first report of intraspecific genetic variation using the CO1 in Paramecium [START_REF] Barth | Intraspecific genetic variation in Paramecium revealed by mitochondrial cytochrome c oxidase I sequences[END_REF], Strüder-Kypke and Lynn (2010) provided a primer pair, which successfully amplifies the CO1 from five ciliate classes. Considering its greater intraspecific variation than any nuclear genes, the CO1 can shed light on population genetics and biogeography. However, it should be noted that CO1 does not satisfy the CBOL (Consortium for the Barcode of Life)'s criteria as a barcode for all ciliates. [START_REF] Stoeck | The D1-D2 region of the large subunit ribosomal DNA as barcode for ciliates[END_REF] noted that some ciliates from anoxic habitats lack active mitochondria and CO1. Furthermore, not all ciliate CO1 can be amplified using the existing CO1 primers (e.g. Litostomatea, Spirotrichea). Rather than as a ciliate barcode, the CO1 can be utilized for the seperation of closely related species or populations targeting specific ciliates. [START_REF] Park | Utility of mitochondrial CO1 sequences for species discrimination of Spirotrichea ciliates (Protozoa, Ciliophora)[END_REF] provided a new primer pair, which can amplify the class Spirotrichea, and they successfully amplified Tintinnopsis cf. cylindrica. Here these new primers were applied to resolve the problematic situation among Parafavella forms.

To clarify the relationships among morphotypes of Parafavella, we employed single cell sequencing allowing us to investigate the genetic relationships of some more or less distinct morphotypes of Parafavella. We sequenced both nuclear ribosomal genes (28S rDNA, 18S rDNA, 5.8S rDNA, internal transcribed spacer regions: ITS1, ITS2) as well as mitochondrial CO1 sequences to group or distinguish 4, JUNG ET AL. morphotypes of Parafavella. We also investigated the co-occurrence of the more or less distinct morphotypes in natural populations by examining material from sites sampled in a transect from the North Pacific into the High Arctic in the summer of 2017.

MATERIALS AND METHODS

Sampling

For the genetic analysis, plankton net-tow samples were obtained from three distinct seas using 20-µm mesh plankton net (for details, see Figure 1; Table 1). The samples were fixed with neutral Lugol's solution (final concentration of 0.04% iodine [w/v]; 0.6% potassium iodide [w/v]) or ethanol (final concentration of 70%).

For the analysis of Parafavella forms found in natural assemblages of tintinnids, plankton net tow samples (20-mesh plankton net) obtained from stations along a transect from the North Pacific (44°N) to the High Arctic (78°N) sampled in July 2017 from aboard the Korean Icebreaker Araon were analyzed. Net tows sampled the water column from 100 m depth to the surface except at some shallow water stations. Aliquots of net tow material for microscopic analysis were fixed with Lugol's Iodine (as above) solution.

Morphology

For genetic analysis, individual cells were first isolated using a micro-capillary under a stereomicroscope (SZ11, Olympus, Japan), and micrographs were captured and measured under an inverted microscope (Eclipse Ti-U, Nikon, Japan) or standard microscopes (Zeiss Axio Imager 2 Carl Zeiss, Germany; BX53, Olympus, Japan) to record the morphology of the cells sequenced.

For analysis of the occurrence of forms of Parafavella in natural assemblages,, aliquots of net tow material from the 2017 sampling were diluted and examined in 3 ml settling chambers. Material in the settling chambers was examined using an Olympus inverted microscope (IX71) equipped with differential interference contrast optics and an Olympus DP71 camera piloted using Olympus Biocell imaging software.

Initially all tintinnids in a single aliquot were enumerated. For samples in which 5, JUNG ET AL.

Parafavella were found (samples from 9 stations), multiple aliquots were examined.

For these latter samples, all Parafavella encountered which were correctly oriented for morphological diagnoses were photographed and multiple aliquots examined until a total of 20-30 cells from the sample were photographed. Only lorica containing a ciliate cell were recorded as empty lorica are not a reliable evidence of living cells (Dolan & Yang, 2017). Parafavella "species" identification (based upon lorica morphology) were made relying as much as possible on illustrations and dimensions as given in the original descriptions (for details see "Parafavella "species" designations -taxonomic notes" in the Supplemental Files).

DNA extraction, PCR, and DNA sequencing

Each individual was transferred to 1.5 ml tube and its genomic DNA was extracted using RED-Extract-N-Amp Tissue PCR Kit (Sigma, St. Louis, MO, USA) according to the manufacturer's instruction. Two PCRs were conducted to amplify nuclear ribosomal genes and mitochondrial CO1 gene. The conditions for the ribosomal genes were as follows: denaturation at 94 for 1 min 30 s, followed by 40 cycles of denaturation at 98 for 10 s, annealing at 58.5 for 30 s, and extension at 72 for 3 min, and a final extension step at 72 for 7 min. The PCR product that covered the genes consisting of 18S rDNA, internal transcribed spacer regions (ITS1, ITS2), 5.8S rDNA, and 28S rDNA (D1, D2) was amplified using slightly modified two primers (New Euk A and LSU rev2) that were described by [START_REF] Sonnenberg | An evaluation of LSU rDNA D1-D2 sequences for their use in species identification[END_REF]. DNA sequencing was performed using three additional internal primers (18SF790v2: 5´-AAA TTA KAG TGT TYM ARG CAG-3´ and 18SR300: 5´-CAT GGT AGT CCA ATA CAC TAC-3´, 18SF1470: 5´-TCT GTG ATG CCC TTA GAT GTC-3´) and an ABI 3700 sequencer (Applied Biosystems, Foster City, CA, USA). For the conditions of CO1, it was as follows: denaturation at 94 for 1 min 30 s, followed by 40 cycles of denaturation at 98 for 10 s, annealing at 53 for 30 s, and extension at 72 for 1 min, and a final extension step at 72 for 7 min (for details and primers see [START_REF] Park | Utility of mitochondrial CO1 sequences for species discrimination of Spirotrichea ciliates (Protozoa, Ciliophora)[END_REF]. All PCR products were purified using a MEGAquick-spin Total Fragment DNA Purification Kit (iNtRON, Korea). Sequence fragments were assembled using Geneious 9.1.6 [START_REF] Kearse | Geneious Basic: an integrated and extendable desktop 17, JUNG ET AL. software platform for the organization and analysis of sequence data[END_REF].

Molecular data analyses

6, JUNG ET AL.

To infer the phylogenetic position of the four tintinnid species, 18S rRNA gene sequences of 199 ciliates were retrieved from the NCBI database, including 197 tintinnids and two oligotrichs, Novistrombidium orientale (FJ422988) and Strombidium stylifer (DQ631805), as outgroups. The sequences were aligned using Muscle alignment [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF] in Geneious, and both ends of the alignment were manually trimmed using the Geneious. The best-fit model of substitution for phylogenetic analysis was selected using jModelTest 2.1.10 ( [START_REF] Darriba | jModelTest 2: more models, new heuristics and parallel computing[END_REF]. We selected the model TIM2 + I (0.5690) + G (0.4550) for 18S rDNA and HKY + I (0.6480) for CO1 based on the Akaike information criterion (AIC). As an outgroup of CO1 tree, Tintinnopsis cf. cylindrica (MG594903) was selected. IQ-TREE 1.5.3 was used to render maximum likelihood (ML) trees, with 1000 bootstrap replicates [START_REF] Nguyen | IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF]. A consensus ML tree was annotated using ggtree 1.8.2 in R (Yu et al., 2017). In addition, the tree using concatenated 28S rDNA, 5.8S rDNA and ITS regions was inferred with the best-fit model TrN + G (0. 1040). Laackmanniella prolongata JQ924056 was selected as an outgroup. Uncorrected pairwise distances were measured using the Geneious.

To infer the median-joining network of CO1 haplotype, POPART 1.7 was used [START_REF] Leigh | POPART: full-feature software for haplotype network construction[END_REF]. Input file includes the CO1 nucleotide sequences with its sampling location (latitude, longitude) to plot them in a map. The program assigned the CO1 sequences into each haplotype and plotted in default map.

In R, the ggplot2 [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF] was used to plot the 30 specimens annotated with its CO1 haplotype and species name to examine the relationship between the haplotypes and morphology. To represent sequence difference matrices, eight sequence data sets were set up to discern a resolution of species identification for Parafavella and they were as follows: 18S rDNA (1,745 bp), ITS1 (101 bp), 5.8S (152 bp), ITS2 (194 bp), 28S (D1, D2; 753 bp), 28S (D1, D2; 594 bp, shorter than previous one to make the sequences correspond with the reference sequence P. parumdentata), CO1 nucleotide (478 bp), and CO1 amino acids (159 aa) sequences.

The sequence difference matrices measured by the Geneious were used to construct heatmaps using lattice in R [START_REF] Sarkar | Utility of genetic markers and morphology for species discrimination within the order Tintinnida (Ciliophora, Spirotrichea)[END_REF].

RESULTS

Lorica morphology of sequenced specimens

7, JUNG ET AL.

Four morphotypes of Parafavella were collected from three distinct seas and here their lorica morphologies are described briefly (Figure 2; Table 1). See the Supplemental Document 1 for further details on distinguishing the morphotypes encountered.

Parafavella jorgenseni was collected from Bering Sea. Three specimens were sequenced: Lorica conical with reticulate surface, denticulate tilted outward on oral rim, 43.1-48.0 µm in lorica oral diameter (LOD), 96.7-106.3 µm in length, aboral end sharply pointed.

Parafavella gigantea was collected from the Bering Sea, Greenland Sea and East/Japan Sea. Twenty one specimens were sequenced: Lorica cylindrical/conical with reticulate surface, denticulate tilted outward on oral rim, 49.3-62.9 µm in LOD, 269.0-511.1 µm in length, aboral end pointed; in Bering and Greenland Sea specimens lorica usually cylindrical without indentation on oral rim, long cylindrical lorica with relatively shorter aboral end (vs. lorica more conical, aboral end longer with the indentation on oral rim in the East/Japan Sea specimens).

Parafavella greenlandica was collected from Greenland Sea. Two specimens were sequenced: Lorica conical with reticulate surface, denticulate tilted outward on oral rim, 33.2-34.8 µm in LOD, 78.4-86.4 µm in length, aboral end somewhat bluntly pointed.

Parafavella subrotundata was collected from Bering Sea. Four specimens were examined and their morphology are as follows: Lorica cylindrical with reticulate surface and indentation on oral rim, denticulate tilted outward on oral rim, 57.4-63.4 µm in LOD, 163.2-250.8 µm in length, aboral end short and pointed; one of the examined specimens with aboral end being toward inside of lorica.

Two small morphotypes, P. greenlandica and P. jorgenseni, were easily separated from each other by the LOD (Table 1). In contrast, the other larger morphotypes, P. gigantea and P. subrotundata, share a similar size of LOD. In particular, P. gigantea was collected from all seas and they can be split into two morphotypes. The specimens from Bering and Greenland Sea usually have the largely cylindrical lorica without the indentation on oral rim, while the Korean population has the slightly conical lorica with the indentation on the oral rim. However, these 8, JUNG ET AL.

characteristics are not distinct across whole specimens. See the specimens P. gigantea Arc07 and Arc09 had slightly conical lorica.

Nuclear ribosomal genes

A total of six data sets were analyzed to examine genetic relationships of the five morphotypes including P. parumdentata retrieved from GenBank (Figure 3). 18S and 5.8S gene sequences for the five species were completely identical (Figure 4A). It should be noted that one nucleotide difference in 18S of P. parumdentata in Figure 3 denotes an ambiguous one (n).

ITS1 and 2 yielded two clusters based on one (1.0%) and five (2.6%) nucleotide differences, respectively (Figure 3,4B). The two clusters consist of P. greenlandica + P. parumdentata and P. jorgenseni + P. gigantea + P. subrotundata, that is the intra/inter-specific variations were 0% and 0-2.6%, respectively.

28S data sets consist of two alignments because of the sequence length difference of P. parumdentata. The other sequences obtained newly here are longer than the P. parumdentata sequence in GenBank. Thus, the main differences between the data sets are the sequence length and absence/presence of P. parumdentata and so '28S' includes 753 bp of 31 sequences and '28S trimmed' includes 594 bp of 30 sequences. Both data sets showed the same 3 clusters these clusters: 1) P. greenlandica + P. parumdentata, 2) P. gigantea (Greenland Sea), 3) P. gigantea (Bering Sea, East/Japan Sea) + P. subrotundata + P. jorgenseni. The intra/interspecific variations were 0-0.1% and 0-1.9% for '28S', and 0-0.2% and 0-2.4% for '28S trimmed'. The cluster 1, 'P. greenlandica and P. parumdentata', share completely identical nucleotide sequences except for the gap/deletion of a nucleotide in P. parumdentata. Except for the gap, this cluster differed at 12 (2.0%) and 13

(2.2%) nucleotide positions from P. gigantea (Greenland Sea) and P. gigantea (Bering Sea, East/Japan Sea) + P. subrotundata + P. jorgenseni, respectively. A single nucleotide of the cluster P. gigantea (Greenland Sea) differed from P. gigantea (Bering Sea, East/Japan Sea) and P. subrotundata and P. jorgenseni (Figure 3,4B).

The data from concatenated 28S, 5.8S, ITS regions showed the same clusters as the 28S data set because the cluster P. greenlandica + P. parumdentata based on ITS regions was completely identical to one of the three clusters based on the 9, JUNG ET AL. concatenated data and the 5.8S sequences have no variations (Figure 4B). The three clusters from the concatenated did not reflect distinct geographical distribution.

Mitochondrial CO1 gene

CO1 nucleotide data showed eight haplotypes while the amino acids set were all completely identical to each other (Figures 3,4C,D). Based on morphological identity, the intra/inter-specific variations of the nucleotide sequences were 0-9.4% and 0-18.2%, respectively. The CO1 nucleotide sequence of Tintinnopsis cf. cylindrica as an outgroup showed inter-specific variations of 20.5-23.4% while the amino acid sequence showed 1.3% from those of Parafavella.

Parafavella gigantea consisted of four haplotypes that clustered in two groups: 1) Par07 + Par08 and 2) Par04 + Par05. These two groups had 43-45 nucleotide differences (9.0%-9.4%) while each group had a single nucleotide difference (0.2%) within them. The two sequences of P. greenlandica were completely identical and assigned to Par06. The three sequences of P. jorgenseni split into three haplotypes. Of these, Par02 showed nine (1.9%) nucleotide differences from Par03 while Par01 was highly distinct from Par02 + Par03 with 41 (8.6%) or 42 (8.8%) differences. Four sequences of P. subrotundata showed completely identical sequences with P. gigantea (East/Japan Sea) that assigned into Par07.

These eight haplotypes clustered into five groups based on the maximum likelihood tree and median-joining network using 97% similarity threshold (Figure 4): Par01 (P. jorgenseni); Par02 + Par03 (P. jorgenseni); Par04 + Par05 (P. gigantea from Greenland Sea); Par06 (P. greenlandica); and Par07 + Par08 (P. gigantea from all three seas and P. subrotundata). As shown above, the barcode gap between intraand inter-specific variations was not observed.

Biogeography of Parafavella based on genes

Of the nuclear ribosomal genes, 28S and ITS regions showed nucleotide variations only within and among Parafavella morphotypes (Figure 4 B). One nucleotide difference was observed between the East/Japan Sea + the Bering Sea and the Greenland Sea for P. gigantea. The concatenated nucleotide sequences of Parafavella were completely identical within each population of the East/Japan and Bering Seas, respectively. Within the Greenland Sea population, two clusters were identified and 10, JUNG ET AL.

each belongs to P. greenlandica and P. gigantea morphotype separated by 18 nucleotide differences.

Of the eight haplotypes in CO1, the Par07 was a common one and occurred from the all seas and groups P. gigantea with P. subrotundata (Figure 4D). The other haplotypes occurred only in one sea. Considering the number of haplotypes in each sea, the Greenland and Bering Sea showed the highest number of haplotypes as four (vs. 2 and 1 in ribosomal genes), respectively. However, the East/Japan Sea was represented by just two haplotypes, which have only one nucleotide difference (vs. identical in ribosomal genes). It should be noted that the nuclear ribosomal genes were too conserved to allow discrimination among populations compared to the CO1 as expected (Figure 4B). However, the eight haplotypes clustering into five groups resulted in neither distinct biogeography nor gradient distribution pattern (Figure 4D).

See Figure 4D, the eight haplotypes (or the five clusters) did not show certain biogeography among haplotypes/clusters.

Co-occurrence of distinct Parafavella morphotypes in natural communities

Parafavella were found in 9 samples from the 2017 cruise of the Araon. None contained only a single morphotype. The number of different Parafavella "species" (based on differences in the loricas) ranged from 2 to 6 within individual samples.

Table 2 data suggests that there was no apparent relationship between the number of Parafavella morphotypes found and latitude. Nor was a relationship between with the overall concentration of tintinnid ciliates or with the concentration of Parafavella cells and the number of morphotypes of Parafavella in a population (for details of the numbers and identities of Parafavella encountered see the Supplemental Plates 1-9).

DISCUSSION

Nuclear ribosomal genes

Previously, of all the species in Parafavella, nuclear gene data was available only for P. parumdentata [START_REF] Santoferrara | Phylogeny, classification and diversity of Choreotrichia and Oligotrichia (Ciliophora, Spirotrichea)[END_REF]. As in their 18S tree, here all five species of Parafavella including four species from this study clustered together because all sequences of Parafavella are completely identical with the exception of a single nucleotide in P. parumdentata (ambiguous nucleotide n).and showed a sister 11, JUNG ET AL. relationship with a cluster " Dadayiella + Xystonella" Considering the five morphotypes in Parafavella, all nuclear ribosomal genes failed to distinguish among them. Even though the 28S showed three clusters, it is still less than the numbers of morphotypes (5) based on lorica morphology and the CO1 nucleotide haplotypes (8).

Of the nuclear genes, the 28S has been considered as a useful barcoding tool for ciliates (Santoferrara et al., 2013;[START_REF] Stoeck | The D1-D2 region of the large subunit ribosomal DNA as barcode for ciliates[END_REF]. Santoferrara et al. (2013) separated tintinnids at 1% cutoff and [START_REF] Stoeck | The D1-D2 region of the large subunit ribosomal DNA as barcode for ciliates[END_REF] identified that <0.6% sequence divergence was ideal to separate Paramecium species. If we follow these thresholds, only two clusters were identified as follows: cluster 1 containing P. greenlandica and P. parumdentata the other cluster containing P. gigantea, P. subrotundata and P. jorgenseni. Among tintinnids, Helicostomella is one of the bestknown genera in terms of genetic data (Santoferrara et al., 2013[START_REF] Santoferrara | Discrimination of closely related species in tintinnid ciliates: New insights on crypticity and polymorphism in the genus Helicostomella[END_REF][START_REF] Xu | Species boundaries in tintinnid ciliates: A case study -morphometric variability, molecular characterization, and[END_REF]. According to [START_REF] Santoferrara | Discrimination of closely related species in tintinnid ciliates: New insights on crypticity and polymorphism in the genus Helicostomella[END_REF], use of these multi-gene barcodes (i.e., 28S, 5.8S-ITS) is better than using a single gene to distinguish closely related species.

However we found that these multi-genes performed little better than a single gene marker in Parafavella yielding only two clusters. Nonetheless we agree with the comment of [START_REF] Santoferrara | Discrimination of closely related species in tintinnid ciliates: New insights on crypticity and polymorphism in the genus Helicostomella[END_REF] that because CO1 is much more variable than these nuclear genes it results in lower resolution at higher level of phylogeny [START_REF] Park | Utility of mitochondrial CO1 sequences for species discrimination of Spirotrichea ciliates (Protozoa, Ciliophora)[END_REF]. In addition, genetic diversity of the CO1 can originate from heteroplasmy (i.e., intracellular polymorphism of mitochondrial genomes), nuclear encoded mitochondrial pseudogenes, and hybrids [START_REF] Tautz | A plea for DNA taxonomy[END_REF]Zhao et al., 2013) The recent development of high-throughput sequencing technology has focused on attention on the variable regions (i.e., V4) in 18S to capture eukaryote diversity [START_REF] Filker | Deep sequencing uncovers protistan plankton diversity in the Portuguese Ria Formosa solar saltern ponds[END_REF][START_REF] Jung | Patchy-distributed ciliate (Protozoa) diversity of eight polar communities as determined by 454 amplicon pyrosequencing[END_REF][START_REF] Santoferrara | Global and local DNA (meta)barcoding reveal new biogeography patterns in tintinnid ciliates[END_REF][START_REF] Stoeck | Multiple marker parallel tag environmental DNA sequencing reveals a highly complex eukaryotic community in marine anoxic water[END_REF]. However, our data supports the view that the nuclear ribosomal genes are of insufficient resolution below the genus level of Parafavella as shown in other tintinnids in previous reports (Santoferrara et al., 2013;[START_REF] Xu | Species boundaries in tintinnid ciliates: A case study -morphometric variability, molecular characterization, and[END_REF]. In addition, there is evidence of low inter-specific variations in hypotrichs, benthic ciliates assigned to the ciliate class containing tintinnids [START_REF] Park | Morphology, 18, JUNG ET AL. morphogenesis, and molecular phylogeny of Pleurotricha oligocirrata nov. spec. (Ciliophora: Spirotrichea: Stylonychinae)[END_REF]Yi et al., 2008).

Considering the 97-98% similarity threshold commonly used by the high-throughput approaches for the species richness, it would appear that only a part of eukaryote diversity is unveiled at least with regard to ciliates.

Mitochondrial CO1 gene

12, JUNG ET AL.

The CO1 of Parafavella was successfully amplified using the CO1 primers targeting spirotrichean ciliates [START_REF] Park | Utility of mitochondrial CO1 sequences for species discrimination of Spirotrichea ciliates (Protozoa, Ciliophora)[END_REF]. Even though they provided only the CO1 of Tintinnopsis cf. cylindrica among tintinnids, our study suggests that other tintinnids could be sequenced using the primers. Strüder-Kypke and Lynn (2010) established new primers to amplify the CO1 of ciliates, however, as they mentioned, they did not reliably amplify Litostomatea and Spirotrichea. The CO1 of Parafavella overlaps 260-348 bp with those region from Strüder-Kypke and [START_REF] Strüder-Kypke | Comparative analysis of the mitochondrial cytochrome c oxidase subunit I (COI) gene in ciliates (Alveolata, Ciliophora) and evaluation of its suitability as a biodiversity marker[END_REF]. Across all sequences of Parafavella, none of indels occurred that resulted in frame shifts or be considered as nuclear encoded mitochondrial pseudogenes (Zhao et al., 2013).

The CO1 sequences of Parafavella showed higher genetic variation than the nuclear ribosomal genes, up to 18.2% in inter-specific variations. Considering the maximum of the inter-specific variations in the 28S, the CO1 showed the variations higher than seven times in Parafavella. A total of eight haplotypes were identified based on the CO1 nucleotide sequences that resulted in five clusters using 97% similarity threshold. If we follow the <1% intra-specific divergence as shown by [START_REF] Lynn | Species of Tetrahymena identical by small subunit rRNA gene sequences are discriminated by mitochondrial cytochrome c oxidase I gene sequences[END_REF] and [START_REF] Chantangsi | Barcoding ciliates: a comprehensive study of 75 isolates of the genus Tetrahymena[END_REF], six clusters can be identified because P. jorgenseni Ara02 and Ara03 have nine nucleotide differences between them (1.9%).

Lorica morphology with CO1 haplotypes

The species designations of Parafavella have been contentious because of suspected lorica variability [START_REF] Davis | Variations of the lorica in the genus Parafavella (Protozoa: Tintinnida) in northern Norway waters[END_REF][START_REF] Kofoid | A conspectus of the marine and freshwater Ciliata belonging to the suborder Tintinnoinea, with descriptions of new species principally from the AGASSIZ expedition to the eastern tropical Pacific, 1904-1905[END_REF][START_REF] Schulz | Hydrographie und Oberflächenplankton des Westlichen Barentsmeeres im Sommer 1927[END_REF]. [START_REF] Kofoid | A conspectus of the marine and freshwater Ciliata belonging to the suborder Tintinnoinea, with descriptions of new species principally from the AGASSIZ expedition to the eastern tropical Pacific, 1904-1905[END_REF] established the genus and assigned species status to 23 forms, most of which had been described originally as variants of other species. The most recent compendium of Zhang et al. ( 2012) lists 32 species found in the literature.

In the WoRMS database [START_REF] Warren | World Ciliophora Database. Parafavella Kofoid & Campbell[END_REF], 35 species are listed including 6 as 'taxon inquirendum', suspected of being variants of other species. [START_REF] Davis | Variations of the lorica in the genus Parafavella (Protozoa: Tintinnida) in northern Norway waters[END_REF] tried to distinguish Parafavella species collected from northern Norway waters, however he failed to distinguish them based on lorica morphology alone. Even though some previous works considered them as a single polymorphic species (P. denticulata), he emphasized the need for further investigation of this genus using methods such as silver staining, electron microscopy, and etc. Here, we found six clusters based the 1% cutoff value from the eight CO1 haplotypes. The results 13, JUNG ET AL.

clearly showed that 1) P. subrotundata has identical nucleotide sequence with P. gigantea; 2) P. gigantea includes two clusters; 3) P. jorgenseni can be split into three clusters; and 4) P. greenlandica is distinct from the others.

Parafavella subrotundata collected from the Bering Sea showed completely identical sequences with P. gigantea in Par07, which is the most common haplotype across the three seas. [START_REF] Kofoid | A conspectus of the marine and freshwater Ciliata belonging to the suborder Tintinnoinea, with descriptions of new species principally from the AGASSIZ expedition to the eastern tropical Pacific, 1904-1905[END_REF] considered that P. subrotundata 'differs from P. cylindrica in stouter proportions and from P. dilatata in less conical bowl'. Based on our CO1 haplotypes, the results suggest that the cylindrical forms in Parafavella might be assigned to P. gigantea.

Of the specimens in our samples, we have P. gigantea with slightly posteriorly dilated lorica like P. ventricosa (see Arc02, Arc05, Arc08, Ara08). Even though the dilated form was not as distinct as shown in [START_REF] Kofoid | A conspectus of the marine and freshwater Ciliata belonging to the suborder Tintinnoinea, with descriptions of new species principally from the AGASSIZ expedition to the eastern tropical Pacific, 1904-1905[END_REF], it suggests they might be conspecific.

Based on the tintinnids records reported in the Arctic Seas by Dolan et al. (2017), Parafavella denticulata is the most commonly reported species in this genus.

Even though we did not obtain the CO1 sequence from P. denticulata, it might be assigned to P. gigantea CO1 clusters like P. subrotundata.

Interestingly, the 2 species with short conical lorica, P. greenlandica and P. jorgenseni, did not cluster together and P. greenlandica was distinct from the other species in Parafavella. This finding suggests that the size and shape of the lorica are poor predictors of genetic relatedness in Parafavella. Considering the utility of oral denticulation in distinguishing species, we observed some P. gigantea with the denticulation, which was almost disconnected from the oral rim during the examination using the scanning electron microscopy (Figure S2). We also found P. gigantea without denticulation in Lugol's-fixed samples (see supplemental Plate 3 ARA08 St4A) in agreement with the observations of [START_REF] Davis | Variations of the lorica in the genus Parafavella (Protozoa: Tintinnida) in northern Norway waters[END_REF]. Thus, oral denticulation does not appear to be reliable character among Parafavella. To date, data on the infraciliature of Parafavella are available only for P. denticulata and P. gigantea [START_REF] Pierce | Morphology and infraciliature of selected species of Tintinnina with a phylogenetic analysis of the Tintinnina based on infraciliature[END_REF]. As mentioned by [START_REF] Santoferrara | Updating biodiversity studies in loricate protists: The case of the tintinnids (Alveolata, Ciliophora, Spirotrichea)[END_REF], infraciliature patterns can shed light on their biodiversity. Unfortunately, commercial production of protargol, the key chemical for protargol staining to examine infraciliature, has ended.

We have tried to stain Parafavella using laboratory-synthesized protargol (Kim & 14, JUNG ET AL. Jung, 2017), but the stain quality was poor for the cells stored for a lengthy time in Bouin's solution. We recommend a focus on the infraciliature of representative haplotypes after screening the CO1 haplotypes rather than examining the infraciliature of all morphotypes.

Co-occurrence of distinct Parafavella morphotypes in natural communities

We found as many as 6 morphotypes, classically distinguished as distinct species, in a single sample. It is important to note that while many previous studies documented the co-occurrence of distinct lorica types of Parafavella (e.g. [START_REF] Ostenfeld | Marine plankton from the East-Greenland Sea (W. of 6° W. long. and N. of 73°30' N. lat.) collected during the Danmark Expedition 1906-1908: II Protozoa[END_REF][START_REF] Schulz | Hydrographische und planktologische Ergebnisse der Fahrt des Fischereischutzbootes "Zeiten" in das Barentsmeer im August-September 1926[END_REF][START_REF] Warren | World Ciliophora Database. Parafavella Kofoid & Campbell[END_REF][START_REF] Burkovsky | Variability of Parafavella denticulata in the White Sea[END_REF][START_REF] Davis | Variations of the lorica in the genus Parafavella (Protozoa: Tintinnida) in northern Norway waters[END_REF] ours is the first study to rely solely on lorica containing ciliate cells. Empty tintinnid lorica are poor indicators of the occurrence and morphology of living tintinnids (Dolan & Yang, 2017;[START_REF] Kato | Tintinnid ciliates as indicator species of different water masses in the western North Pacific Polar Front[END_REF]. However, it unclear what factors may influence the occurrence of distinct morphotypes of Parafavella. There was no obvious relationship between the number of morphotypes found and location in the form of latitude. Notably there was no relationship between the number of morphotypes and either the concentration neither of Parafavella cells nor with concentrations of total tintinnids (Table 1).

Our findings with regard to the occurrence of Parafavella morphotypes contrast to findings with regard to Cymatocylis in Antarctic waters. Like Parafavella of boreal waters, [START_REF] Kim | Morphological and ribosomal DNA-based characterization of six Antarctic ciliate 5 morphopecies from the Amundsen Sea with phylogenetic analyses[END_REF] found different morphotypes of Cymatocylis of the Southern Ocean microzooplankton, previously considered as distinct species, to be variants of a single species. The number of different morphotypes found in populations in the Amundsen Sea was unrelated to the concentration or composition of their phytoplankton prey but was positively related to the concentration of Cymatocylis cells [START_REF] Dolan | Tintinnid ciliates of the Amundsen Sea (Antarctica) Plankton Communities[END_REF]. Culture work with Favella showed that lorica construction by a newly divided cell is rapid (2-10 minutes) and the length and shape of the lorica constructed appears to depend on the amount lorica material within the newly-divided ciliate cell (reviewed in [START_REF] Agatha | Biology and Ecology of Tintinnid Ciliates: Models for Marine Plankton[END_REF]. It was then hypothesized that co-occurrence of multiple morphotypes of Cymatocylis reflected periods of rapid population growth with variable morphotypes formed by rapidly reproducing cells. This does not appear to be the case in Parafavella, as the number of morphotypes found in a population was not related to the concentration of Parafavella. At this time 15, JUNG ET AL.

we can offer no explanation for the apparent polymorphism of certain haplotypes of Parafavella.

Figure legends

Fig. 1 .

 1 Fig. 1. Map of the locations sampled from three different seas. The numbers in parenthesis denote the number of specimens sequenced in this study. The haplotype represents CO1 nucleotide diversity.

Fig. 2 .

 2 Fig. 2. Micrographs of four Parafavella species all at the the same magnification.

Fig. 3 .

 3 Fig. 3. Heatmaps based on genetic dissimilarities.

Fig. 4 .

 4 Fig. 4. Maximum likelihood phylogenetic trees and median-joining network. The network was constructed using CO1 nucleotide sequences.

Fig. 5 .

 5 Fig. 5. Morphological characteristics and CO1 nucleotide haplotypes of four Parafavella species.

ACKNOWLEDGEMENTS

This work was supported by grants from the National Research Foundation of Korea (NRF) funded by the Korea government (MSIP; Ministry of Science, ICT & Future Planning) (No. NRF-2017R1C1B5017183), and the Korea Polar Research Institute (KOPRI)[grant number PE17900]. Partial support was provided by grants from the K-AOOS program (KOPRI; PM18040 & 20160245) funded by the MOF, Korea.

Greenland Sea (n = 12) East/Japan Sea (n = 10) Bering Sea (n = 8) Par01 Haplotype Par02 Par03 Par04 Par05 Par06 Par07 Par08 Ara01 Ara02 Ara03 Ara04 Ara05 Ara06 Ara07 Ara08 Arc01 Arc02 Arc03 Arc04 Arc05 Arc06 Arc07 Arc08 Arc09 Arc10 Kor01 Kor02 Kor03 Kor04 Kor05

Arc11 Arc12 0 2 4 6 8 10 12 14 Ara01 Ara02 Ara03 Ara04 Ara05 Ara06 Ara07 Ara08 Arc01 Arc02 Arc03 Arc04 Arc05 Arc06 Arc07 Arc08 Arc09 Arc10 Kor01 Kor02 Kor03 Kor04 Kor05 Kor06 Kor07 Kor08 Kor09 Kor10 Arc11 Arc12 P_par Ara01 Ara02 Ara03 Ara04 Ara05 Ara06 Ara07 Ara08 Arc01 Arc02 Arc03 Arc04 Arc05 Arc06 Arc07 Arc08 Arc09 Arc10 Kor01 Kor02 Kor03 Kor04 Kor05 Kor06 Kor07 Kor08 Kor09 Kor10 Arc11 Arc12 P_par 0

20, JUNG ET AL.

temporal distribution of Helicostomella species (Ciliophora, Tintinnina). Journal of Eukaryotic Microbiology,59,[351][352][353][354][355][356][357][358] Yi, Z., Chen, Z., Warren, A., Roberts, D., Al-Rasheid, K., Miao, M., Gao, S., Shao, C., & Song, W. (2008). Molecular phylogeny of Pseudokeronopsis (Protozoa, Ciliophora, Urostylida), with reconsideration of three closely related species at inter-and intra-specific levels inferred from the small subunit ribosomal RNA gene and the ITS1-5.8 S-ITS2 region sequences. Journal of Zoology,275,[268][269][270][271][272][273][274][275] Yu, G., Smith, D. K., Zhu, H., Guan, Y., & Lam, T. T.-Y. (2017). ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution, 8, 28-36 Zhang, W., Feng, M., Yu, Y., Zhang, C., & Xiao, T. (2012). An illustrated guide to contemporary tintinnids in the world. Science Press, Beijing Zhao, Y., Gentekaki, E., Yi, Z., & Lin, X. (2013). Genetic differentiation of the mitochondrial cytochrome oxidase c subunit I gene in genus Paramecium (Protista, Ciliophora). PLoS ONE, 8, e77044