Manas Mejari

Mihály Petreczky

Consistent and computationally efficient estimation for stochastic LPV state-space models: realization based approach

The article presents an identification algorithm for stochastic Linear Parameter-Varying State-Space Affine (LPV-SSA) representations, where the dependency of state-space matrices on scheduling signals is affine. Based on stochastic realization theory, a computationally efficient and statistically consistent identification algorithm is proposed to estimate the LPV model matrices, which are computed from the empirical covariance matrices of outputs and scheduling signal observations. The effectiveness of the proposed realization algorithm is shown via a numerical case study.

I. INTRODUCTION

The Linear Parameter-Varying (LPV) modeling paradigm can be viewed as a natural extension of Linear Time-Invariant (LTI) model class, having linear dynamic relation between input and output signals. However, unlike LTI systems, this relation can change over time according to a measurable time-varying signal, the so called scheduling variable. By the virtue of scheduling variables, non-linear and time-varying dynamics can be embedded via scheduling variable dependent LPV model coefficients. In this way, LPV models can describe the behavior of many time-varying and non-linear systems such as automobiles [START_REF] Cerone | Set-membership LPV model identification of vehicle lateral dynamics[END_REF], [START_REF] Novara | Direct identification of optimal SM-LPV filters and application to vehicle yaw rate estimation[END_REF], aircrafts [START_REF] Marcos | Development of linear-parameter-varying models for aircraft[END_REF], and distillation columns [START_REF] Bachnas | A review on datadriven linear parameter-varying modeling approaches: A high-purity distillation column case study[END_REF] etc. Naturally, identification and control of LPV models has attracted significant attention.

Over the past few years, many methods have been proposed for the identification of LPV models, in input-output (IO) [START_REF] Bamieh | Identification of linear parametervarying models[END_REF], [START_REF] Mejari | A bias-correction method for closed-loop identification of Linear Parameter-Varying systems[END_REF], [START_REF] Piga | LPV system identification under noise corrupted scheduling and output signal observations[END_REF] as well as state-space (SS) representations [START_REF] Felici | Subspace identification of MIMO LPV systems using a periodic scheduling sequence[END_REF], [START_REF] Tanelli | Identification of LPV state space models for autonomic web service systems[END_REF], [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF], [START_REF] Verdult | Kernel methods for subspace identification of multivariable LPV and bilinear systems[END_REF], [START_REF] Tóth | Modeling and identification of linear parameter-varying systems[END_REF]. From a control perspective, the IO models are not well suited for controller synthesis. The controller design approaches (e.g., [START_REF] Scherer | Mixed H 2 /H∞ control for time-varying and linear parametrically-varying systems[END_REF], [START_REF] Zhou | Essentials of Robust Control[END_REF]) often require the LPV models to be in SS representation with an affine dependency on the scheduling signal. To this end, realization theory of LPV models plays a key role. Realization theory aims at understanding the conditions under which the observed behavior of a system can be realized by a state-space representation. It also allows to formulate identification algorithms for estimating (preferably minimal) state-space representation from a finite set of observations. In fact, this paper formulates such an algorithm. The realization theory for deterministic Linear Parameter-Varying State-Space Affine (LPV-SSA) representation has been developed in [START_REF] Petreczky | Realization theory for LPV state-space representations with affine dependence[END_REF]. The results of [START_REF] Petreczky | Realization theory for LPV state-space representations with affine dependence[END_REF] have been used for model reduction in [START_REF] Batu | Moment matching based model reduction for LPV state-space models[END_REF] and to derive LPV-SS identification algorithm in [START_REF] Cox | Estimation of LPV-SS models with static dependency using correlation analysis[END_REF]. Recently, a maximum likelihood based LPV-SS identification algorithm has been proposed in [START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF] which combines LPV-IO identification methods with an LPV-IO to LPV-SS realization scheme. These methods are focused mainly on deterministic realizations, which for certain control and filtering problems are too restrictive. This shortfall can be overcome by modeling the behavior of the system using stochastic state-space representations.

Intuitively, in a stochastic setting, the LPV-SSA representation can be decomposed into deterministic and stochastic parts. Hence, intuitively, identification of stochastic LPV-SSA representations can be decomposed into two independent problems: identification of the deterministic part, and identification of the stochastic part. As it was pointed out above, there are provenly correct algorithms for identifying the deterministic part [START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF]. In this paper, we will focus on the identification of the stochastic part, i.e., on the identification of stochastic LPV-SSA representations without inputs, which only depend on noise processes. The contribution of this paper can be seen as the first step towards identification of the more general class of stochastic LPV-SSA representation with external input processes. Note that the argument of the decomposition above is not a mathematically strict one. Making this argument strict calls for developing realization theory of stochastic LPV representations with inputs. The latter was highly non-trivial even for linear time-invariant systems [START_REF] Lindquist | Linear Stochastic Systems: A Geometric Approach to Modeling, Estimation and Identification[END_REF], and we leave it for future research.

The identification algorithm presented in this paper is a non-parametric one and it is based on stochastic covariance realization algorithm for LPV-SSA representations [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF]. The main idea is to estimate the unknown model matrices by applying a modification of the covariance realization algorithm [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF] to empirical covariances computed from the observed dataset. Furthermore, in order to reduce the computational complexity, the idea of selections from [START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF] is adapted to the stochastic realization algorithm [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF]: instead of using the full Hankel-matrix, a suitable sub-matrix of the full Hankel matrix is selected for computing the estimates of system matrices and noise covariances. The proposed identification algorithm is both provenly consistent and computationally efficient.

Subspace based identification of stochastic LPV-SSA representations has a rich literature, e.g., [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF], [START_REF] Santos | Identification of bilinear systems with white noise inputs: An iterative deterministicstochastic subspace approach[END_REF], [START_REF] Favoreel | Subspace identification of bilinear systems subject to white inputs[END_REF]. The algorithm proposed in this paper differs significantly from [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF], [START_REF] Santos | Identification of bilinear systems with white noise inputs: An iterative deterministicstochastic subspace approach[END_REF], [START_REF] Favoreel | Subspace identification of bilinear systems subject to white inputs[END_REF]. In contrast to subspace identification algorithms, the proposed algorithm is provenly consistent and it does not require the observability assumption of [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF], [START_REF] Santos | Identification of bilinear systems with white noise inputs: An iterative deterministicstochastic subspace approach[END_REF], [START_REF] Favoreel | Subspace identification of bilinear systems subject to white inputs[END_REF]. However, this comes at a price, as the proposed algorithm is provenly consistent only for a very specific class of scheduling signals and stochastic LPV-SSA representations. The proposed algorithm avoids the curse of dimensionality, but this comes at the price of either using some prior knowledge on the system to determine the correct selection of the rows and columns of a Hankel-matrix or using an exhaustive search to find such a selection.

Outline: The terminology and notation used throughout the paper is introduced in Section ??. The identification problem for stochastic LPV-SSA representations is formalized in Section II. A computationally efficient stochastic realization algorithm is formulated in Section III, based on which an identification algorithm is developed in Section IV. A numerical example is reported in Section V to asses the performance of the proposed algorithm. Finally, concluding remarks are given in Section VI.

II. PROBLEM FORMULATION

We consider the following, discrete-time Linear Parameter-Varying State-Space Affine (LPV-SSA) representation

x(t + 1) = A(µ(t))x(t) + K(µ(t))v(t), ȳ(t) = Cx(t) + Dv(t), (1a)
where x, ȳ and v denote state process, output process and noise process taking values in R n , R p , and R m respectively, µ denotes the scheduling signal process taking values in R d , and at every time instance t ∈ Z,

µ(t) = [1, µ 2 (t), • • • , µ d (t)]
. The matrices A(•) and K(•) are affine functions of the scheduling signal µ(t) given by,

A(µ(t)) = A 1 + d i=2 A i µ i (t), K(µ(t)) = K 1 + d i=2 K i µ i (t). (1b)
The matrices A i ∈R n×n , K i ∈R n×m , ∀i=1, . . . , d, C ∈ R p×n and D ∈ R p×m are real, constant matrices. We define the noise covariance matrices Q i for each R-valued scheduling signal component µ i as,

Q i = E v(t)v (t)µ 2 i (t)
, for all i = 1, . . . , d. The LPV-SSA representation (1) is said to be a realization of the observed stochastic process y, if for all t ∈ Z, y(t) = ȳ(t).

The identification problem can be stated as follows: Problem 1: Assume that ȳ : Z→ R p is a sample path of the process y and μ : Z→R d is a sample path of the process µ. Given a dataset {ȳ(t), μ(t)} N t=1 consisting of samples of the output and scheduling process, compute from this dataset the estimates

{{ Âi , Ki , Qi } d i=1 , Ĉ, D}, such that as N → ∞, the estimated matrices {{ Âi , Ki , Qi } d i=1 , Ĉ, D} converge (in a suitable sense) to the matrices {{A i , K i , Q i } d i=1 , C
, D} of a true system, i.e., the LPV-SSA of the form [START_REF] Bachnas | A review on datadriven linear parameter-varying modeling approaches: A high-purity distillation column case study[END_REF] with

Q i = E v(t)v (t)µ 2 i (t) , i = 1, . . . , d
, is a realization of output process y. In particular, for a large enough N , the output ŷ obtained from the estimated matrices {{ Âi , Ki } d i=1 , Ĉ, D}, is close to the observed output y in a suitable sense.

A. Technical preliminaries

Let us define the set Σ = {1, . . . , d}. The set Σ will be referred to as the alphabet. The elements of Σ are termed as letters. A non empty word over Σ is a finite sequence of letters, i.e., w

= σ 1 σ 2 • • • σ k , where 0 < k ∈ Z, σ 1 , . . . , σ k ∈ Σ.
The set of all nonempty words is denoted by Σ + . We denote an empty word by . Let In this paper, we consider zero-mean independent identically distributed (i.i.d.) scheduling process µ = [1, µ 2 , . . . , µ d ] such that for Σ = {1, . . . , d}, we have µ 1 = 1, and for each i = 2, . . . , d, µ i is zero mean i.i.d. process, and for all t ∈ Z, we have E µ 2 σ (t) = p σ , ∀σ ∈ Σ. In particular, p 1 = 1.

Σ * = ∪ Σ + . The concatenation of two nonempty words v = a 1 a 2 • • • a m and w = b 1 b 2 • • • b n is defined as vw = a 1 • • • a m b 1 • • • b n for some m, n > 0. Note that if w = or v = , then v = v
Notations 1: For each word w ∈ Σ * , we define the number p w as follows: for an empty word w = , we define p = 1, and if w = vσ, then p w = p v p σ for v ∈ Σ * and σ ∈ Σ.

For every word w ∈ Σ + where w

= σ 1 σ 2 • • • σ k , k ≥ 1, σ 1 , . . . , σ k ∈ Σ, we define the input process µ w as µ w (t) = µ σ1 (t-k+1)µ σ2 (t-k+2) • • • µ σ k (t), ∀t ∈ Z. (2)
For a stochastic process y ∈ R k , we define for each w ∈ Σ + , the process z y w as follows

z y w (t) = y(t -|w|)µ w (t -1) 1 √ p w , ∀t ∈ Z. (3)
Notations 2 (Matrix Product): Consider a collection of square matrices

A σ ∈ R n×n , σ ∈ Σ. For any word w ∈ Σ + of the form w = σ 1 σ 2 • • • σ k , k > 0 and σ 1 , . . . , σ k ∈ Σ,
we define the following matrix product notation:

A w = A σ k A σ k-1 • • • A σ1 . For an empty word , A = I n .
We remark that in this paper, the notion of Zero Mean Wide Sense Stationary w.r.t. Inputs (ZMWSSI) process will be a central notion for the mathematical framework of stochastic LPV-SSA representations. We refer the reader to [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF]Definition 2] for the detailed definition of ZMWSSI processes. All the process considered in this paper will be assumed to be ZMWSSI process w.r.t. scheduling µ.

Definition 1 (Stationary LPV-SSA): The LPV-SSA representation (1) is stationary for {µ σ } σ∈Σ , if 1. The joint state and noise process

x v is ZMWSSI, [16, Definition 2]. 2. ∀w ∈ Σ + , t ∈ Z: E v(t -|w|)µ w (t -1)v (t) = 0. 3. ∀w ∈ Σ + , t ∈ Z: E x(t -|w|)µ w (t -1)v (t) = 0
and for all σ ∈ Σ, t ∈ Z: E x(t)µ 2 σ (t)v (t) = 0. 4. Stability: The eigenvalues of the matrix σ∈Σ p σ A σ ⊗ A σ are inside the open unit circle. Stationary LPV-SSA representations correspond to stationary generalized bilinear systems in the terminology of [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF], if µ σ (t) is viewed as an input of the bilinear system for all σ ∈ Σ. It can be proved (see [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF]) that if a stationary LPV-SSA is a realization of y, then y is ZMWSSI. Note that for stationary LPV-SSA representations, the state process x is uniquely determined by the system matrices and the noise process v satisfying Definition 1. This prompts us to use the following short-hand notation for stationary LPV-SSA representations.

Notations 3: A stationary LPV-SSA representation of the form (1) is identified with the tuple ({A σ , K σ } σ∈Σ , C, D, v).

B. Forward innovation form

In the sequel, we will especially be interested in LPV-SSA in the so called forward innovation form.

Notations The LPV-SSA representation of the form (1), is said be a realization of y in forward innovation form, if it is stationary, it is a realization of y and D = I p and v = e, i.e.,

x(t+1)= d i=1 µ i (t) (A i x(t)+K i e(t)) , y(t)=Cx(t)+e(t). (4)
A stationary LPV-SSA is a minimal realization of y if it has the minimal dimension among all possible stationary realizations of y. It turns out that we can restrict attention to LPV-SSA realizations in forward innovation form. To state the result formally, we need to introduce the following definitions. A process y is said to be square integrable w.r.t. {µ σ } σ∈Σ (SII process), if ∀w ∈ Σ + , the following process z y+ is square integrable. ∀t ∈ Z : z y+ (t) = y(t + |w|)µ w (t + |w| -1) 1 √ pw . Let us call a process y full rank, if its innovation process satisfies the following condition: for all σ ∈ Σ, the covariance matrix E e(t)e T (t)µ 2 σ (t) is invertible.

Theorem 1 ( [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF]): Assume that process y is SII.

• If y has a realization by stationary LPV-SSA of the form (1), then it has a minimal LPV-SSA realization ({ Âσ , Kσ } σ∈Σ , Ĉ, I, e) in forward innovation form. • The matrices { Âσ , Kσ } σ∈Σ , Ĉ of this LPV-SSA realization in forward innovation form and the covariances { Qi =E e(t)e (t)µ 2 i (t) } d i=1 of the innovation process can be computed from the matrices Problem 2: Assume that ȳ : Z → R p and μ : Z → R d are sample paths of the processes y and µ respectively. Given a dataset {ȳ(t), μ(t)} N t=1 , compute from this dataset the estimates { Âσ , Kσ , Qσ , Ĉ} σ∈Σ , such that as N → ∞, the estimated matrices {{ Âσ , Kσ , Qσ } σ∈Σ , Ĉ} converge (in a suitable sense) to the system matrices {{A σ , K σ , Q σ } σ∈Σ , C }, of a minimal realization of output process y in forward innovation form (4). In particular, for a large enough N , the output ŷ obtained from the estimated matrices {{ Âσ , Kσ } σ∈Σ , Ĉ}, is close to the observed output y in a suitable sense. Minimal LPV-SSAs in forward innovation form are unique up to isomorphism, which renders the problem formulation above well-posed. The well-posedness of a similar problem formulation was implicitly assumed in the LPV subspace identification literature.

{{A i , K i , Q i = E[v(t)v T (t)µ 2 i (t)]} d i=1 ,

III. REALIZATION ALGORITHM FOR STOCHASTIC LPV-SSA REPRESENTATION

In this section, we present a computationally efficient realization algorithm for stochastic LPV-SSA representation in forward innovation form (4). This algorithm is a modification of Algorithm 2 presented in [16, page 10], the main difference being the size of the Hankel-matrix used. In order to formulate the announced realization algorithm, we define the covariance sequence and Hankel matrices.

Let y ∈ R p be a ZMWSSI stochastic process. Let us define the covariance matrices Λ y w ∈ R p×p , w ∈ Σ + , and T y σ,σ ∈ R p×p , σ ∈ Σ, as follows Λ y w = E y(t)(z y w (t)) , T y σ,σ = E z y σ (t)(z y σ (t)) , (5) where z y w is as defined in (3). The covariance sequence of y is defined as the function Ψ y : Σ * → R p×pd , such that

Ψ y (w) = [Λ y 1w , Λ y 2w , . . . , Λ y dw ] . (6)
The covariance sequence can be seen as a collection of higher order moments which can determine the matrices of an LPV-SSA representation which realizes y. Conversely, the values of Ψ y (w) can be expressed in terms of matrices of a LPV-SSA representation (1). To this end, we define the following matrices

B σ = 1 √ p σ (A σ P σ C + K σ Q σ ) ∀σ ∈ Σ, B = [B 1 , B 2 , . . . , B d ], (7)
P σ = E x(t)x (t)µ 2 σ (t) , Q σ = E e(t)e (t)µ 2 σ (t)
Note that, by [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF] the state-covariance P σ is the unique solution of the following equation:

P σ = p σ σ1∈Σ A σ1 P σ1 A σ1 + K σ1 Q σ1 K σ1 . (8)
Thus, P σ is uniquely determined by the matrices {A σ , K σ , Q σ }. Then from [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF], the covariances ( 5), ( 6) satisfy the following equation:

Ψ y (w) = √ p w CA w B, ∀w ∈ Σ * and T y σ,σ = 1 pσ (CP σ C + Q σ ), ∀σ ∈ Σ.
Here, we used Notation 2 for A w .

Next, we introduce the notion of r-selection. Let us define the set Σ r as the set of all words w ∈ Σ * of length less than or equal to r, i.e., Σ r = {w ∈ Σ * | |w| ≤ r}. We define r-selection as a pair (α, β) such that [START_REF] Petreczky | Realization theory of stochastic jumpmarkov linear systems[END_REF]

, 1. α ⊆ Σ r × {1, 2, • • • , p} and β ⊆ Σ r × {1, 2, • • • , p • d}; 2. card(α) = card(β) = r,
where card denotes cardinality of the set. Thus, α and β are sets such that each element of these sets is a word-index pair. We fix the following ordering of α and β.

α = {(u i , k i )} r i=1 , u i ∈ Σ r , k i ∈ {1, 2, • • • , p}, β = {(v j , l j )} r j=1 , v j ∈ Σ r , l j ∈ {1, 2, • • • , p • d}, (9)
Example 1: Consider r = 2, number of outputs p = 2, and scheduling signal dimension d = 2, we have, Σ r = { , 1, 2, 11, 12, 21, 22}. Then, one of the r-selection pair (α, β) can be chosen as, for e.g.,

α = {(u 1 , k 1 ) , (u 2 , k 2 )} = {( , 1) , (11, 2)} and β = {(v 1 , l 1 ) , (v 2 , l 2 )} = {(21, 3) , (22, 4)}.
Let us now define the Hankel matrix H Ψ,α,β ∈ R n×n as follows: i, j = 1, . . . , n, the (i, j)th element of

H Ψ,α,β is of the form [H Ψ,α,β ] i,j = [Ψ y (v j u i )] ki,lj , (10)
where Ψ y is defined in ( 6) and [Ψ y (v j u i )] ki,lj denotes the entry of Ψ y (v j u i ) on the k i th row and l j th column, and (u i , k i ) ∈ α, (v j , l j ) ∈ β are as in the ordering of ( 9). Intuitively, the rows of H Ψ,α,β are indexed by word-index pair (u i , k i ) ∈ α, and similarly, the columns of Hankel matrix H Ψ,α,β are indexed by word-index pair (v j , l j ) ∈ β, and the element of H Ψ,α,β with the row indexed (u i , k i ) and column index (v j , l j ) is the (k i , l j )th entry of Ψ y (v j u i ).

In addition, we define the σ-shifted Hankel-matrix Z σ,Ψ,α,β ∈ R n×n as follows: its i, j-th entry is given by 10)); shifted Hankel-matrix {Z Ψ,σ,α,β } σ∈Σ (defined in [START_REF] Marcos | Development of linear-parameter-varying models for aircraft[END_REF]); {Ψ y (u i )} n i=1 , {Ψ y (v i )} n i=1 where {u i , v i } n i=1 as in ( 9) with r = n, and T y σ,σ (defined in ( 5)), maximum number of iterations M > 0.

[Z σ,Ψ,α,β ] i,j = [Ψ y (v j σu i )] ki,lj . (11
1. Compute Âσ ∈ R n×n as linear least-square solution of

Âσ H Ψ,α,β = 1 √ p σ Z σ,Ψ,α,β .
2. Compute Bσ ∈ R n×p as follows: Let B ∈ R n×p•d be matrix such that its i, j-th element is given by 9), r = n, and j = 1, 2, . . . , p • d Then, Bσ ∈ R n×p for each σ ∈ Σ is obtained by appropriate partitioning of B as B = B1 , B2 , . . . , Bd .

B i,j = [Ψ y ( u i )] ki,j = [Ψ y (u i )] ki,j where
(u i , k i ) ∈ α is as in (
3. Compute Ĉ ∈ R p×n as follows: Let C such that its i, j-th element is given by C

i,j = [Ψ y (v j )] i,lj = [Ψ y (v j )] i,lj
where (v j , l j ) ∈ β as in ( 9), r = n, and i = 1, 2, . . . , p. Then, Ĉ = C(H Ψ,α,β ) -1 . 4. Set P 0 σ = 0 and for i = 1, . . . , M do

P i+1 σ = σ1∈Σ p σ Âσ1 P i σ1 Â σ1 + Ki σ1 Qi σ1 ( Ki ) σ1 Qi σ = p σ T y σ,σ -Ĉ P i σ1 Ĉ Ki σ = Bσ √ p σ -Âσ P i σ Ĉ Qi σ -1 Output: Matrices ({ Âσ , Bσ , KM σ , QM σ , P M σ } σ∈Σ , Ĉ) ({ Âσ , Kσ , } σ∈Σ , Ĉ, I, e
) is a minimal LPV-SSA realization of y in forward innovation form with Qσ = E e(t)e (t)µ 2 σ (t) and the unique state x of this LPV-SSA realization satisfies Pσ = E x(t)x (t)µ 2 σ (t) , σ ∈ Σ. Proof: [Sketch of the proof of Theorem 2] The proof is a combination of the proof of the correctness of the basis reduced Ho-Kalman like realization algorithm in [START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF] and the proof of correctness of Algorithm 2 in [16, page 10].

IV. IDENTIFICATION ALGORITHM FOR STOCHASTIC LPV-SSA REPRESENTATION

In this section, we formulate a computationally efficient identification algorithm based on stochastic realization and selections, using N -length observation sequence of outputs and scheduling signals as detailed in Algorithm 2. In Algorithm 2, the following assumptions are made:

Assumption 1: (1) The process y is SII, full rank, and y has a realization by stationary stochastic LPV-SSA representation. (2) The n-selection pair (α, β) is such that rank H Ψ,α,β = n, where n is the state-space dimension of a minimal LPV-SSA realization of y. (3) The process (y, {µ w } w∈Σ + ) is ergodic and there exist sample paths ȳ : Z → R p and μ : Z → R d of the processes y and µ respectively such that {ȳ(t), {μ σ (t)} σ∈Σ } N t=1 is observed and for all v, w ∈ Σ + the following holds,

E y(t)(z y w (t)) = lim N →∞ 1 N N t=|w| ȳ(t)(z w (t)) E z y v (t)(z y w (t)) = lim N →∞ 1 N N t=max(|v|,|w|) zv (t)z w (t)
where for all } σ∈Σ , ĈN ) of forward innovation form (4), for a given number of iterations M , from the finite time series data {ȳ(t), {μ σ (t)} σ∈Σ } N t=1 . Algorithm 2 computes the approximations of true system matrices as well as noise covariances such that these approximations converge to the true values as N, M → ∞.

w = σ 1 σ 2 • • • σ r ∈ Σ + , r > 0, μw (t) = μσ1 (t -k + 1)μ σ2 (t -k + 2) • • • μσr (t) zw (t) = ȳ(t -|w|)μ w (t -1) 1 √ p w (12)

Algorithm 2 Identification of stochastic LPV-SSA models

Input: Observations sequence {ȳ(t), {μ σ (t)} σ∈Σ } N t=1 , and n-selection (α, β); maximum number of iterations M > 0.

1. Compute approximate covariances from observations:

Λ N w = 1 N N t=|w| ȳ(t)z w (t); T N σ,σ = 1 N N t=|w| zσ (t)z σ (t)
for all σ ∈ Σ and for every w ∈ Σ + , such that w = ivu or w = iv or w = iu or w = ivσu for some words v, u ∈ Σ * , i, σ ∈ Σ, (u, k) ∈ α, (v, l) ∈ β for some k = 1, . . . , p, l = 1, . . . , pd, and for all w ∈ Σ + , zw (t) is as in [START_REF] Mejari | A bias-correction method for closed-loop identification of Linear Parameter-Varying systems[END_REF]. 2. Construct the empirical covariance sequences Ψ N y (w) by replacing Λ w by Λ N w in [START_REF] Cox | Estimation of LPV-SS models with static dependency using correlation analysis[END_REF], and construct the empirical Hankel matrix H N Ψ,α,β ∈ R n×n and the empirical σ shifted Hankel matrix Z N σ,Ψ,α,β ∈ R n×n , σ ∈ Σ, by replacing Ψ y by Ψ N y in ( 10) and ( 11) respectively, i.e., for all i, j = 1, . . . , n,

Ψ N y (w) = Λ N 1w , Λ N 2w , . . . , Λ N dw H N Ψ,α,β i,j = Ψ N y (v j u i ) ki,lj Z N σ,Ψ,α,β i,j = Ψ N y (v j σu i ) ki,lj 3

V. NUMERICAL EXAMPLE

In this section, we present a numerical example to test the effectiveness of our algorithm. The proposed algorithm is compared to the Predictor Based Subspace IDentification (PBSID) algorithm reported in [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF]. We remark that the PBSID algorithm, requires the observability assumption of the first local model. The realization algorithm proposed in this contribution relaxes this assumption. However, it relies on choosing the correct selections (α, β). Specifically, in the following example, we show that the observability assumption is not satisfied by the true data generating system, yet, the proposed algorithm is able to estimate the model matrices such that the estimated output matches closely with the observed outputs. The quality of the match between estimated and true outputs is quantified on a noise-free validation data of length N val via Best Fit Rate (BFR) and Variance Accounted For (VAF) criterion defined for each output channel y i , i = 1, . . . , p,

as BFR yi = max 1 - N val t=1 (yi(t)-ŷi(t)) 2 N val t=1 (yi(t)-ȳi) 2 , 0 × 100%, VAF yi = max 1 -var(yi-ŷi)
var(yi) , 0 × 100%, where ŷi denotes the simulated one-step ahead model output and ȳi denotes the sample mean of the output over the validation set. The operator var(•) denotes the variance of its argument.

The LPV-SSA representation in forward innovation form (4) is used for data generation with following matrices: 4) by considering the scheduling signal process µ = [µ 1 µ 2 ] such that µ 1 (t) = 1 and µ 2 (t) is a white-noise process with uniform distribution U(-1.5, 1.5). This corresponds to the parameter values {p σ } σ∈{1,2} to be p 1 = E µ 2 1 (t) = 1 and p 2 = E µ 2 2 (t) = 0.75. The standard deviation of the white Gaussian noise e corrupting the training output is 1, i.e., e ∼ N (0, 1).

A 1 =   0.4 0.4 0 0 0 0 0 0 0   , A 2 =   0 0 0 0 0.4 0.4 0 0.4 0.4   , K 1 =   -0.036 0 1   , K 2 =   0 0.015 1.17
We run the Algorithm 2 with M = 50 iterations and with the following r-selection pair (α, β), with r = 3, α = {( , 1), (1, 1), (21, 1)}, β = {( , 1), (2, 1), (21, 1)}, which is used to choose corresponding entries of the Hankel matrix. The mean time taken to run the algorithm is 1.5 sec. 1The validation result using one-step ahead predicted outputs ŷ are reported in Table I. We compare our algorithm with the PBSID approach reported in [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF] and implemented using PBSID MATLAB toolbox. For PBSID, we consider past and future window length 8 with zero external input signal as well as with a white Gaussian noise external input u ∼ N (0, 1) with input matrices B 1 = B 2 = [0 0 1] . The proposed approach significantly outperforms the PBSID method. As the data generating LPV-SSA representation does not satisfy the underlying assumptions of PBSID method, namely non-observability of the first local model A 1 , the PBSID method fails.

VI. CONCLUSIONS AND FUTURE WORKS

In this paper, an identification algorithm for stochastic linear parameter-varying state-space representations has been presented. The proposed algorithm provides a computationally efficient alternative to the parametric subspace approaches which are prone to computational and dimensionality problems as the size of the dataset or the dimension of state-space grows. Moreover, the proposed algorithm is provenly consistent. As a future research direction, we would like to extend the results to a more general LPV model class with external inputs processes.

M.

  Mejari is with IDSIA Dalle Molle Institute for Artificial Intelligence SUPSI, 6928 Manno, Switzerland manas.mejari@idsia.ch M. Petreczky is with the Centre de Recherche en Informatique, Signal et Automatique de Lille, UMR CNRS 9189, Ecole Centrale de Lille, Villeneuve dAscq 59651, France : mihaly.petreczky@ec-lille.fr This work was partially funded by CPER Data project, which is cofinanced by European Union with the financial support of European Regional Development Fund (ERDF), French State and the French Region of Hauts-de-France and by CNRS project PEPS blanc BayesRealForRNN.

  and w = w, moreover, = . The length of the word w ∈ Σ * is denoted by |w|, and | | = 0. Eg: for alphabet Σ = {1, 2}, Σ * = { , 1, 2, 11, 12, 21, 22 . . .}, for the word w = 111 ∈ Σ * , |w| = 3.

  C, D} of (1).• If y is full rank, then the minimal LPV-SSA realizations of y in forward innovation form are isomorphic, i.e. if ({ Âσ , Kσ } σ∈Σ , Ĉ, I, e) and ({A σ , K σ } σ∈Σ , C, I, e) are two minimal realizations of y in forward innovation form, then there exists a non-singular square matrix T , such that T A σ T -1 = Âσ , T K σ = Kσ , σ ∈ Σ, CT -1 = Ĉ. Theorem 1 prompts us to aim at identifying the LPV-SSA representation in a forward innovation form (4).



  , C = 1 0 0 ,

  4 (Orthogonal projection E l ): We will use the orthogonal projection E l [•] from [16, Notation 3]. We define the forward innovation process e of ZMWSSI y as e(t) = y(t) -E l [y(t) | {z y w (t)} w∈Σ + ]. Informally, E l [y(t) | {z y w (t)} w∈Σ + ] denotes projection of y on its past values z y

w and e is the prediction error of the best linear predictor of y which is computed based on the products between past outputs and past scheduling signal values.

  reduced Hankel matrix H Ψ,α,β , and the matrices {Z σ,Ψ,α,β } σ∈Σ and covariances {Ψ y (u i ), Ψ y (v i )} n i=1 , {T y σ,σ } σ∈Σ . These matrices can be constructed from finitely many (O(n 2 )) values of Ψ y . Intuitively, the matrices { Âσ , Bσ } σ∈Σ , Ĉ are such that Ψ y (w) = Ĉ Âw B, B = [ B1 , B2 , . . . , Bd ], and for Kσ = lim M →∞ KM σ , ({ Âσ , Kσ } σ∈Σ , Ĉ, I, e) Assume y is SII, and y has a realization by stationary stochastic LPV-SSA representation and y is full rank. Assume that the n-selection pair (α, β) be such that rank H Ψ,α,β = n, where n is the state-space dimension of a minimal realization. Then Algorithm 1 is well-defined, and the matrices ({ Âσ , Bσ , KM Matrix computations using H Ψ,α,β and nselection Input: n-selection (α, β); Hankel matrix H Ψ,α,β (defined in (

						)
	Algorithm	1,	computes	the	matrices
	({ Âσ , Bσ , KM σ , QM σ , P M σ } σ∈Σ , Ĉ),	given	n-selection
	(α, β), is	a	minimal	realization
	of y in forward innovation form whose state	x
	satisfies Pσ = E x(t)x (t)µ 2 σ (t) = lim M →∞ Qσ =E e(t)e (t)µ 2 σ (t) =lim M →∞ QM σ .	P M σ , and
	This intuition is formalized in the following theorem.
	Theorem 2 (Correctness of Algorithm 1):

σ , QM σ , P M σ } σ∈Σ , Ĉ) returned by Algorithm 1 have the following property: Kσ = lim M →∞ KM σ , Qσ = lim M →∞ QM σ , QM σ > 0, Pσ = lim M →∞ P M σ , and Algorithm 1

  →∞ ĈN , where { Âσ , Bσ } σ∈Σ , Ĉ are the matrices computed in Steps 1,2,3 of Algorithm 1 usingH Ψ,α,β , Z σ,Ψ,α,β , {Ψ y (u i )} n i=1 , {Ψ y (v i )} n i=1 where {u i , v i } n i=1 as in (9) with r = n. It can be shown by induction on M that for all M , KM,N M σ } σ∈Σ are the matrices calculated in Step 4 of Algorithm 1. The statement of Theorem 3 is then a direct consequence of Theorem 2.

	Theorem 3 (Consistency): Under Assumption 1, the result
	of Algorithm 2 satisfies the following for all σ ∈ Σ:
	Kσ = lim M →∞	lim N →∞	KM,N σ	, Qσ = lim M →∞	lim N →∞	QM,N σ	,
	Pσ = lim M →∞	lim N →∞	P M,N σ	, Âσ = lim N →∞	ÂN σ , Ĉ = lim N →∞	ĈN ,
	and the LPV-SSA ({ Âσ , Kσ } σ , Ĉ, x, e) is a minimal real-
	ization of y in forward innovation form.
	Proof: [Sketch of the proof of Theorem 3] By the
	Part (3) of Assumption 1, lim N →∞ H N Ψ,α,β = H Ψ,α,β and
	lim N →∞ Z N σ,Ψ,α,β = Z σ,Ψ,α,β . Since ÂN σ , ĈN , BN σ are con-
	tinuous functions of the entries of H N Ψ,α,β , Z N σ,Ψ,α,β , σ ∈ Σ,
	it follows that Âσ = lim N →∞ Ĉ = lim N σ , QM,N σ , P M,N ÂN σ , Bσ = lim N →∞ are continuous functions BN σ , σ of ÂN , σ QM σ = lim N →∞ QM,N σ , P M σ = lim N →∞ P M,N , σ ∈ Σ, σ where { KM
	. Run Algorithm 1, with: (1) the empirical Hankel
	matrix H N Ψ,α,β instead of H Ψ,α,β ; (2) the empirical
	shifted Hankel matrix Z N σ,Ψ,α,β instead of Z σ,Ψ,α,β ; (3)
	the empirical {Ψ N y (u i )} n i=1 , {Ψ N y (v i )} n i=1 instead of {Ψ y (u i )} n i=1 , {Ψ y (v i )} n i=1 where {u i , v i } n i=1 as in (9) with r = n; (4). the empirical covariance T N σ,σ instead of T y σ,σ , ∀σ ∈ Σ.
	Output: Estimates { ÂN

σ , BN σ , KM,N σ , QM,N σ , P M,N σ } σ∈Σ , ĈN σ , BN σ , ĈN , and in fact KM σ = lim N →∞ KM,N σ , QM σ , P

  corresponds to state-dimension n=3, output dimension p = 1, and scheduling signal dimension d = 2 with Σ = {1, 2}. Note that, the system corresponding to the first local model Ã1 = A 1 -K 1 C is not observable, i.e., rank([C (C Ã1 ) . . . (C Ãl-1 1 ) ] ) = 2 < n. Training and noise free validation output sequences of length N = 500, 000 and N val = 20000, respectively, are generated from (

			TABLE I
	BEST FIT RATE (BFR) AND VARIANCE ACCOUNTED FOR (VAF) ON A
		NOISE-FREE VALIDATION DATA
		Algorithm 2 PBSID [21] PBSID (external inputs)
	BFR	93.98 %	1.58 %	0 %
	VAF	99.74 %	3.16 %	0 %

which

All computations are carried out on an i5 1.8-GHz Intel core processor with 8 GB of RAM running MATLAB R2018a.


S. J. Conway , CNRS Laboratoire de Planétologie et

 Géodynmaique, Nantes, France. (susan.conway@univnantes.fr), C. Morino, EDYTEM, Chambéry, France. J. Wright, School of Physical Sciences, Open University. F.E.G. Butcher, Department of Geography, University of Sheffield. J. Collins-May, Geography Department Newcastle University.

Introduction: Molards are cones of debris that result from the disaggregation of ice-cemented blocks transported by mass movements [1,2] (Figure 1). Recently, we linked the origin of molards directly to permafrost degradation on Earth [3]. Permafrost conditions are a pre-requisite for forming the ice cement within the blocks, and its subsequent degradation triggers the mass movement leading to the formation of the final molards. In [3], we show that molards are recognisable in remote sensing data. In addition, we report that the distribution and morphology of molards give insights into landslide dynamics. This new terrestrial work led us to search for such landforms on other planets where the role of volatiles in landscape evolution is under debate. Volatiles are substances that easily change phase, for this work we refer to ices and solids (H2O, CO2, N2, S, Cl) becoming gas or liquid. Our preliminary work has revealed molard-like landforms on Mars and on Mercury related with impact craters.

Mars -Hale Crater: Candidate molards are found in the ejecta flows of the one billion year old Hale Crater, which are similar in morphology and spatial distribution to molards found in the Mount Meager debris-avalanche deposits on Earth [5] (Figure 1). Hale impact has already been proposed to have occurred into ice-rich materials [6] and hence it is not unlikely that the ejecta may also have contained ice. The presence of molards suggests that not all the ice was melted by the impact and that the ice itself was contained within the pores of the soil, ra-ther than massive sheets of glacial ice. In future work we aim to constrain both the initial ground ice content at the Hale impact, but also the conditions during ejecta emplacement.

Mercury -Caloris Ejecta: Conical mounds are found around the Caloris Basin on Mercury. We mapped their locations and found that their distribution is consistent with them being part of the Caloris ejecta blanket. Their shape and the fact that debris from their flanks buries surrounding materials (Figure 1) suggests they continued to evolve long after the basin was formed. Because their flank-slopes and shapes are similar to molards on Earth we suggest a similar formation via volatile loss. In the last few years Mercury has been revealed to be much more volatile-rich than previously thought [e.g. 4], and our discovery of candidate molards in the ejecta of the Caloris Basin suggests that the crust and mantle has been enriched in volatiles since the basin formed ~3.8 billion years ago.

  Bibliography: [1] D.M. Cruden (1982) Can. J. Earth Sci., 19, 975-981. [2] J.P. Milana (2015) Permafr. Periglac. Process. 27, 271-284. [3] C. Morino et al. (2019) Earth Planet. Sci. Lett., 516, 136-147. [4] D.T. Blewett et al. (2011) Science, 333, 1856-1859. [5] G. Roberti et al. (2017) Geosphere, 13, 369-390. [6] A.P. Jones et al. (2011) Icarus, 211, 259-272. Acknowledgments: SJC is grateful for the financial support of CNES in support of her HiRISE and BepiColombo work. We acknowledge funding by the ANR for the Permolards project. This work acknowledges the financial support from Région Pays de la Loire, project GeoPlaNet (convention N° 2016-10982). JW is funded by Planmap (EU Horizon 2020 grant no. 776276).

Figure 1 .

 1 Figure 1. Left to right: A molard in the Móafellshyrna landslide, which occurred in 2012, height ~20 m. Planview comparison between candidate molards in Hale crater ejecta, and molards in the landslide deposits at the foot of Mount Meager, where red arrows indicate similar features (flow boundary and molards embedded in flow deposits). Conical knob on Mercury where the yellow arrow indicates crater on Caloris Planitia partially obscured by knob material and a subtle notch in the knob is visible above this crater


Gabrielle Woronoff

Philippe Nghe

Jean Baudry

Laurent Boitard

Erez Braun

Andrew D Griffiths
email: andrew.griffiths@espci.fr

Jérôme Bibette
email: jerome.bibette@espci.fr

Metabolic Cost of Rapid Adaptation of Single Yeast Cells

Cells can rapidly adapt to changing environments through non-genetic processes; however, the metabolic cost of such adaptation has never been considered. Here we demonstrate metabolic coupling in a remarkable, rapid adaptation process (10 -3 cells/hour) by simultaneously measuring metabolism and division of thousands of individual Saccharomyces cerevisiae cells using a droplet microfluidic system. Following a severe challenge, most cells, while not dividing, continue to metabolize, displaying a remarkably wide diversity of metabolic trajectories from which adaptation events can be anticipated. Adaptation requires the consumption of a characteristic amount of energy, indicating that it is an active process. The demonstration that metabolic trajectories predict a priori adaptation events provides the first evidence of tight energetic coupling between metabolism and regulatory reorganization in adaptation.

One Sentence Summary: Demonstration of the tight coupling between metabolic activity and regulatory processes during rapid adaptation at the single-cell level.

the histidine biosynthesis system from its native regulatory system and placed it under the control of the GAL system, which is highly induced in the presence of galactose and strongly repressed in glucose [START_REF] Stolovicki | Synthetic gene recruitment reveals adaptive reprogramming of gene regulation in yeast[END_REF]. Switching cells from galactose to glucose in a medium lacking histidine presents the yeast with the challenge of reinitiating histidine synthesis in order to resume growth and prevent extinction. The population dynamics of this system have previously been studied in detail [START_REF] Stolovicki | Synthetic gene recruitment reveals adaptive reprogramming of gene regulation in yeast[END_REF][START_REF] Moore | Population dynamics of metastable growth-rate phenotypes[END_REF][START_REF] Moore | Induced mutations in yeast cell populations adapting to an unforeseen challenge[END_REF][START_REF] Stern | Genome-wide transcriptional plasticity underlies cellular adaptation to novel challenge[END_REF][START_REF] Stolovicki | Collective dynamics of gene expression in cell populations[END_REF][START_REF] David | Inherited adaptation of genome-rewired cells in response to a challenging environment[END_REF][START_REF] Braun | The unforeseen challenge: from genotype-to-phenotype in cell populations[END_REF]. After switching to glucose, growth continues for ~1 day (Phase I) then, after a few days with little or no growth (Phase II), normal growth is resumed at the population level (Phase III) with a doubling time similar to that before the challenge, indicating that after only a few days the population is fully adapted (Fig. 1).

To track adaptation at the single-cell level, we used a droplet microfluidic system that allows simultaneous measurements of growth and metabolism of several thousands of single yeast cells over time [START_REF] Boitard | Monitoring single-cell bioenergetics via the coarsening of emulsion droplets[END_REF]. Individual cells, harvested from batch cultures 3 h after the beginning of Phase II (Fig. 1; green triangle), were compartmentalized in 30 pL aqueous droplets in an inert carrier oil that were immobilized in a 2D array in a closed glass observation chamber and incubated at 30°C for 65 to 70 h (Fig. 2A). A small fraction of the droplets contained single cells (typically 6%), and these were surrounded by empty droplets. The consumption of nutrients (glucose) in a cell-containing droplet creates an osmotic imbalance, resulting in water flux, inducing the shrinkage of the droplet and the swelling of neighboring cell-free droplets (Movie S1). Images were taken every 20 or 30 minutes over 3 days and used to determine the change in volume of each droplet as a function of time, reflecting cell metabolism [START_REF] Boitard | Monitoring single-cell bioenergetics via the coarsening of emulsion droplets[END_REF], while simultaneously imaging the cells to measure division (Materials and Methods). Three independent batch cultures were analyzed, however, only the consolidated data from these three independent experiments is plotted and analyzed (see Table S1 for analysis of individual experiments).

Metabolic activity was detectable for 88% of cells (final droplet volume <0.95 of initial volume, see Materials and Methods, Fig. S1). The volume variation of droplets containing cells displaying metabolic activity was plotted for both the control experiment, with histidine in the medium (Fig. 2B), and the adaptation experiment, in the absence of histidine (Fig. 2C). Cells in medium lacking histidine showed remarkably diverse metabolic trajectories (CV = 0.60 at 10 hours), while in the presence of histidine the diversity of the metabolic trajectories was more restricted (CV = 0.24 at 10 hours, Fig. S2) and similar to wild-type yeast [START_REF] Boitard | Monitoring single-cell bioenergetics via the coarsening of emulsion droplets[END_REF].

Analysis of the change of droplet volumes over time (Fig. S1-5) revealed 3 different classes of metabolizing cells (Fig. 2D, Materials and Methods): (i) steadily metabolizing cells (35% of the total population) (Fig. 2E); (ii) cells that underwent metabolic arrest (49% of the total population), the volume curve flattening before nutrient exhaustion, i.e. at a value larger than the final volume of cells grown with histidine (Fig. 2F,G), and (iii) cells that underwent metabolic recovery (4% of the total population), the shrinkage rate increasing to a similar level as for cells under non-stressed conditions (in the presence of histidine) (Fig. 2H,I). The time of metabolic arrest (T arr ) and time of metabolic recovery (T rec ) were determined as the minimum and the maximum of the second derivative of the drop-shrinkage curves, respectively (arrowheads in Fig. 2D, Materials and Methods).

We next addressed the question of the correlation between metabolic profiles and adaptation -the resumption of division. Each metabolizing cell within the shrinking drops was visualized to determine if division occurred at least once and at what time; the time of division T div being defined as the time of appearance of a first bud eventually leading to a division event.

We found that within the subset of cells showing steady metabolism, only 9% of the cells resumed division, whereas within the metabolic recovery class, 73% resumed division (Fig. 3A,B). Plotting the distribution of the time difference between the onset of cell division (T div ) and metabolic recovery (T rec ) shows that the two events are strongly correlated (R 2 =0.79, p<10 -5 , Fig. S6), with a mean value T div -T rec ≈ -1.35 ± 5.5 h (Fig. 3C), which is significantly smaller than the average time of recovery (29 h). This indicates that division starts when the maximum metabolic acceleration is reached.

The cumulated fraction of metabolic recoveries (Fig. 4A, green) and cells re-commencing division (Fig. 4A, red) both displayed a sigmoidal shape, with a quadratic increase up to 30 hours, followed by a leveling off. The instantaneous rate of adaptation also increased up to 30 h, reaching a maximum of ~10 -3 cells/h (Fig. 4B). Plotting the recovery time (T rec ) as a function of the inverse initial metabolic rate (R 0 -1 ) revealed a clear linear correlation (R 2 =0.33, p<10 -3 , Fig. 4C). The droplet volume on recovery peaks at ~74% of the initial droplet volume (Fig. 4D, Table S1), indicating that, on average, a characteristic amount of glucose (156 pg) is required for a cell to adapt. The distribution of adaptation times is explained by the distribution of the initial metabolic rates which peaks at low values (Fig. 4E, Materials and Methods). In contrast, the cumulated fraction of metabolic arrests increases quadratically throughout the course of the experiments (Fig. 4F), leading to an instantaneous death rate that continuously increases with time (Fig. 4G). The time of metabolic arrest (T arr ) is, however, poorly correlated with the initial metabolic rate (R 2 =0.12, p<10 -3 , Fig. 4H), in agreement with an age driven mechanism essentially independent of the initial metabolic rate.

We finally examined the metabolism of single cells in the adapted populations taken 2, 4 and 7 days after the batch population entered Phase III (Fig. 1; blue arrows). The overall dynamics of the adapted population (Fig. S7A), is highly comparable with cells in the absence of stress (Fig. 2B), with a similar spread of final droplet volumes (Fig. S7B) and metabolic rates (Fig. S7C,CV=0.35 for adapted compared to CV=0.24 for non-stressed cells). None of the 174 dividing cells analyzed underwent metabolic arrest, demonstrating stability of the adapted state over days, consistent with results at the population level [START_REF] Stolovicki | Synthetic gene recruitment reveals adaptive reprogramming of gene regulation in yeast[END_REF][START_REF] David | Inherited adaptation of genome-rewired cells in response to a challenging environment[END_REF].

What mechanism might explain this rapid adaptation? Consistent with previous results [START_REF] Stern | Genome-wide transcriptional plasticity underlies cellular adaptation to novel challenge[END_REF], the adaptation process observed here appears to differ from known S. cerevisiae stress responses, including the response to amino-acid starvation [START_REF] Natarajan | Transcriptional Profiling Shows that Gcn4p Is a Master Regulator of Gene Expression during Amino Acid Starvation in Yeast[END_REF], in a number of ways. Known stress responses are rapid and transient (with changes in gene expression typically peaking after <2 hours) [START_REF] Gasch | Genomic expression programs in the response of yeast cells to environmental changes[END_REF][START_REF] Zhang | High-throughput single-cell analysis for the proteomic dynamics study of the yeast osmotic stress response[END_REF], and all cells respond broadly similarly [START_REF] Zhang | High-throughput single-cell analysis for the proteomic dynamics study of the yeast osmotic stress response[END_REF]. In contrast, the characteristic timescale of adaptation of single cells that we observe is much longer (~10 3 h, see Fig. 4B) and only a fraction of cells adapt. However, the rate of adaptation (~10 -3 h -1 ) is also several orders of magnitude higher than rates reported based on genetic mechanisms [START_REF] Zeyl | Capturing the adaptive mutation in yeast[END_REF][START_REF] Hall | Selection-induced mutations occur in yeast[END_REF][START_REF] Kondrashov | Measurements of spontaneous rates of mutations in the recent past and the near future[END_REF], for example, the spontaneous reversion rates by point mutation in quiescent histidine auxotrophic yeast cells in the absence of histidine is only 10 -5 to 10 -9 h -1 [START_REF] Hall | Selection-induced mutations occur in yeast[END_REF]. This may point to an epigenetic mechanism, as recently suggested for a similar engineered yeast system [START_REF] Freddolino | Stochastic tuning of gene expression enables cellular adaptation in the absence of pre-existing regulatory circuitry[END_REF].

In conclusion, this work establishes that, in yeast, rapid adaptation necessitates consumption of a characteristic amount of energy, causing certain cells that metabolize more efficiently to adapt more rapidly. The adapted state is stable, backing up previous observations in bulk, which indicated that the adapted state is stably inherited across generations [START_REF] Stolovicki | Synthetic gene recruitment reveals adaptive reprogramming of gene regulation in yeast[END_REF][START_REF] David | Inherited adaptation of genome-rewired cells in response to a challenging environment[END_REF], showing that this is a genuine adaptation process. It is also an active process, requiring the consumption of energy, which implies exploration of different states, and fixation of the solution(s) that allow adaptation.

Finally, it is tempting to speculate that a related system may also play a role in other important processes, including cellular differentiation and development. Indeed, stochastic fluctuations at the single-cell level, have been proposed to play an important role in early stage embryonic stem cell differentiation [START_REF] Kalmar | Regulated fluctuations in nanog expression mediate cell fate decisions in embryonic stem cells[END_REF][START_REF] Macarthur | Statistical mechanics of pluripotency[END_REF][START_REF] Stumpf | Stem Cell Differentiation as a Non-Markov Stochastic Process[END_REF] and epigenetic mechanisms are increasingly believed to play an important role the development, progression and emergence of drug resistance in cancer [START_REF] Flavahan | Epigenetic plasticity and the hallmarks of cancer[END_REF][START_REF] Salgia | The Genetic/Non-genetic Duality of Drug 'Resistance' in Cancer[END_REF]. Furthermore, it may have implications for important biotechnologies, such as reprogramming of human somatic cells to generate induced Pluripotent Stem Cells (iPSCs) [START_REF] Buganim | Mechanisms and models of somatic cell reprogramming[END_REF].

Figures:

Fig. 1. Population growth in bulk cultures during adaptation. The culture of rewired cells started with galactose as the sole carbon source (first red arrow). Then, the medium was switched to glucose (second red arrow). The culture continued to grow during the "initial growth phase" (phase I), before entering the "latent phase" (phase II) where growth stopped. Finally, the population grew again, entering the "adapted phase" (phase III). Green and blue triangles indicate time points when samples were taken for single-cell measurements in droplets (green triangle: 3 hours after the beginning of phase II; blue triangles: 2, 3 and 6 days after the onset of phase III).

Fig. 2. Single-cell analysis of metabolic dynamics in droplets.

A) Time-lapse sequences of droplets containing a metabolizing dividing cell, a metabolizing non-dividing cell and an arrested cell. B) Traces of droplet volumes V normalized by the initial volume V 0 , in the absence of metabolic challenge (with histidine). C) Traces of normalized droplet volumes starting 3 hours after switching from galactose to glucose. Lines are randomly color coded. D) Examples of traces of normalized droplet volumes of each category of metabolic response. E) Subset of steady metabolisms from panel C. F) Distribution of final normalized droplet volume (V f /V 0 ) in the presence (red) and absence (blue) of histidine. V e indicates the largest final droplet volume with histidine. G) Subset of arrested metabolisms from panel C. H) Distributions of the difference ΔR (metabolic rate increase) between the initial (at 2 hours) and maximum metabolic rates in the presence (red) and absence (blue) of histidine. ΔR* indicates the cross-over between the two distributions. I) Subset of metabolic recoveries from panel C.

Materials and Methods

All chemicals were purchased from Sigma-Aldrich unless otherwise mentioned.

Strains and culture conditions. The Saccharomyces cerevisiae YPH499 strain carries a deletion of the endogenous chromosomal HIS3 gene and a plasmid containing the divergent GAL1/GAL10 promoters. There is a HIS3 gene under sole regulation of GAL1 promoter, and a GFP reporter gene under control of the GAL10 promoter (4). Cells were grown from colonies on plates in 25 mL of culture medium comprising: 1.7 g/L of yeast nitrogen base without amino acids and ammonium sulfate, 5 g/L ammonium sulfate, 1.4 g/L amino acid dropout powder (without Trp, His, Leu, Ura), 0.004 g/L L-tryptophan, and 0.002 g/L uracil, 20 g/L of galactose or glucose; and incubated at 30°C, shaking at 200 r.p.m. in 100 mL Erlenmeyer flasks. Cultures were diluted every 12h to maintain OD<1.0 (never transferring less than 3.10 6 cells). Rewired yeast cells, were first grown in a batch culture with galactose medium lacking histidine. The optical density (OD 600nm ) of the culture was monitored and it was maintained in the exponential phase of growth by serial dilutions for 2 days. The medium was then switched from galactose to glucose. As previously reported (4, 9), naïve cells (i.e. rewired cells that had never been grown in glucose before) were able to grow for about 20 hours in glucose media lacking histidine with a doubling time, t D , of 1.9 h (phase I). Following this phase, growth slowed considerably (t D 17.1 h), leaving the OD 600nm approximately constant (phase II, the "latent phase"). At the end of this phase, after durations that varied across repeated experiments from (1 day to 6 days), the population adapted and started growing and proliferating (phase III, the "adapted phase"). The doubling time of the adapted population reduced from 3.1 h 1 day after entering phase III, to 2.1 h after 8 days in phase III, reflecting the fact that the population continues to adapt after resuming growth [START_REF] Moore | Population dynamics of metastable growth-rate phenotypes[END_REF]. Three independent batch cultures were analyzed and the growth curve for one of these experiments is shown in Fig. 1.

Microfluidic device fabrication.

Microfluidic devices were obtained using conventional soft lithography methods (25), as described [START_REF] Frenz | Reliable microfluidic on-chip incubation of droplets in delay-lines[END_REF]. Molds were prepared using SU8-2015 or SU8-2075 photoresist (MicroChem Corp.) and used to pattern 20 and 75 µm-deep channels onto silicon wafers (Siltronix). The channels of the devices were passivated with Aquapel in HFE7100 (3M) and subsequently flushed with compressed nitrogen gas.

Formation and imaging of droplet arrays. Single-cell metabolic and growth dynamics of yeast cells from batch cultures were measured as in Boitard et al. [START_REF] Boitard | Monitoring single-cell bioenergetics via the coarsening of emulsion droplets[END_REF]. Yeast cells from batch culture were individually compartmentalized in monodisperse 30 pL volume (38 µm diameter) aqueous droplets containing fresh glucose medium lacking histidine by hydrodynamic flow-focusing [START_REF] Anna | Formation of dispersions using "flow focusing" in microchannels[END_REF] with a fluorinated oil phase containing fluorosurfactant [START_REF] Holtze | Biocompatible surfactants for water-in-fluorocarbon emulsions[END_REF]. For each experiment, about 200,000 droplets were incubated at 30°C in a 2D array and images were taken every 30 minutes over three days. Briefly, a flow-focusing device [START_REF] Anna | Formation of dispersions using "flow focusing" in microchannels[END_REF] was used for droplet generation and flow rates were controlled using standard-pressure infuse/withdraw PHD 22/2000 syringe pumps (Harvard Apparatus Inc., Holliston, MA). Syringes (Hamilton) connected to the microfluidic device using 0.6 x 25 mm Neolus needles (Terumo Corporation) and PTFE tubing with an internal diameter of 0.56 mm and an external diameter of 1.07 mm (Fisher Bioblock Scientific). The aqueous phase, comprising cells suspended in culture medium, was injected at 300 µL/h and dispersed in a continuous phase consisting of HFE-7500 fluorinated oil (3M) containing 2% (w/w) EA surfactant (Raindance Technologies), a PFPE-PEG-PFPE amphiphilic tri-block copolymer [START_REF] Holtze | Biocompatible surfactants for water-in-fluorocarbon emulsions[END_REF], injected at 150 µL/h, forming 30 pL volume (38 µm diameter) droplets at 2,800 droplets s -1 . Droplets were produced in a compact manner and were directly incubated in a glass chamber of 3.5 x 1.5 cm to form a compact 2D droplet array [START_REF] Boitard | Monitoring single-cell bioenergetics via the coarsening of emulsion droplets[END_REF]. Compaction of the emulsion prevents droplet movement to enable their tracking over long time scales. The chamber was maintained at 30°C on a Nikon T300 inverted microscope with a Thorlabs MAX202 XY stage. Images of the droplet array were taken every 30 min using a Hamamatsu Orca-ER camera. Custom-made Labview software was used to automate image acquisition and microscope control.

Data processing. For a schematic representation see Fig. S5.

Droplet tracking and normalization. Droplets were detected with an in-house Matlab segmentation routine and tracked over time by a nearest neighbour criterion, knowing that displacements from image to image are much smaller than the characteristic droplet diameter. Images in which large displacements occured were easily detected as they displayed strongly discontinuous volume traces and were discarded. In order to minimize volume fluctuations caused by defocusing in time, volumes were normalized by the time course of the average of more than 10 empty droplets for each image. Images for which normalization failed were discarded (Fig. S5). The first and second derivatives of the normalized volume over time were computed over 5 time point sliding windows and provide the metabolic rate and metabolic acceleration and were used to automatically determine the maximum metabolic rate, time of metabolic acceleration and time of metabolic arrest (Fig. S3).

Significantly shrinking droplets.

A droplet was considered to be significantly shrinking (i.e. the compartmentalized cell showed detectable metabolism) if its final normalized volume V end was smaller than V s =1-3σ, where σ is the coefficient of variation of empty droplet final volumes taken over all experiments. We measured V s = 0.95. Shrinking but empty droplets that passed this filter were removed by visual inspection of the droplet content.

Poisson parameter and initial fraction of dead cells. The Poisson parameter λ (the mean number of cells per droplet) was determined from the frequency of occupied droplets on a sample of 10 initial images (~1500 droplets) per experiment. The number of droplets containing cells arrested from the start was estimated by calculating the difference between λ and the number of significantly shrinking droplets (Fig. S1).

Metabolic arrest.

Metabolically arrested cells were defined as cells: (i) which do not consume all the nutrients available in the droplets within the timeframe of the experiment, as defined by V end > V e , where V end is the final droplet volume and V e = 0.55 is the maximal normalized volume reached by cells in the control experiment (with histidine), and (ii) for which the volume trace significantly flattens during the time course of the experiment. The latter was determined as the existence of an inflection in the volume curve leading to reach a rate R < R max /10, where R max is the maximum rate observed for this cell. The metabolic arrest point is defined as the time at which R reaches R max /10 (Fig. S3). In the presence of histidine, an abrupt decrease in metabolic rate was reached when glucose in the droplet was exhausted by the cells [START_REF] Boitard | Monitoring single-cell bioenergetics via the coarsening of emulsion droplets[END_REF], corresponding to a mean final droplet volume of 0.49±0.002 of the initial value, and always <0.55 of the initial value (Fig. 2F,red). In contrast, in the adaptation experiment in the absence of histidine (Fig. 2F, blue), 54% of cells had a final volume >0.55 of the initial value, indicating that they had not consumed all the glucose in the droplet. These cells fall into two categories. In the first category are cells (49%) that arrested their metabolism during the experiment (Fig. 2G), as they were initially metabolically active and showed an abrupt metabolic deceleration and a stable volume, indicating that the cells were not metabolically active from that point in time (T arr ). A second category of cells (5%) neither arrested, nor recovered, but had a metabolism slow enough that glucose was not exhausted at the end of the experiment.

Metabolic recovery. Metabolic recoveries in the adaptation experiment were determined as those cells whose metabolic rate increased to levels comparable to the control with histidine, as computed by the difference (ΔR) between the maximum rate (R max ) and the initial rate (R 0 ) (2 hours after compartmentalization) obtained from the volume time course of each droplet (Fig. 2C). For this we determined the lower and upper bounds of the crossing values for the histograms of ΔR (Fig. 2H) for different binning intervals in the presence and absence of histidine (Fig. S4). The threshold value ΔR* = 0.013 h -1 , which indicates the cross-over between the two distributions, was obtained when these bounds converged. Cells with ΔR > ΔR* were classed as metabolic recoveries. In the control, with histidine, the metabolic rate of cells steadily increased due to cell growth and division before glucose exhaustion (Fig. 2B) [START_REF] Boitard | Monitoring single-cell bioenergetics via the coarsening of emulsion droplets[END_REF], with a mean rate increase of 3.6±0.7 10 -2 h -1 (Fig. 2H,red). In the adaptation experiment, 4% of cells accelerated their metabolic rates (N=120, Fig. 2H, blue) to values similar to those observed in the control with histidine (>1.3 10 -2 h -1 ).

Computation of instantaneous rates of arrest and metabolic recovery. Three distinct cell fates were categorized from the metabolic curves: metabolic arrest, metabolic recovery and steady metabolism. Calling N the total number of cells observed during the experiment, N a and N r the respective cumulative number of arrested and recovering cells, the instantaneous rate of arrest and recovery were respectively computed as 𝑑𝑁 ! /𝑑𝑡 (𝑁 -𝑁 ! -𝑁 ! ) and 𝑑𝑁 ! /𝑑𝑡 (𝑁 -𝑁 ! -𝑁 ! ).

Effective coefficients of determinations. The correlations of Fig. 4C,H display regions of experimentally inaccessible parameter values, due to the finite amount of resource available in each droplets and the noise in droplet volume measurements. The coefficient of determination R a 2 computed directly from the data thus comprises a contribution of these boundary constraints. This contribution R b 2 to the total covariance was estimated using a Monte Carlo approach preserving the marginal distributions: 10 5 random permutations of the measured values were generated under the boundary constraints, taking R b 2 as the mean coefficient of determination of the last 10 4 realizations. Assuming independence between the experimental constraints and cell fates in the absence of resource exhaustion, the additive contribution of these phenomena to the covariance leads to effective coefficients of determinations The metabolic increase ΔR is the difference between the initial metabolic rate (taken at 2 hours) and the maximum measured metabolic rate R max (see Fig. S3). For each bin size, we determined the lower and upper bound for the cross-over between the metabolic increase distributions (see Fig. 2H), in the presence or absence of metabolic challenge. These bounds converge at a certain bin value which optimally defines the cross-over ΔR* between these two distributions.

R 2 = R b 2 -R a 2 .

Fig. S5. Classification tree of metabolic responses.

For the details of the treatment, refer to the Materials and Methods. V end is the final droplet volume; V s =1-3σ, where σ is the coefficient of variation of empty droplet final volumes taken over all experiments; V e is the maximal normalized volume reached by cells in the control experiment (with histidine); λ is the Poisson parameter corresponding to the average number of cells encapsulated per droplet; R is the metabolic rate; R max is the maximum metabolic rate reach for each cell; T arr is the time at which R reaches R max /10 in arresting cells; ΔR is the metabolic rate increase computed as the difference between R max and the initial rate R 0 (at 2 hours); ΔR* is the cross-over value between low and high metabolic rate increases (Fig. S4); A is the metabolic acceleration, which is the time derivate of R; A max is the maximum metabolic acceleration (Fig. S3); T rec is the time at which A max is reached within the population of cells that significantly increase their metabolism (ΔR > ΔR*); T div is the time at which the first bud appears in droplets where the final number of cells is larger than 2. and dividing cells (89±2%). Consistent with adaptation observed at the population level, the fraction of dividing cells is >10-fold higher than at the beginning of Phase II (5).

Fig. 3 .

 3 Fig. 3. Relation between onset of division and metabolic recovery. A) Fraction of dividing cells in steady metabolisms. B) Fraction of dividing cells in metabolic recoveries. C) Distribution of the difference between division time (T div ) and metabolic recovery time (T rec ).

Fig. 4 .

 4 Fig. 4. Characterization of adaptation. A) Cumulative fraction over time for metabolic recoveries (green) and division recoveries (red). Black lines are a quadratic fit over the first 30 hours. B) Instantaneous frequency of metabolic recoveries of cells. At each time point, frequencies are computed over 7 consecutive points and error bars are standard deviations over the values obtained over these points. C) Plot of metabolic recovery time (T rec ) as a function of inverse initial metabolic rate (R 0 -1 ). The parameter region above the gray line corresponds to cells that have exhausted droplet resources before being able to recover (Material and Methods). D) Distribution of normalized droplet volume at T rec . E) Distribution of R 0 for the subset of metabolic recoveries. F) Cumulative fraction over time for metabolic arrests (black). The red line is a quadratic fit. G) Instantaneous frequency of arrests of cells that have neither recovered nor arrested their metabolism. Frequencies and error bars were computed as in B). H) Plot of time of metabolic arrest (T arr ) as a function of R 0 -1 . The parameter region above the upper gray line corresponds to cells that have exhausted droplet resources before arresting. The parameter region below the lower gray line corresponds to curves whose final volume is indistinguishable from curves of empty droplets (overall volume variation <5%). The contribution of the constraints in panels C) and H) has been accounted for to compute the coefficients of determination R 2 (Material and Methods).

Fig. S1 .

 S1 Fig. S1. Threshold for significant shrinking. The threshold V s = 0.95 for maximal final volume indicating a significant volume decrease is defined as 3 standard deviation away from the mean of volume distribution of empty droplets.

Fig. S2 .

 S2 Fig. S2. Distribution of metabolic rates at 10 hours after encapsulation of yeast in droplets.Red: in the absence of metabolic challenge (with histidine). Blue: in response to a metabolic challenge (no histidine).

Fig. S3 .

 S3 Fig. S3. Automated analysis of metabolic kinetics from the time-lapse measurement of droplet volumes. A) Example of 20 droplet volumes V over time normalized by their initial volume, V 0 . The point of metabolic acceleration (A max ) is shown by triangles, of maximum rate (R max ) by squares and metabolic arrest (T arr ) by circles (see below). B) Metabolic rate R computed as minus the time derivative of the normalized volume. R max (square) is determined as the maximum of each curve. Metabolic arrest is determined as the point at which R reaches 1/10 of R max , and generally corresponds to the minimum of the acceleration curves in panel c, although the latter yields a less reliable estimate. C) Metabolic acceleration A computed as the derivative of the metabolic rate. The time of acceleration (triangle) is determined from the maximum A max of this curve. The time of recovery (T rec ) is the time at which A max is reached within the population of cells that significantly increase their metabolism (ΔR > ΔR*) (see Fig. 2H and Fig. S4).

Fig. S4 .

 S4 Fig.S4. Determination of the metabolic increase threshold ΔR* for significant metabolic recovery. The metabolic increase ΔR is the difference between the initial metabolic rate (taken at 2 hours) and the maximum measured metabolic rate R max (see Fig.S3). For each bin size, we determined the lower and upper bound for the cross-over between the metabolic increase distributions (see Fig.2H), in the presence or absence of metabolic challenge. These bounds converge at a certain bin value which optimally defines the cross-over ΔR* between these two distributions.

Fig. S6 .

 S6 Fig. S6. Correlation between the times of division (T div ) and metabolic recovery (T rec ).Each dot corresponds to the analysis of a single volume trace measured in the absence of histidine. The diagonal corresponds to T div = T rec .

Fig. S7 .

 S7 Fig. S7. Adapted population. A) Volume traces of droplets containing single cells taken from the culture 4 days after the challenge. Droplet volume V is normalized by the initial volume V 0 .Lines are randomly color-coded. B) Distribution of final volumes (at 70 hours) 4 days after the challenge (blue) and for the control with histidine (red). The vertical dashed line represents V e , the maximum normalized volume reached by cells in the histidine positive control. C) Distribution of metabolic rates at 10 hours, 4 days after stress (blue) and for the control with histidine (red). The variance in metabolic rate of the dividing cells is similar in the samples taken 2, 4 and 7 days after entering the adapted phase (CV = 0.35) and close to that observed for control cells switched to glucose plus histidine (CV = 0.24). D) Fraction of each cell type 2, 4 and 7 days after the challenge. The proportions of the different classes of cells are very similar at the three time points: non-metabolizing cells (1.6±1%), non-dividing metabolizing cells (10±1%) and dividing cells (89±2%). Consistent with adaptation observed at the population level, the fraction of dividing cells is >10-fold higher than at the beginning of Phase II[START_REF] Moore | Population dynamics of metastable growth-rate phenotypes[END_REF].

Acknowledgments:

The authors thank Paul M. Chaikin (New York University) for helpful

support from a NSF-BSF grant (#2014713).

Table S1. Parameters measured and computed for each individual experiment after stress, after adaptation, for the control with histidine, and over groups of experiments. In all cases, the symbol "N" refers to a number of droplets. Movie S1. Droplets were incubated at 30°C in a sealed glass chamber imaged over 70 h. Images were taken every 30 min, as indicated by the time stamp. Individual cells were compartmentalized in aqueous droplets dispersed in an oil phase and parked in a 2D array in a closed glass chamber. About 6% of the droplets contain a single cell. The consumption of nutrients (glucose) in a cell-containing droplet creates an osmotic imbalance, resulting in an osmotically driven water flux, which induces the shrinkage of this droplet and the swelling of neighboring cell-free droplets The change in volume of the droplet containing the cell is extracted by image analysis to measure cell metabolism as a function of time, while imaging enables measurement of cell division.