Neil Cammardella

Ana Bušić

Sean Meyn

Kullback-Leibler-Quadratic Optimal Control in a Stochastic Environment

This paper presents advances in Kullback-Leibler-Quadratic (KLQ) optimal control: a stochastic control framework for Markovian models. The motivation is distributed control of large networked systems. The objective function is composed of a control cost in the form of Kullback-Leibler divergence plus a quadratic cost on the sequence of marginal distributions. With this choice of objective function, the optimal probability distribution of a population of agents over a finite time horizon is shown to be an exponential tilting of the nominal probability distribution. The same is true for the controlled transition matrices that induce the optimal probability distribution.

However, one limitation of the previous work is that randomness can only be introduced via the control policy; all uncontrolled processes must be modeled as deterministic to render them immutable under an exponential tilting. In this work, only the controlled dynamics are subject to tilting, allowing for more general probabilistic models.

Numerical experiments are conducted in the context of power networks. The distributed control techniques described in this paper can transform a large collection of flexible loads into a 'virtual battery' capable of delivering the same grid services as traditional batteries. Additionally, quality of service to the load owner is guaranteed, privacy is preserved, and computation and communication requirements are reduced, relative to alternative centralized control techniques.

I. INTRODUCTION

The setting of this paper is optimal control of Markov Decision Processes (MDPs). The state space S and input space U are assumed to be finite. A finite time horizon is considered, indexed by {k : 1 ≤ k ≤ K}. The controlled transition matrix T k defines the statistics of the state process S with input process U :

T k (x, s ′ ) = P{S k+1 = s ′ | S i , U i , 0 ≤ i ≤ k; S k = s, U k = u}
The policies {ϕ k } are assumed to be Markovian:

ϕ k (u | s) = P{U k = u | S i , U i , 0 ≤ i < k; S k = s} (1)
As in [START_REF] Cammardella | Kullback-Leibler-Quadratic optimal control of flexible power demand[END_REF], [START_REF] Cammardella | Simultaneous allocation and control of distributed energy resources via Kullback-Leibler-Quadratic optimal control[END_REF], the Kullback-Leibler-Quadratic (KLQ) optimization criterion is based on convex functions of the marginal probability mass functions (pmfs) of the joint stateinput process X k = (S k , U k ):

ν k (x) = P{S k = s, U k = u} , x = (s, u) ∈ X (2)
Funding from the National Science Foundation under award EPCN 1935389 and French National Research Agency grant ANR-16-CE05-0008 is gratefully acknowledged.

where X = S × U is the state-input space.

A sequence {C k } of cost functions is given, and the control objective is to obtain the solution to the optimization problem

J * (ν 0 ) = min U K k=1 C k (ν k ) (3)
where the minimum is over all randomized policies. Two classes of constraints are imposed. First, the initial pmf ν 0 for X 0 is fixed. Second are the dynamics, which can be expressed as a sequence of linear constraints on the marginals:

u ′ ν k (s ′ , u ′ ) = s,u ν k-1 (s, u)T k (x, s ′ ) , s ′ ∈ S (4)
One standard application is variance penalized MDPs. Our motivation is applications to mean field control as an approach to distributed control

A. Distributed control

The general optimization problem (3) falls outside of textbook stochastic control problems. It is inspired by meanfield game theory [START_REF] Lasry | Mean field games[END_REF], [START_REF] Huang | Large-population cost-coupled LQG problems with nonuniform agents: Individual-mass behavior and decentralized ε-Nash equilibria[END_REF], [START_REF] Huang | Large population stochastic dynamic games: closed-loop Mckean-Vlasov systems and the Nash certainty equivalence principle[END_REF], [START_REF] Caines | Mean field games[END_REF], [START_REF] Guéant | Mean Field Games and Applications[END_REF], [8] (see [START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications I[END_REF], [START_REF]Probabilistic Theory of Mean Field Games with Applications II: Mean Field Games with Common Noise and Master Equations, ser. Probability Theory and Stochastic Modelling[END_REF] for recent surveys), and motivated in particular by applications to distributed control [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF].

The control objective emerges from the approximation of a particular distributed control problem: a central authority wishes to shape the aggregate behavior of N ≫ 1 homogeneous agents, each modeled by the transition kernel T k , with state-input denoted {X i k = (S i k , U i k ) : 1 ≤ i ≤ N }. The sequence of empirical distributions is denoted

ν N k (x) = 1 N N i=1 I{S i k = s , U i k = u}, x = (s, u) ∈ X
The optimization criterion of interest is (3), but with ν k replaced by ν N k . A mean-field control approximation is justified by applying the law of large numbers: fix a sequence of randomized policies {ϕ k : 0 ≤ k ≤ K}, and consider N as a variable. The empirical distribution {ν N k } converges as N → ∞ for each k, and the limit satisfies the linear constraints [START_REF] Huang | Large-population cost-coupled LQG problems with nonuniform agents: Individual-mass behavior and decentralized ε-Nash equilibria[END_REF].

Section IV focuses on a homogeneous population of residential refrigerators for the creation of "virtual energy storage" for power grid applications. The goal is to shape the power usage of the population of loads. Let Y : X → R + denote power consumption as a function of state, so that the average power consumption of the population tracks the reference signal r, with acceptable error. The temperature of the ith load at time k is denoted Θ i k , and the power mode (0 or 1) is denoted M i k . A typical linear model is given by

Θ i k+1 = Θ i k + α[Θ a -Θ i k ] -βM i k (5)
where α > 0, β > 0 and Θ a ∈ R denotes ambient temperature. This is a (deterministic) MDP model with stateinput given by

X i k = (Θ i k , M i k ).
In this work, our control design allows for the inclusion of disturbances in the model [START_REF] Huang | Large population stochastic dynamic games: closed-loop Mckean-Vlasov systems and the Nash certainty equivalence principle[END_REF].

B. Literature review

Our primary motivation is application to distributed control of power systems, specifically Demand Dispatch. The term was introduced in the conceptual article [START_REF] Brooks | Demand dispatch[END_REF] to describe the possibility of distributed intelligence in electric loads, designed so that the population would help provide supplydemand balance in the power grid. Contributing to this science has been a focus of the authors for the past decade [START_REF] Hao | How demand response from commercial buildings will provide the regulation needs of the grid[END_REF], [START_REF] Meyn | Ancillary service to the grid from deferrable loads: The case for intelligent pool pumps in Florida[END_REF], [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF], and many others (see [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF] for a recent bibliography).

The goal in much of this prior work is to modify the behavior of loads so that their aggregate power consumption tracks the reference signal r that is broadcast by a Balancing authority (BA), based on distributed control, with local randomized decision rules. Randomized control techniques for Demand Dispatch have been proposed in [START_REF] Mathieu | State estimation and control of electric loads to manage real-time energy imbalance[END_REF], [START_REF] Tindemans | Decentralized control of thermostatic loads for flexible demand response[END_REF], [START_REF] Almassalkhi | Packetized energy management: asynchronous and anonymous coordination of thermostatically controlled loads[END_REF], [START_REF] Benenati | A tractable formulation for multi-period linearized optimal power flow in presence of thermostatically controlled loads[END_REF] based on entirely different control architectures.

The following control strategy is common to all of the approaches described in [START_REF] Meyn | Ancillary service to the grid from deferrable loads: The case for intelligent pool pumps in Florida[END_REF], [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF]. It is assumed that a family of transition matrices {P ζ : ζ ∈ R} is available at each load. A sequence {ζ 0 , ζ 1 , . . . } is broadcast from the BA, based on measurements of the grid, and at time k an individual load transitions according to this law:

P{X k+1 = x ′ | X k = x, ζ k = ζ} = P ζ (x, x ′ )
The paper [START_REF] Meyn | Ancillary service to the grid using intelligent deferrable loads[END_REF] re-interprets the control solution of [START_REF] Todorov | Linearly-solvable Markov decision problems[END_REF] as a technique to create the family {P ζ } through the solution to the nonlinear program:

P ζ := max π,P ζ⟨ π, Y ⟩ -K(P ∥P 0 ) , ζ ∈ R , (6)
where K denotes the infinite-horizon relative entropy rate for two Markov chains:

K(P ∥P 0 ) := x,x ′ π(x)P (x, x ′ ) log P (x, x ′ ) P 0 (x, x ′ ) (7)
in which π is the invariant pmf for P . The maximum in [START_REF] Caines | Mean field games[END_REF] is over all (π, P ) subject to the invariance constraint πP = π [START_REF] Meyn | Ancillary service to the grid using intelligent deferrable loads[END_REF], [START_REF] Bušić | Distributed randomized control for demand dispatch[END_REF], [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF].

The optimal control formulation here is also based on (7), with similar motivation, but the control approaches are entirely different. Further discussion is postponed to Section II-C.

Organization. The remainder of this paper is organized as follows: Section II contains a description of our control objective and a characterization of our main result. Computation of the optimal policy {ϕ * k } is described in Section III. Results from numerical experiments are collected together in Section IV. Conclusions and directions for future research are contained in Section V. The appendix contains abbreviated proofs of some of the results of this paper. Complete proofs can be found in [START_REF] Cammardella | Kullback-Leiblerquadratic optimal control[END_REF].

II. KULLBACK-LEIBLER-QUADRATIC CONTROL

A. Optimal control formulation KLQ is designed to balance two objectives: (i) ν k ∼ ν 0 k , where {ν 0 k } models nominal behavior. (ii) ⟨ν k , Y⟩ ≈ r k , where {r k } is a reference signal, and Y : X → R. A pmf p 0 on X K+1 defines the nominal model:

p 0 (⃗ x) = ν 0 0 (x 0 )P 0 0 (x 0 , x 1 )P 0 1 (x 1 , x 2 ) • • • (8)
where ⃗ x denotes the elements of X K+1 , {P 0 k } are the nominal Markov transition matrices with transition probabilities

P 0 k (x, x ′ ) := P{X k+1 = x ′ | X k = x}
, and the nominal marginal pmfs are:

ν 0 k (x k ) = xi:i̸ =k p 0 (⃗ x) , ⃗ x ∈ X K+1
The nominal transition probabilities are represented as a product of two conditional pmfs:

P 0 k (x, x ′ ) = T k (x, s ′ )ϕ 0 k+1 (u ′ | s ′ ) , x, x ′ ∈ X (9)
where {ϕ 0 k } denotes the nominal randomized policy. Also, any randomized policy {ϕ k } produces Markov transition matrices {P k } with transition probabilities

P k (x, x ′ ) = T k (x, s ′ )ϕ k+1 (u ′ | s ′ ) , x = (s, u) ∈ X (10)
The marginal pmfs evolve according to linear dynamics, which are shown to be equivalent to (4):

ν k = ν k-1 P k-1 , 1 ≤ k ≤ K (11)
where the kth marginal ν k is interpreted as a d-dimensional row vector. The two control objectives motivate the cost function considered in this paper:

C k (ν) = D(ν, ν 0 k ) + κ 2 ⟨ ν, Y ⟩ -r k 2
in which κ > 0 is a penalty parameter, and D penalizes deviation from nominal behavior. The finite-horizon optimal control problem is thus

J * (ν 0 0 ) = min K k=1 D(ν k , ν 0 k ) + κ 2 ⟨ ν k , Y ⟩ -r k 2 ( 12
)
where the initial pmf ν 0 0 is given. The relative entropy rate (7) will be adopted as the cost of deviation. Under our assumptions, this reduces to

D(ν k , ν 0 k ) := s,u ν k (s, u) log ϕ k (u | s) ϕ 0 k (u | s) (13)
The terminology is justified through the following steps. First, we have seen that any randomized policy gives rise to a pmf p that is Markovian:

p(⃗ x) = ν 0 0 (x 0 )P 0 (x 0 , x 1 )P 1 (x 1 , x 2 ) • • •
where P k is defined in [START_REF]Probabilistic Theory of Mean Field Games with Applications II: Mean Field Games with Common Noise and Master Equations, ser. Probability Theory and Stochastic Modelling[END_REF], and the initialization ν 0 0 is specified. The relative entropy is the mean log-likelihood:

D(p∥p 0 ) = L(⃗ x) p(⃗ x) (14)
where L = log(p/p 0 ) is an extended-real-valued function on X K+1 . The expression for P k in [START_REF]Probabilistic Theory of Mean Field Games with Applications II: Mean Field Games with Common Noise and Master Equations, ser. Probability Theory and Stochastic Modelling[END_REF] and the analogous formula for P 0 k using ϕ 0 k+1 gives

L(⃗ x) = log p(⃗ x) p 0 (⃗ x) = K k=1 log ϕ k (u k | s k ) ϕ 0 k (u k | s k ) (15)
Consequently, D(p∥p

0 ) = K k=1 D(ν k , ν 0 k ).
The optimal control problem ( 12), subject to the constraint [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF], can be expressed

J * (ν 0 0 ) := min ν,γ K k=1 D(ν k , ν 0 k ) + κ 2 K k=1 γ 2 k (16a) s.t. γ k + r k -⟨ ν k , Y ⟩ = 0 (16b) u ′ ν k (s ′ , u ′ ) - s,u ν k-1 (s, u)T k (x, s ′ ) = 0 (16c)
Prop. 2.1 asserts that the objective function is convex. It is also evident that the constraints (16b) and (16c) are linear in ν k ; hence, the optimization problem ( 16) is convex. Proposition 2.1: The optimization problem ( 16) is jointly convex in {ν k , γ k : 1 ≤ k ≤ K}. Furthermore, the constraint (16c) is equivalent to [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF].

⊓ ⊔

An abbreviated proof can be found in Appendix A.

B. Dead-beat control

In Section 5 of the book chapter [START_REF] Chertkov | Ensemble control of cycling energy loads: Markov Decision Approach[END_REF] a similar optimal control formulation is proposed:

min p D(p∥p 0 ) subject to E p Y(X k ) = r k , 1 ≤ k ≤ K (17)
Relative entropy is a useful measure of cost of deviation from nominal behavior because the optimizer has a simple form: a "tilting" (or "twisting") of the nominal model [START_REF] Todorov | Linearly-solvable Markov decision problems[END_REF], [START_REF] Meyn | Ancillary service to the grid from deferrable loads: The case for intelligent pool pumps in Florida[END_REF]. This is motivation for the use of relative entropy in this prior work, which leads to the solution to [START_REF] Almassalkhi | Packetized energy management: asynchronous and anonymous coordination of thermostatically controlled loads[END_REF]:

p * (⃗ x) = p 0 (⃗ x) exp K k=1 β k Y(x k ) -Λ(β)
in which β ∈ R K are Lagrange multipliers corresponding to the average power constraints, and Λ(β) a normalizing constant.

The optimization criterion ( 17) is a form of dead-beat control, which might cause concern: is the optimization problem feasible? are there stability issues, as is well known for dead-beat control of linear systems?

For the KLQ formulation described in this paper, the tracking constraint in ( 17) is replaced by a quadratic loss function. As κ → ∞ we recover the solution to the deadbeat control problem (however, our final result is different than the solution in [START_REF] Chertkov | Ensemble control of cycling energy loads: Markov Decision Approach[END_REF] wherein ν 0 is not constrained).

The convex program formulation ( 16) has many advantages. First, ( 16) is always feasible, while feasibility of ( 17) requires conditions on p 0 and r. Theorem 3.1 requires no assumptions on the model or reference signal. Second is the value of flexibility in choice of κ, so that we can learn what is an "expensive" reference signal. It is anticipated that the penalty parameter κ can be used to make tradeoffs between tracking performance and robustness to modeling error: robustness and sensitivity analysis will be a topic of future research.

C. KLQ and IPD

The finite-horizon version of ( 6) is also considered in [START_REF] Meyn | Ancillary service to the grid from deferrable loads: The case for intelligent pool pumps in Florida[END_REF], [START_REF] Bušić | Distributed randomized control for demand dispatch[END_REF], similar to the KLQ formulation:

p ζ := arg max p ζE p K k=1 Y(x k ) -D(p∥p 0 ) . (18)
This and ( 6) are versions of the Individual Perspective Design (IPD) [START_REF] Bušić | Distributed randomized control for demand dispatch[END_REF].

The IPD design [START_REF] Benenati | A tractable formulation for multi-period linearized optimal power flow in presence of thermostatically controlled loads[END_REF] has the following alternative interpretation. For a scalar r 0 ∈ R, consider the constrained optimization problem

max p -D(p∥p 0 ) subject to E p K k=1 Y(x k ) = Kr 0 , 1 ≤ k ≤ K (19)
The dual function φ * : R → R is defined by

φ * (λ) = max p λE p K k=1 Y(x k ) -D(p∥p 0 ) -λKr 0
where λ ∈ R is a Lagrange multiplier. It is evident that the optimizer p * λ is an IPD solution for each λ. Consequently, for each ζ, the IPD solution (18) also solves [START_REF] Meyn | Ancillary service to the grid using intelligent deferrable loads[END_REF] for some scalar r 0 (ζ).

Note however that the IPD approach is designed to construct a family of transition matrices {P ζ }. This is then used in a pure feedback control design in which {ζ k } is broadcast to the loads from the BA. The present work is focused on feed-forward control.

III. DUALITY

Structure for the solution of ( 16) will be obtained by consideration of a dual, in which λ ∈ R K and g ∈ R K×|S| denote the vectors of Lagrange multipliers for the first and second set of constraints, respectively. The Lagrangian is thus

L(ν, γ, λ, g) = K k=1 D(ν k , ν 0 k ) + κ 2 K k=1 γ 2 k + K k=1 λ k b k + K k=1 s ′ g k (s ′ )c k (20)
where b k and c k refer to the left-hand-side of equations ( 16b) and (16c), respectively. The dual function is the minimum:

φ * (λ, g) := min ν,γ L(ν, γ, λ, g)
and the dual of the optimization problem ( 16) is defined as the maximum of the dual function φ * over λ and g. We will see that there is no duality gap, so that for a quadruple (ν * , γ * , λ * , g * ),

J * (ν 0 0 ) = L(ν * , γ * , λ * , g * ) = φ * (λ * , g * ) .
In the following subsections we obtain a representation of the dual function that is suitable for optimization, and in doing so we obtain a representation for the optimal policy. Properties of the dual function are contained in Theorem 3.1 and Prop. 3.2 that follow. The statement of these results requires additional notation: define a function

T λ k : R |S| → R |S| , for f : S → R, λ ∈ R K , and s ∈ R |S| via T λ k (f ; s) = log u ϕ 0 k (u | s) exp s ′ T k (x, s ′ )f (s ′ ) + λ k Y(s, u)
The maximum of the dual function over g is denoted

φ * (λ) := max g ϕ * (λ, g) = φ * (λ, g λ )
where g λ is a maximizer:

g λ ∈ arg max g ϕ * (λ, g)
We will show that the sequence of functions g λ is given by the recursion

g λ k = T λ k (g λ k+1 ) , 1 ≤ k ≤ K
, where g λ K+1 ≡ 0 (21) and denote:

G λ k (x k-1 ) = s T u k-1 (s k-1 , s)g λ k (s) (22)
Theorem 3.1: There exists a maximizer {λ * k , g * k : 1 ≤ k ≤ K} for φ * , and there is no duality gap:

φ * (λ * , g * ) = J * (ν 0 0 )
The optimal policy is obtained from {g * k } via:

ϕ * k (u | s) =ϕ 0 k (u | s) exp s ′ T k (x, s ′ )g * k+1 (s ′ ) + λ * k Y(s, u) -g * k (s) where g * k (s) =T λ k (g * k+1 ; s) and g * K+1 ≡ 0 (23)
⊓ ⊔ Proposition 3.2: The following hold for the dual of ( 16): for each λ ∈ R K , (i) A maximizer g λ is given by [START_REF] Bušić | Distributed randomized control for demand dispatch[END_REF] (ii) The maximum of the dual function over g is the concave function

φ * (λ) = λ T r - 1 2κ ∥λ∥ 2 -⟨ ν 0 , G λ 1 ⟩ (24)
(iii) The function ( 24) is continuously differentiable, and

∂ ∂λ k φ * (λ) = r k - 1 κ λ k -⟨ ν λ k , Y ⟩ (25)
where {ν λ k } is the sequence of marginals obtained from the randomized policy defined in [START_REF] Chertkov | Ensemble control of cycling energy loads: Markov Decision Approach[END_REF], substituting {g * k } by {g λ k } defined in (i).

⊓ ⊔

To conclude this section, we provide representations of the log-likelihood ratio L(⃗ x), relative entropy D(p λ ∥p 0 ), and primal objective function,

J(p λ , ν 0 0 ) := D(p λ ∥p 0 ) + κ 2 K k=1 ⟨ ν λ k , Y ⟩ -r k 2 ( 26
)
where p λ is the pmf obtained from the randomized policy defined in [START_REF] Chertkov | Ensemble control of cycling energy loads: Markov Decision Approach[END_REF], substituting {g * k } by {g λ k } defined in Prop. 3.2, part (i).

Corollary 3.3: The following hold for all {λ k , g λ k : 1 ≤ k ≤ K}:

(i) The log-likelihood ratio can be expressed:

L(⃗ x) = K k=1 {∆ k (x k-1 , s k ) + λ k Y(x k )} -G λ 1 (x 0 ) (27)
where for each k (recalling

x k = (s k , u k )), ∆ k (x k-1 , s k ) = G λ k (x k-1 ) -g λ k (s k ) (28)
(ii) The relative entropy is given by

D(p λ ∥p 0 ) = K k=1 λ k ⟨ ν λ k , Y ⟩ -⟨ ν 0 , G λ 1 ⟩ (29)
(iii) The value of the primal is

J(p λ , ν 0 0 ) = -⟨ ν 0 , G λ 1 ⟩ + K k=1 λ k ⟨ ν λ k , Y ⟩ + κ 2 ⟨ ν λ k , Y ⟩ -r k 2 ⊓ ⊔
The stochastic process {∆ k (X k-1 , S k )} is a martingale difference sequence; it vanishes when nature is deterministic, reducing to the solution obtained in [START_REF] Cammardella | Kullback-Leibler-Quadratic optimal control of flexible power demand[END_REF], [START_REF] Cammardella | Simultaneous allocation and control of distributed energy resources via Kullback-Leibler-Quadratic optimal control[END_REF].

IV. NUMERICAL EXPERIMENTS

Numerical experiments are conducted in the context of Demand Dispatch. The goal here is to modify the behavior of flexible loads such that their aggregate power consumption tracks a reference signal {r k } that is broadcast by a BA. Previous work has demonstrated the potential for pool pumps [START_REF] Meyn | Ancillary service to the grid from deferrable loads: The case for intelligent pool pumps in Florida[END_REF], HVAC systems [START_REF] Hao | How demand response from commercial buildings will provide the regulation needs of the grid[END_REF], water heaters [START_REF] Mathias | Demand dispatch with heterogeneous intelligent loads[END_REF], and refrigerators [START_REF] Mathias | Smart fridge / dumb grid? Demand Dispatch for the power grid of 2020[END_REF] to provide grid services. The following numerical experiments demonstrate distributed control of a collection of homogeneous residential refrigerators.

A. Algorithms

We have found in examples that using gradient ascent on the dual function curve may be slow to converge, likely due to a large "overshoot" when applying standard first-order methods [START_REF] Cammardella | Kullback-Leibler-Quadratic optimal control of flexible power demand[END_REF]. In the numerical results that follow we opt for proximal gradient methods [START_REF] Parikh | Proximal Algorithms, ser. Foundations and Trends in Optimization[END_REF]. Monte Carlo methods have also been used to estimate λ * . This is motivated by the representation of the gradient in terms of the first-order statistics of the random variable {Y(X)} when X ∼ ν λ k :

E k [Y(X)] = x k ∈X ν λ k (x k )Y(x k ) = ⟨ ν λ k , Y ⟩ (30)
Lemma 4.1 follows from (25) combined with (30): Lemma 4.1:

For any λ ∈ R K and 1 ≤ k ≤ K, ∂ ∂λ k φ * (λ) = r k - 1 κ λ k -E k [Y(X)] , X ∼ ν λ k . (31)
⊓ ⊔

B. Designing the nominal model

The nominal model is carefully designed to ensure the uncontrolled dynamics remain unchanged. The state-input space is the cartesian product of two state-input spaces: X = X n × X u , where X n contains the uncontrollable components and X u contains the the controllable components. The nominal state transition matrices (9) that define the nominal model (8) are products of T k , which represents the uncontrollable components, and ϕ 0 , which represents the nominal control policy. In order to satisfy this decomposition, the state at time k for a refrigerator is defined to be

S k = (θ k , θ k-1 , U k-1 ) (32)
where θ k ∈ R denotes the temperature inside the refrigerator and U k ∈ {0, 1} is the power mode. The uncontrollable component T k can not be modified and is derived from the linear model ( 5), with the addition of disturbances to represent randomness from nature. The nominal policy ϕ 0 k is designed to approximate the deterministic control that is standard for a refrigerator [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF].

24 hours . . . . . . . . . GW Truth Pred. N ν k ,Y

C. Designing the Experiment

The output process as a function of the state is defined to be the power consumption of a single refrigerator

Y(x) = uϱ , x = (s, u) ∈ X .
where ϱ is its rated power consumption. In the plots shown below, the expected power consumption ⟨ν k , Y⟩ and the reference signal r k are multiplied by the number of refrigerators N to show the aggregate response. The reference signal is obtained from the solution to a separate optimization problem [START_REF] Cammardella | Balancing California's grid without batteries[END_REF] which was designed to calculate optimal power trajectories for distinct classes of flexible loads. These power trajectories are optimal with respect to minimizing peaks and ramps in net load while satisfying quality of service constraints such as maintaining a given temperature range inside a refrigerator.

To simulate prediction error, the reference signal is created by adding random noise to the optimal power trajectory. Data from the California Independent System Operator (CAISO) show that hour-ahead predictions are fairly accurate, which inspires our use of model predictive control (MPC): fix two time periods T 0 and T , where T 0 ≪ T ; we choose T 0 = 1 hour and T = 12 hours. At the initial time t 0 , the marginal pmf ν 0 is estimated, the reference signal is predicted over the time window [t 0 , t 0 + T ], where its value at time t 0 is observed and has zero prediction error, and a solution is computed over the time window [t 0 , t 0 + T ]. Then, the solution is implemented, but restricted to the smaller time window [t 0 , t 0 + T 0 ]. At time t 0 + T 0 , the marginal pmf is estimated, the reference signal is predicted over the time window [t 0 + T 0 , t 0 + T 0 + T ], where its value at time t 0 + T 0 is observed and has zero prediction error, a solution is computed over [t 0 + T 0 , t 0 + T 0 + T ], and the process continues.

Truth

Pred.

N ν k ,Y . . . . . . . . . GW 24 hours

D. Results

The first experiment, shown in Fig. 1, depicts an ideal case, where the reference signal is feasible, with respect to satisfying strict temperature bounds for each refrigerator. The variance of the prediction error is chosen to match the variance of the prediction error in hour-ahead forecasts observed by CAISO on February 25, 2021 [START_REF]CA 95763-9014[END_REF]. The reference signal is an optimal power trajectory for 20 million refrigerators on a typical sunny day in California. Notice how the refrigerators are asked to consume more during the day when solar output is high and to ramp down during the evening peak. Tracking is nearly perfect in this ideal case.

The second experiment, shown in Fig. 2, depicts a case where the reference signal is made infeasible by reducing the number of participating loads to 7 million. Satisfying the energy requirements of this signal would require refrigerators to violate their temperature bounds. This is not possible since quality of service is guaranteed with our distributed control control architecture, where switching decisions are still made at the load level [START_REF] Cammardella | Kullback-Leibler-Quadratic optimal control of flexible power demand[END_REF]. Notice how the collective power consumption is gracefully truncated at the peak and valley of the reference signal while tracking remains nearly perfect at all other times.

V. CONCLUSIONS

A Lagrangian decomposition separates the finite-horizon optimal control problem into K separate convex programs, one for each time step. By applying a few well-known concepts from information theory, the optimal policy at each time step is found to be an exponential tilting of the nominal policy. The main novelty is to allow for more general Markovian models. Numerical experiments demonstrate the usefulness of this distributed control technique in a power system setting: a collection of flexible loads can be controlled such that their aggregate power consumption tracks a reference signal.

Plans for future research include: (i) Evaluate robustness and sensitivity (ii) Investigate techniques to reduce computation (iii) Extend this stochastic KLQ formulation into a continuous-time setting and search for more effective relaxation techniques; initial work may be found in [START_REF] Bušić | Distributed control of thermostatically controlled loads: Kullback-Leibler optimal control in continuous time[END_REF].

(iv) Consider other cost functions, e.g., the Wasserstein distance.

(v) Investigate the relationship between optimality and coupling of the pmfs, and the implications to control design.

(vi) Careful design of a terminal cost function may result in better performance for smaller time horizons [START_REF] Chen | Value iteration and optimization of multiclass queueing networks[END_REF].

APPENDIX

This appendix contains abbreviated proofs of some of the results of this paper. Complete proofs of every result can be found in [START_REF] Cammardella | Kullback-Leiblerquadratic optimal control[END_REF].

A. Proof of Prop. 2.1 Proof: Clearly, ( 16) is convex in γ k , a quadratic function of ν k . Convexity in ν k is established by constructing a sub-gradient for D(ν k , ν 0 k ): the function

g µ,k (s, u) = log µ(s, u) ν 0 k (s, u) -log μ(s) ν0 k (s) , (33)
satisfies the sub-gradient property

D(ν k , ν 0 k ) ≥ D(µ, ν 0 k ) + ⟨ ν k -µ, g µ,k ⟩ (34)
where μ(s) = u ′ µ(s, u ′ ) and νk (s) = u ′ ν k (s, u ′ ). Hence, ( 16) is jointly convex in ν k and γ k since it is a sum of convex functions.

The equivalence claim is proven by multiplying both sides of (16c) by ϕ k (u ′ | s ′ ), yielding

P{X k = x ′ } = x P{X k-1 = x}P{X k = x ′ | X k-1 = x}
which is identical to [START_REF] Chen | Distributed control design for balancing the grid using flexible loads[END_REF]. The proof of the implication ( 11) =⇒ (16c) is similar. A complete proof can be found in [START_REF] Cammardella | Kullback-Leiblerquadratic optimal control[END_REF].

B. Relative entropy and duality

The proofs of Theorem 3.1 and Prop. 3.2 make use of the following four lemmas. The first is based on a well known result regarding relative entropy. For any function h :

X K+1 → R denote Ł 0 (h) := sup p ⟨p, h⟩ -D(p∥p 0 ) (35)
Lemma 1.1 (Convex dual of relative entropy): For each p 0 and function h : X K+1 → R, the (possibly infinite) value of (35) coincides with the log moment generating function: Ł 0 (h) = log⟨ p 0 , e h ⟩ Moreover, provided Ł 0 (h) < ∞, the supremum in ( 35) is uniquely attained with p * = p 0 exp(h -Ł 0 (h)). That is, the log-likelihood L * = log(dp * /dp 0 ) is given by

L * (⃗ x) = h(⃗ x) -Λ 0 (h) ⊓ ⊔ Lemma 1.2:
The dual function can be expressed

φ * (λ, g) =λ T r - 1 2κ ∥λ∥ 2 -⟨ ν 0 , G λ 1 ⟩ + K k=1 min s g k (s) -T λ k (g k+1 ; s) (36)
where g K+1 ≡ 0.

⊓ ⊔

Proof: First, substitute ν k (s, u) = νk (s)ϕ k (u | s) in the Lagrangian [START_REF] Todorov | Linearly-solvable Markov decision problems[END_REF] so that its minimization requires obtaining each of the K minimizers {ν λ,g k : ν λ,g k (s, u) = νλ,g k (s)ϕ λ,g k (u | s)}. It follows from Lemma 1.1 that the minimizers are given by

ϕ λ,g k (u | s) =ϕ 0 k (u | s) exp s ′ T k (x, s ′ )g k+1 (s ′ ) + λ k Y(s, u) -Λ k (s) (37)
with Λ k (s) = T λ k (g k+1 ; s). Lemma 1.1 also gives the value:

arg min ϕ k L = -T λ k (g k+1 ; s) (38)
resulting in min ν L(ν, γ, λ, g)

= K k=1 κ 2 γ 2 k + λ k γ k + λ k r k - s,u ν 0 0 (s, u)G λ 1 (s, u) + K k=1 min νk ⟨ νk , g k -T λ k (g k+1 ) ⟩ (39)
Next, observe that the minimizer νλ,g k is obtained when the support of each νk satisfies supp νk (s) ⊆ arg min

s g k (s) -T λ k (g k+1 ; s) so that min s g k (s) -T λ k (g k+1 ; s) = ⟨ νλ,g k , g k -T λ k (g k+1 ) ⟩ Also, the minimizer γ λ k is γ λ k = - 1 κ λ k (40)
Substituting the minimizers {ν λ,g k , γ λ k } into (39) results in (36).

C. Duality

Lemma 1.3: The maximum of the dual function over g is

φ * (λ) := max g φ * (λ, g) = λ T r - 1 2κ ∥λ∥ 2 -⟨ ν 0 , G λ 1 ⟩ (41) with G λ 1 (s, u) = s ′ T k (x, s ′ )g λ 1 (s ′ ).
A maximizer g λ is given by the recursive formula:

g λ k = T λ k (g λ k+1 ) , 1 ≤ k ≤ K , where g λ K+1 ≡ 0 , (42) ⊓ ⊔
Proof: Adding a constant to any of the (g 1 , g 2 , . . . , g K ) does not change the value of L or φ * (this follows from ( 21)), so without loss of generality we assume, for each k,

min s g k (s) -T λ k (g k+1 ; s) = 0 (43)
and consequently

g k ≥ T λ k (g k+1 ) for each k . (44)
Thus, in view of (36), φ * (λ) =

λ T r - 1 2κ ∥λ∥ 2 -min g1 s,u ν 0 0 (s, u) s ′ T k (x, s ′ )g 1 (s ′ ) , (45)
where the minimum is subject to the constraint (44). Next, observe that T λ k is a monotone operator, so that for each k ≤ K,

g k ≥ T λ k • T λ k+1 • • • • • T λ K (g K+1 ) . = g λ k ,
where g K+1 ≡ 0

Based on the expression (45), we now show that the maximum arg max g ϕ * (λ, g) is obtained by choosing each g k to reach this lower bound, giving (42). Indeed, g λ 1 achieves the minimum in (45), since g λ 1 ≤ g 1 for any g 1 for which (44) holds. This result along with (43) yields (41).

Lemma 1.4: The maximizers {g λ k } have at most linear growth in ∥λ∥:

|g λ k (s)| ≤ ∥Y∥ ∞ K i=k |λ i | ≤ √ K∥Y∥ ∞ ∥λ∥ (46)

⊓ ⊔

Proof: The proof is by induction, starting with the base case k = K. Then, we assume the hypothesis is true for k ≤ K and show it holds for k -1.

Next, we present a proof of Theorem 3.1.

Proof: We prove the existence of a maximizer λ * by showing that ϕ * (λ) is an anti-coercive function, i.e., ϕ * (λ) → -∞ as ∥λ∥ → ∞. By Lemma 1.4,

φ * (λ) = λ T r - 1 2κ ∥λ∥ 2 - s,u ν 0 0 (s, u) s ′ T k (x, s ′ )g λ 1 (s ′ ) ≤ ∥λ∥∥r∥ - 1 2κ ∥λ∥ 2 + max s ′ |g λ 1 (s ′ )| ≤ ∥λ∥∥r∥ - 1 2κ ∥λ∥ 2 + √ K∥Y∥ ∞ ∥λ∥
Since ϕ * (λ) is upper-bounded by an anti-coercive function, ϕ * (λ) itself is an anti-coercive function. Thus a maximizer λ * exists, and (λ * , g * ) = (λ * , g λ * ) by ( 42).

The primal is a convex program, as established in Prop. 2.1. To show that there is no duality gap it is sufficient that Slater's condition holds [START_REF] Boyd | Convex Optimization[END_REF]Section 5.3.2]. This condition holds: the relative interior of the constraint-set for the primal is non-empty since it contains {ν 0 k }. Optimality of ( 23) is established by substituting g * k+1 into (37) and by making the substitution g * k = T λ k (g * k+1 ) implied by (42). Next, we present a proof of Prop. 3.2.

Proof: This proof has three parts:

(i) ( 21) is proven by Lemma 1.3.

(ii) ( 24) is proven by Lemma 1.3.

(iii) The representation of the derivative in part (iii) is standard (e.g., Section 5.6 of [START_REF] Boyd | Convex Optimization[END_REF]), but we provide the proof for completeness. The representation [START_REF] Todorov | Linearly-solvable Markov decision problems[END_REF] implies that φ * is concave in (λ, g), since it is the infimum of linear functions. This representation also gives a formula for a derivative:

∂ ∂λ k φ * (λ, g) = r k - 1 κ λ k -⟨ ν λ,g k , Y ⟩
where ν λ,g k is any optimizer in (39). Using φ * (λ) = φ * (λ, g λ ) then gives

∂ ∂λ k φ * (λ) = r k - 1 κ λ k -⟨ ν λ k , Y ⟩ + ∂ ∂g φ * (λ, g λ ) • ∂ ∂λ k g λ
The first order condition for optimality gives ∂ ∂g φ * (λ, g λ ) = 0, which completes the proof of the representation. It is evident that φ * is continuously differentiable since ν λ k is continuously differentiable for each k by construction.

Next, we present a proof of Corollary 3.3.

Proof: This proof has three parts:

(i) Application of ( 15) and ( 23) results in the loglikelihood ratio:

L(⃗ x) = K k=1 s T k (x k , s)g λ k+1 (s) + λ k Y(x k ) -g λ k (s k ) = K k=1 G λ k+1 (x k ) + λ k Y(x k ) -g λ k (s k )
where the second identity follows from the definition [START_REF] Cammardella | Kullback-Leiblerquadratic optimal control[END_REF].

We have from the definitions, G λ K+1 ≡ 0, which results in

L(⃗ x) = -G λ 1 (x 0 ) + K k=1 G λ k (x k-1 ) + λ k Y(x k ) -g λ k (s k )
This combined with (28) yields [START_REF] Cammardella | Balancing California's grid without batteries[END_REF].

(ii) Applying the definition of relative entropy as the mean log-likelihood, and noticing that

E p λ ∆ k (X k-1 , S k ) = 0 for 1 ≤ k ≤ K, results in ⃗ x p λ (⃗ x)L(⃗ x) = ⃗ x p λ (⃗ x) K k=1 λ k Y(x k ) -G λ 1 (x 0 ) = K k=1 x k xi,i̸ =k p λ (⃗ x)λ k Y(x k ) - x0 xi,i̸ =0 p λ (⃗ x)G λ 1 (x 0 ) = K k=1 λ k ⟨ ν λ k , Y ⟩ -⟨ ν 0 , G λ 1 ⟩
(iii) Substitution of ( 29) into (26) results in (3).

Fig. 1 .

 1 Fig. 1. 20 million refrigerators tracking a feasible signal

Fig. 2 .

 2 Fig. 2. 7 million refrigerators tracking an infeasible signal

Department of Electrical and Computer Engineering, University of Florida, Gainesville, FL 32611. SM's research also supported through an INRIA International Chair.

Inria and DI ENS, École Normale Supérieure, CNRS, PSL Research University, Paris, France


PhD Olivier M Girard
email: olivier.girard@univ-amu.fr

Virginie Callot

Valentin H Prevost

Benjamin Robert

Manuel Taso

Guilherme Ribeiro

Gopal Varma

Novena Rangwala

David C Alsop

Guillaume Duhamel

Magnetization transfer from inhomogeneously broadened lines (ihMT): Improved imaging strategy for spinal cord applications

Keywords: Spinal cord, inhomogeneous magnetization transfer, CSF pulsatility, ECG synchronization, motion correction, myelin, white matter

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Inhomogeneous magnetization transfer (ihMT) shows great promise for central nervous system (CNS) imaging thanks to its apparent high myelin specificity compared to conventional magnetization transfer (MT) [START_REF] Varma | Magnetization Transfer from Inhomogeneously Broadened Lines: A Potential Marker for Myelin[END_REF][START_REF] Girard | Magnetization transfer from inhomogeneously broadened lines (ihMT): Experimental optimization of saturation parameters for human brain imaging at 1.5 Tesla[END_REF]. Although important milestones, such as the development of a theoretical model [START_REF] Varma | Interpretation of magnetization transfer from inhomogeneously broadened lines (ihMT) in tissues as a dipolar order effect within motion restricted molecules[END_REF], have been achieved, thorough understanding of the ihMT effect is still in an early phase. The high specificity of ihMT toward myelin is presumably attributed to its unique architecture, consisting of multiple sheaths of lipid bilayers, and specific molecular composition, notably high lipid and relatively low protein content [START_REF] Baumann | Biology of oligodendrocyte and myelin in the mammalian central nervous system[END_REF][START_REF] Horch | Origins of the ultrashort-T21H NMR signals in myelinated nerve: A direct measure of myelin content?[END_REF]. Whereas proton-proton dipolar interaction usually results in homogeneous broadening [START_REF] Portis | Electronic Structure of $F$ Centers: Saturation of the Electron Spin Resonance[END_REF] of the NMR spectrum for large macromolecules, the rapid axial rotation of the lipid chains composing the myelin results in inhomogeneous broadening for biological membranes [START_REF] Davis | High resolution 1H nuclear magnetic resonance of a transmembrane peptide[END_REF][START_REF] Chen | Highresolution MAS NMR spectroscopy detection of the spin magnetization exchange by crossrelaxation and chemical exchange in intact cell lines and human tissue specimens[END_REF]. In addition multiple mechanisms such as spin diffusion, chemical exchange and molecular motion, which usually tend to homogenize broad semi-solid NMR lines in tissue, seem highly constrained or inefficient in the case of myelin sheath environment. Molecular motions within lipid bilayers are restricted as compared to free lipids, water exchange from the myelin water to the bulk water is relatively slow [START_REF] Vavasour | Different magnetization transfer effects exhibited by the short and long T(2) components in human brain[END_REF], and spin diffusion within lipid membranes is a slow process [START_REF] Huster | Membrane Protein Topology Probed by 1 H Spin Diffusion from Lipids Using Solid-State NMR Spectroscopy[END_REF].

Altogether, these evidences strongly support the myelin spectrum being inhomogeneous on the time scale of several milliseconds, in compliance with the recent measure and sensitivity demonstration of the ihMT effect in vivo.

So far the ihMT technique has been successfully applied for brain applications at various field strengths, including clinical 1.5T [START_REF] Girard | Magnetization transfer from inhomogeneously broadened lines (ihMT): Experimental optimization of saturation parameters for human brain imaging at 1.5 Tesla[END_REF] and 3T (1) as well as preclinical 11.75T [START_REF] Prevost | Magnetization Transfer from Inhomogeneously Broadened Lines (ihMT): Qualitative Evaluation of ihMT Specificity toward Myelinated Structures[END_REF]. Initial reports concerned both methodological developments and applications to demyelinating pathologies such as multiple sclerosis [START_REF] Varma | Analysis of Inhomogeneous Magnetization Transfer Applied In Patients with Multiple Sclerosis[END_REF][START_REF] Duhamel | Magnetization Transfer from Inhomogeneously Broadened Lines (ihMT): Application on Multiple Sclerosis (MS)[END_REF]. In addition, recent preliminary studies have also described the feasibility of ihMT for spinal cord (SC) investigations [START_REF] Rangwala | Quantification of myelin in the cervical spinal cord using inhomogeneous magnetization transfer imaging[END_REF][START_REF] Girard | Magnetization Transfer from Inhomogeneously Broadened Lines (ihMT): Improved SC Imaging using ECG Synchronization[END_REF], which is of the utmost importance since the technique may open new perspectives to better describe the SC healthy structure and its pathological alterations. Challenges for SC ihMT studies include a relatively coarse resolution compared to the small size of the cord, the presence of susceptibility effects from surrounding tissues (e.g. bones), and most importantly sensitivity to the pulsatility of the cerebrospinal fluid (CSF), which induces motion in and around the cord [START_REF] Friese | The influence of pulse and respiration on spinal cerebrospinal fluid pulsation[END_REF]. Indeed, in comparison with conventional MT, such motion may have more dramatic impact on ihMT because it has lower sensitivity (ihMT signal on the order of 10% of the unsaturated signal) and requires data subtraction from multiple images potentially acquired at different phases of motion. Hence, as has been observed in DTI studies [START_REF] Summers | A preliminary study of the effects of trigger timing on diffusion tensor imaging of the human spinal cord[END_REF][START_REF] Kharbanda | Effects of cord motion on diffusion imaging of the spinal cord[END_REF], pulsatile motion could have critical effects on SC image quality and potential contamination of SC tissue signal may lead to inaccurate ihMT quantification.

The current report expands on earlier work [START_REF] Girard | Magnetization Transfer from Inhomogeneously Broadened Lines (ihMT): Improved SC Imaging using ECG Synchronization[END_REF] and presents a complete strategy for cervical SC ihMT imaging free of CSF pulsatility artifacts. The proposed method combines electro-cardiogram (ECG) synchronization and retrospective data pairing. This strategy aims at correcting for signal variations of long T1-components, such as CSF, that occur when a variable recovery delay is introduced in the sequence as a consequence of variable cardiac cycle. In comparison with other strategies correcting for such T1-recovery effects [START_REF] Guimaraes | Imaging Subcortical Auditory Activity in Humans[END_REF][START_REF] Shin | Correction for the T1 effect incorporating flip angle estimated by Kalman filter in cardiac-gated functional MRI[END_REF], the proposed method is model-free as it simply matches images acquired with a similar recovery period. A semi-automatic algorithm for data correction is presented along with accurate ihMT ratio quantification of the cord.

Methods

All experiments were performed at 1.5T (Siemens Avanto, Erlangen, Germany) on healthy volunteers, in compliance with guidelines from our institutional committee on clinical investigations.

ihMT magnetization preparation

A pulsed ihMT preparation module was used, as described previously [START_REF] Girard | Magnetization transfer from inhomogeneously broadened lines (ihMT): Experimental optimization of saturation parameters for human brain imaging at 1.5 Tesla[END_REF]. The ihMT contrast generation relies on the comparison of two different MT weighted images (1): one using a single-frequency offset (+Δf) leading to the saturated image MT + and the other using a dual-frequency offset and identical total energy (+Δf and -Δf applied sequentially) leading to the saturated image MT +-. In order to compensate for potential MT asymmetry effects [START_REF] Varma | Magnetization Transfer from Inhomogeneously Broadened Lines: A Potential Marker for Myelin[END_REF][START_REF] Girard | Magnetization Transfer from Inhomogeneously Broadened Lines (ihMT): Improved SC Imaging using ECG Synchronization[END_REF], two additional single (MT -) and dual-frequency (MT -+ ) acquisitions are measured and the final ihMT image is generated by the following combination:

Eq. 1 Note that both dual-frequency saturation weighted images, MT +-and MT -+ , are essentially similar. The ihMT ratio is defined by ihMTR = ihMT/S0, where S0 is the signal measured with RF saturation power set to zero.

+ - - + - + - - + = MT MT MT MT ihMT
The following parameters were used, derived from brain imaging optimization at 1.5T

(2): frequency offsets ±Δf=±7khz, pulsewidth pw=500μs and interpulse repetition time Δt=1ms. The total saturation time (t), optimized at t=700ms for brain studies, had to be adapted to the ECG synchronization required in SC imaging in order to ensure acquisition of the image during the quiescent phase of the cervical spinal cord. This phase starts on average about 300ms after the R-wave and lasts around 50% of the whole cardiac cycle [START_REF] Summers | A preliminary study of the effects of trigger timing on diffusion tensor imaging of the human spinal cord[END_REF] (typically in the range of 800-1100ms). Hence the MT saturation, triggered by the R-wave, was reduced to t=500ms, i.e. leaving about 200-350ms for image readout during the quiescent phase. The time-integrated squared B1 (i.e. the average RF energy absorbed during the saturation period τ, noted ETR) was adjusted so that the RF energy of the saturation pulsed train corresponded to ETR = 17.2μT².s. This value is lower than the optimized value used for brain (ETR ≥ 30μT².s), hence slightly reducing the sensitivity of the technique (see Fig. 7 of reference ( 2)). However this reduced ETR value was chosen to allow scanning under regulatory Specific Absorption Rate (SAR) limitations for all the investigated nominal repetition times (TR) down to 2500ms.

Imaging

A product 12-channel head / 4-channel neck combined coil was used for the entire study.

The ihMT preparation was combined with a single-shot 2D HASTE readout module, to provide low sensitivity to susceptibility effects, high SNR and rapid acquisition time, desirable features for SC imaging optimization. Such readout, however, is only suitable for single-slice imaging because of rather long echo-train duration. The sagittal plane was chosen to acquire simultaneously signals from cervical spinal cord levels C2-C4 and from a reference brain area (i.e. the corpus callosum -CC). The following parameters were used: 25cm Field-Of-View, 192×192 matrix size, 1.3mm in-plane resolution, 789Hz/pixel bandwidth, 120° refocusing flip angles, no parallel acceleration. Partial volume effects due to WM/GM structure distribution within the SC were assessed using two different slice thicknesses: 5mm and 8 mm. In addition, potential blurring effects related to the long HASTE readout echo-train were assessed by applying Phase-Encoding (PE) in the Antero-Posterior (A-P) and the Superior-Inferior (S-I) directions. For the S-I direction, 50% phase oversampling was prescribed to prevent aliasing.

The shortest possible echo times (TEs) were used, corresponding to 31ms in the S-I direction (for a total echo-train length of ~550ms), and to 20ms in the A-P direction (echo-train length of ~365ms). 50% phase oversampling was also tested in the A-P direction to provide a fair comparison (similar TE settings as in the S-I direction).

Image acquisition was synchronized to the ECG and four nominal TR settings were studied: TR=2500, 3500, 4500 and 10000 ms. The latter configuration corresponded to a reference acquisition, for which full recovery of the longitudinal magnetization is expected for all human tissues including CSF. For each scan, ECG physiological data and trigger instances were recorded in a log file throughout the scan in order to allow further retrospective data correction.

In order to increase the signal-to-noise ratio (SNR) of the ihMT images, 30 signal averages (NEX) of each of the four kinds of MT weighted images (see Eq. 1) were acquired for all conditions (except for the reference acquisition -i.e. TR=10000ms, for which 20 NEX were acquired to reduce scan time), and 3 NEX for S0. Overall, the total number of images per studied condition was 123 images (83 for TR=10000ms), leading to ~ 6.5 min, 8.5 min, 10.5 min and 14.5 min total acquisition time for TR of 2500, 3500, 4500 and 10000 ms, respectively. Each aforementioned imaging condition was tested on three healthy volunteers (males, 33 ± 6 y.o.).

Post-Processing

Because of the small size of the SC, only ~ 4-6 pixels covered its antero-posterior extension for the 1.3mm in-plane resolution used in this study, and hence, significant Gibbs' ringing occurred. This well-known MR artifact can have dramatic impact on SC image quality because it adds a non-local signal component corresponding to signal secondary-lobes coming from neighboring tissues. These lobes typically arise 1-2 voxel away from edges, and hence tend to add CSF signal onto the SC which should be avoided in the current study which aims to measure unbiased SC ihMT metrics.

Windowing the k-space data helps to strongly attenuate these confounding signal oscillations and also improves SNR at the same time, although at the cost of a broadening or blurring of the main lobe, causing a loss in spatial resolution. Here a Hamming-shape window was used, offering a fair compromise between secondary lobe attenuation and blurring [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF].

The full width at half maximum corresponds to 1.81× the nominal resolution as compared to 1.21 for unfiltered data, resulting in a reduction of the effective resolution by a factor of 1.5.

Additionally, a zero-filling factor of 2 was used for selection of regions-of-interest (ROIs) and subsequent quantitative measurements. The effect of k-space windowing is illustrated in Fig. 1.

The multiple averages of the different instances of the S0 and MT weighted images were all registered together using the SPM8 realign routine [START_REF] Penny | Statistical Parametric Mapping: The Analysis of Functional Brain Images: The Analysis of Functional Brain Images[END_REF] prior to image combination.

ECG synchronization and retrospective data pairing

Image acquisition synchronized with ECG led to effective TR values that depend on the cardiac cycle. Here, the trigger initiated the ihMT preparation module and the variable delay was placed after the readout, hence lengthening the signal recovery period. Any TR variation from one image acquisition to the other directly translates into different T1-weighting of long T1 components (e.g. CSF) in individual MT images (except for the fully-relaxed acquisition).

The result of combining all the MT images regardless of their actual TR values would induce errors in the ihMT signal because it relies on image subtraction (see Eq. 1). In particular, the CSF compartment may give rise to an artefactual non-zero ihMT signal, which could affect SC tissue conspicuity and ihMT metrics because of partial volume effects.

A custom retrospective pairing algorithm was developed in MATLAB (R2010b, The MathWorks Inc., Natick, MA, USA) and tested to correct for such residual errors in the ECGtriggered ihMT images. The proposed algorithm basically works by extracting the actual TR values from the recorded physiological data, and then applying a lookup loop to generate pairwise TR-matched MT images when processing the ihMT subtraction image from Eq. 1 (i.e. a dual frequency MT image subtracted from a single frequency image with identical TR).

In order to select data from the whole ihMT dataset, all instances of each of the four MT weighted images used in Eq. 1 were considered along with their corresponding actual TR value:

[ , ] would be considered for the first MT weighted image acquisition, [ , ] would be considered for i th acquisition. Then the algorithm considered a single saturation MT image, e.g. [ , ], and searched for a matched TR within all instances of MT +-and MT -+ . When a match was found, e.g. for the j th image index, [TR " , MT " +-], the two images were coupled together and stored aside. The same process was repeated among the remaining single-and dual-saturation MT dataset until a maximum number of image pairs were generated.

After the whole process, images that could not be paired were discarded from the final dataset and were not considered to form the composite ihMT image, hence leading to SNR loss. Thus, in order to prevent a too restrictive matching process, a tolerance margin δTR was used so that the TR matching condition was defined as TR ± δTR.

Semi-automation for optimal subject-and scan-dependent pairing A total number of 67 δTR values were tested according to the following scheme: δTR varied in 2.5ms-step from 0 to 25ms, then in 5ms-step up to 50ms, in 10ms-step up to 100ms, in 25ms-step up to 500ms and in 50ms-step up to 2000ms. Note that 2.5ms corresponded to the temporal resolution of the physiological data recording, whereas δTR=2000ms was appropriate

1 TR + 1 MT i TR - + i MT 1 TR + 1

MT

to cover the entire range of actual TRs (i.e. not to exclude any data through the pairing algorithm) for the majority of experimental scenarios, including the trigger missing an R-wave because of poor ECG signal quality. For each δTR value, the retrospective pairing algorithm was applied and the corresponding ihMT and ihMTR maps were generated and stored along with the number of paired images from the whole dataset, defined as NEXδTR.

From these multiple available ihMT datasets the selection criterion was designed as a trade-off between negligible residual ihMT signal in CSF (defined as ε, such that |ihMTRCSF|< ε) and highest SNR in SC. For that purpose, ROIs were manually selected on the average S0 image in the CSF and SC (covering levels C2 to C4) areas and used to measure the ihMT metrics (i.e. both ihMTR and SNR) across the entire dataset (identical ROIs were used for all δTR values). ROIs were drawn slightly away from the structure edges to prevent partial volume effects due to blurring from neighboring tissues. Although refined modeling accounting for the coil configuration and acquisition strategy would be necessary for accurate SNR calculation, a simple estimation was used in the current study, defined as the mean ihMT signal in the ROI divided by the noise standard deviation measured in a large signal-free area on the same image, located in the neck posterior side. In practice, three values of ε were systematically evaluated: 0.1, 0.3 and 0.5%. For each ε, the algorithm looped over all δTR values and tested for the |ihMTRCSF|< ε condition; and finally among the datasets satisfying the latter condition, the algorithm selected the dataset providing the highest SNR in the SC ROI.

IhMTR measurements were also performed in the genu of the corpus callosum (CC), considered as a control structure because it provides high ihMT signal while being less prone to CSF motion and partial volume than SC.

Results

Typical images obtained with and without ECG synchronization are shown on Fig. 2.

Whereas unsynchronized MTR data (2d and 2e) do not show obvious signs of motion as compared to ECG-triggered acquisitions (2b), the effect of the CSF pulsatility is particularly strong on ihMT images (2f and 2g), where artifacts are clearly noticeable in the cervical area (green arrows). Conversely, the strategy combining ECG-synchronization and retrospective filtering leads to clean SC ihMT signal free of CSF artifacts (2c).

Typical measurements of MT + , MT -, MT +-and MT -+ across the different repetitions are shown in Fig. 3 for CSF, SC and CC. Using no ECG-synchronization, little signal difference between MT weighted signals was present in CSF (3a), consistent with expected zero-ihMT effect in this structure. SC MT signals (3b) demonstrated poor separation of the singlefrequency data from the dual-frequency data, consistent with CSF pulsatility effects.

Conversely, in CC, a clean separation of these data was observed (3c), emphasizing little CSF motion near this structure. Image acquisition synchronized with ECG led to variable TR and hence to different T1-weightings, which translated into additional signal recovery for increasing TR in CSF (3d). Noteworthy little signal difference was observed between the MT signal data clouds, again consistent with no ihMT effect expected in CSF. For SC, the signal separation between single frequency and dual frequency MT images was clearly evident (3e). Of interest, the mild signal recovery with increasing TR in SC suggested partial volume effects with long T1 components, presumably from CSF. Lower signal variations with TR were obtained in CC (3f), consistent with reduced partial volume effects. Besides, the intensity of the ihMT effect, which corresponds to the signal difference between single-and dual-frequency saturation, appeared similar for CC for both triggered and non-triggered acquisitions, again consistent with low pulsatile motion in this area.

Requirement for the TR-matching process may be instinctively appreciated from Fig. 3, and the automation procedure is illustrated in Fig. 4 where two particularly striking cases are displayed. Whereas no ihMT signal is expected in the CSF compartment, ECG rate fluctuations may induce a strong residual ihMT signal in CSF when no TR-matching is used. The quasirandom nature of these fluctuations makes the artifactual CSF signal either positive (as in (4d))

or negative (as in (4f)), in relation to the average TR of the single-and dual-saturation MT image datasets. The impact of the δTR value on the resulting ihMTR in CSF and SC, on the number of preserved NEX and on the ihMT SNR in SC is demonstrated on the curves. In Fig.

(4a) (positive CSF bias), starting from the right side, (δTR=1000ms, no data exclusion: NEXδTR = total number of NEX = 120), the algorithm progressively discards non TR-matched data-pairs as δTR is shortened, making the CSF ihMTR decrease and tend towards zero. The SC ihMTR progressively decreases as well towards a stable value as δTR is shortened, showing a similar trend as CSF and suggesting partial volume effect. For very short δTR (e.g. <50ms here) the SC ihMT SNR drops because of high number of excluded pairs. Using a threshold value of ε = 0.3% for the CSF signal, the algorithm converges toward the optimal dataset corresponding to δTR=50ms and NEXδTR= 110. The corresponding ihMTR map demonstrated improved conspicuity of the SC and brain tissue (4c). Similar behavior of the TR-filtering algorithm occurs for the negative CSF bias case (4b). The CSF ihMT signal progressively increases towards zero as δTR decreases, and the SC ihMTR converge to a baseline value, again consistent with CSF partial volume. For very short δTR (e.g. < 20 ms here) the actual value of ihMTRCSF was not stable, consistent with poor SNR due to excessive data exclusion during the pairing process. Using ε = 0.3%, the optimal dataset corresponded to δTR=35ms and NEXδTR= 92. Similar to the positive case, the ihMTR map resulting from this correction greatly improved the SC delineation (4e) compared to that of the unpaired data, where the strong negative signal from CSF contaminates the cord signal. Finally, one can appreciate the similarities between the two TR-matched datasets in (4c) and (4e) as well as the strong differences between the unpaired ones in (4d) and (4f).

Comparison between acquisitions performed using different nominal TR values further demonstrated the efficiency of the proposed TR-filtering strategy for various T1-weighting conditions (Fig. 5). Regardless of the intensity of CSF T1-recovery effects (strong in (5a), moderate in (5b), and very low in (5c)), the best TR-matched images ((5g)-(5i)) demonstrated a similar quality to that of the non TR-matched full relaxation image (5f), here taken as a reference due to almost constant CSF signal for TR > 10000ms. Conversely, without TRmatching and for short nominal TR, the CSF residual signal greatly decreased the image quality ((5d) and (5e)). Fig. 6 and 7 summarizes the main quantitative findings observed on the 3 healthy volunteers for all the investigated sets of parameters. Results are indicated as mean value ± standard error of the mean, calculated for the selected ROIs over the three volunteers. Panels (6a) and (6b) focus on ihMTR metrics obtained without and with the TR-matching strategy. For all TR cases the ihMTR values measured in SC with 8mm slice thickness were slightly lower (red curves) than that with 5mm slice thickness (blue curves). When no TR-matching was used (6a), the CSF ihMT signal measured for nominal TRs < 4s was not stable, demonstrating strong signal fluctuations (very large error bars, up to 2.7% absolute scale), and causing increased variability of the SC data as well (relatively large error bars, up to 1% absolute scale). Note that the error bars correspond to the dispersion of the ihMTR values across volunteers, hence emphasizing poor robustness. Conversely, with the TR-matching strategy (6b) and the CSF ihMT signal constrained to 0±ε%, variations of the ihMTR signals greatly decreased in the SC (maximum error bar of 0.5%), thus improving ihMTR measurement stability across all nominal TRs.

The SNR efficiency, defined as the SNR of ihMT measurement obtained per total scan time (i.e. SNR/ $ 𝑇 $%& ), is presented in (6c). Despite data-pair exclusions, the TR-matching strategy combined with the use of shorter TR values generally allowed greatly improved SNR efficiency (roughly doubled from TR = 10000ms to TR = 3500ms for 5mm slice thickness, blue curve). Of interest, data showed low inter-subject variations down to a TR of 3500ms, whereas stronger variability occurred for 2500ms nominal TR.

Finally, Fig. (7a) compares ihMT ratios measured for the different phase encoding

conditions (A-P, A-P +50% phase oversampling and S-I) for both SC and the reference CC structure. Whereas CC measurements were similar for the three PE conditions, SC measurements in the S-I direction (ihMTR ≈ 7.8 ± 0.3%) demonstrated a clear increase as compared to those in the A-P direction both without (ihMTR ≈ 6.3 ± 0.4%, p<0.05) and with 50% PE oversampling (ihMTR ≈ 5.8 ± 0.4%, p<0.01).

Discussion

Whereas non ECG-triggered SC acquisitions may be exploitable on some rare occasions (e.g. when the subject heart rate is particularly stable and synchronous with TR over the whole scan time), a prohibitively strong variability of SC signal was often observed over multiple NEX acquisitions, hence driving the need for an alternative SC ihMT imaging strategy. Such a strong sensitivity to motion is a consequence of the image combination (i.e. subtraction) required to generate the composite ihMT images. Overall, the hereby proposed correction method demonstrated good efficiency with ECG-synchronization able to remove CSFpulsatility artifacts (Figs. 2 and3) and TR-matching process able to remove undesirable residual CSF signal occurring for incomplete CSF signal recovery (Figs. 456). Although the latter feature moderately affected the quantification of ihMTR values within cervical SC on average (Figs.

(6a) vs. (6b)), the conspicuity of the SC was strongly affected on some occasions (Fig. 4) and the metric variability increased when no TR-matching was used (large error bars in (6a)), consistent with partial volume with CSF. Besides, an optimal protocol should provide high fidelity signal for both SC and CSF and good robustness, especially in the clinical context where accurate depiction of the SC morphology is required. More generally, the TR-matching approach could be useful for other SC imaging modalities based on data subtraction such as MT asymmetry (e.g. CEST), arterial spin labeling and fMRI.

The proposed method uses a dedicated TR-matching algorithm to choose relevant data pairs from the individually MT weighted images to form the composite ihMT image. A proper record of the ECG trace, containing appropriate sequence landmarks, is necessary to allow retrospective extraction of the actual TR values of each image instances. Looping over multiple tolerance values δTR for TR-matching ensures automatic selection of the ihMT composite dataset that maximizes SNR in SC while ensuring minimal ihMTR in CSF (0±e%). The flexibility offered by this automatic selection, rather than using a fixed δTR value for all subjects, ensures individual best SNR conditions independent of the ECG variations and of the actual CSF signal bias. Hence when the TR-matching correction was not necessary to obtain negligible ihMT signal in CSF (e.g. Fig. (6a), data corresponding to 5mm slice and 3500ms nominal TR), the algorithm used the whole dataset to compose the ihMT image. Conversely when strong residual CSF signal was present on the non TR-matched ihMT image, a rather short δTR would be automatically selected, corresponding to a certain amount of data exclusion, yielding the optimal SNR for this experimental scenario as explicitly searched for by the selection algorithm (Fig. 4).

Besides the three ROIs that have to be defined manually (SC, CSF and noise) the only other user-defined parameter is the value of ε, i.e. the maximum allowed residual CSF signal.

For all acquired data, three different ε values were systematically tested: 0.1%, 0.3% and 0.5%. Whereas ε = 0.1% appeared too restrictive in 33% of the processed datasets (failure in finding a value of δTR fulfilling the condition |ihMTRCSF| < ε), conditions ε = 0.3% and ε =0.5% yielded satisfying results with respectively only 11% and 4% failures. Of interest, these latter failures all occurred for 50% PE oversampling, hence suggesting that signal blurring may be the cause of such outlier behavior. Indeed, if the residual signal in the CSF area corresponds to actual ihMT signal coming from neighboring CNS tissues, it would not be possible to constrain this signal to zero and the only way to force the algorithm to converge would be to increase ε.

To summarize, all presented quantitative measurements were obtained using a default ε of 0.3%, which was possibly raised to 0.5%, or 0.7%, when necessary and we recommend using similar settings for the general use of this algorithm. Surprisingly, the ihMTR values measured in SC were significantly lower for the A-P PE direction as compared to the S-I direction (Fig. 7a). The readout-module weighting difference due to 50% phase oversampling (547ms readout duration and 31ms effective TE compared to 365ms readout and 20ms TE without oversampling) could be discarded as a hypothetical explanation for this difference because similar (within error bars) ihMTRs were measured on SC when comparing the A-P direction with and without oversampling. Conversely CSF partial volume effects could be hypothesized as a plausible cause of the underestimation of ihMTR in the SC for the A-P PE direction. Indeed, even though the CSF ihMT signal is negligible, CSF would increase the reference S0 signal (especially with a T2-weighted HASTE readout), hence leading to underestimation of the ihMT ratio of the contaminated tissue. As already pointed out, partial volume effects could be appreciated on Fig. 4 where the dynamic of ihMTR vs. δTR measured on SC follows the one from CSF. Moreover these effects were also evidenced in the signal recovery observed in SC (see slope in Fig. 3e) since T1-recovery is not expected over this long TR range (TR≈4000ms) for CNS tissue at 1.5T [START_REF] Rooney | Magnetic field and tissue dependencies of human brain longitudinal 1H2O relaxation in vivo[END_REF], and these recovery effects were stronger in the A-P PE direction as compared to S-I (steeper slope in Fig. 7b). Whereas slice profile imperfections could explain some partial volume contribution, it remains unclear why CSF would "blur" onto SC for relatively short readout (300 to 500ms) with A-P PE, especially considering that data have been Hamming-filtered in k-space to prevent non-local signal contamination. However the extension of the cord in the S-I direction likely explains why partial volume with CSF would be more pronounced on data acquired in A-P PE.

Finally, measurements performed in the CC, less prone to CSF contamination, yielded similar ihMTR values for the three PE conditions (7a). Compared to SC, reduced and hardly noticeable T1-recovery effects were observed in CC for S-I and A-P PE, respectively (see slopes in 3f and 7c), again supporting the CSF signal contamination hypothesis in SC.

A nominal TR value of 3500ms appears to be optimal to maximize both SNR efficiency (for reduced scan time) and robustness, characterized by low inter-subject variability. The global trend of increased SNR efficiency with shorter TR (Fig. 6c) is definitely due to a more efficient use of the timing intervals within the sequence, made possible by proper CSF signal suppression. However, for TR=2500ms the SNR efficiency show strong signs of instability evidenced by larger standard errors for both TR-matched and non TR-matched data. Two combined reasons could explain this phenomenon. On one hand, shorter TR leads to more intense quasi-random value of the CSF ihMTR bias signal for the non TR-matched data, which in turn may affect the SC signal through partial volume effects. On the other hand, the TRmatching algorithm preserves only a finite number of data pairs especially for short TRs, which is also in essence an arbitrary process, and this may translate into strong SNR loss.

Overall it was found that the slice thickness had some impact on the derived metrics measured from C2 to C4 (slightly higher ihMTR values obtained for 5mm slice), likely due to the GM/WM distribution within the cord [START_REF] Taso | Construction of an in vivo human spinal cord atlas based on high-resolution MR images at cervical and thoracic levels: preliminary results[END_REF]. Although one would like the ideal single-slice sagittal protocol to cover the whole cord extension (about 10-15mm in R-L dimension), hence prescribing a rather thick slice, it is important to keep in mind that typical HASTE slice profiles are not very sharp. Particular care was deployed to prescribe the slices in this study, using multiple axial 𝑇 ' * -weighted images to center the slice positioning in the middle of the cord throughout the R-L direction. However, it is often difficult to reach a perfect alignment of the cord, either because of the subject positioning within the RF-coil or simply because of the natural curvature of the cord. Hence no conclusions should be drawn here about the amount of WM and GM contributing to the measured ihMT signal and we can simply conclude that both slice thicknesses provided reasonable and stable estimation of the SC ihMTR.

Despite different subject populations and slice settings the ihMTR values measured here were found to be in the same range of previous measurements performed on brain (2) at 1.5T

for similar energy levels (ihMTR of 8.1 ± 0.3% in Internal Capsule and 6.7 ± 0.3% in frontal WM for ETR ≈ 15uT² vs. ihMTR of 7.8 ± 0.3% in SC and 7.5 ± 0.2% in CC for ETR ≈ 17.2uT² in this study). Here the actual ETR value was kept constant and was chosen to allow scanning within regulatory SAR limit for a nominal TR of 2500ms. Higher energy levels could be considered for a longer TR of 3500ms but a large increase of ihMTR is not anticipated [START_REF] Girard | Magnetization transfer from inhomogeneously broadened lines (ihMT): Experimental optimization of saturation parameters for human brain imaging at 1.5 Tesla[END_REF].

Finally, although the proposed methodology has been carefully demonstrated here on SC images acquired in the sagittal plane at 1.5T, the same strategy may be readily applied in the axial plane as recently demonstrated at 3T [START_REF] Taso | Regional and age-related variations of the healthy spinal cord structure assessed by multimodal MRI (diffusion, inhomogeneous magnetization transfer, ihMT)[END_REF], and could be further used to assess the ihMTR distribution within and across the cord. This would for instance be of particular use to further understand the WM-tract microstructure differences (e.g. potential different myelination of sensory and motor tracts) using atlas-based segmentation [START_REF] Lévy | White matter atlas of the human spinal cord with estimation of partial volume effect[END_REF][START_REF] Fonov | Framework for integrated MRI average of the spinal cord white and gray matter: The MNI-Poly-AMU template[END_REF][START_REF] Taso | A reliable spatially normalized template of the human spinal cord -Applications to automated white matter/gray matter segmentation and tensor-based morphometry (TBM) mapping of gray matter alterations occurring with age[END_REF] and to assess in vivo and non-invasively such features. This would also be very appropriate to go further in the description of the WM impairment in degenerative pathologies such as myelopathy [START_REF] Taso | Cervical myelopathy patient follow-up after decompressive surgery using diffusion tensor imaging (DTI) and inhomogeneous magnetization transfer (ihMT): preliminary application and results[END_REF], MS, SLA and to evaluate potential therapeutic strategies.

Conclusion

This study aimed at providing a fast and robust protocol for ihMT imaging of the SC, free of CSF pulsatility artifacts. ECG synchronization allowed acquiring stable images with respect to the CSF pulsatile motion, but resulted in variable TRs, which translated into a biased residual ihMT signal from CSF that in turn affected SC delineation and ihMT metrics. Fully relaxed acquisition (TR>10s) may be used to solve the problem, though at the expense of a poor SNR efficiency. An alternative method based on retrospective data pairing obtained with a TRmatching semi-automatic algorithm was proposed, allowing improved SNR efficiency and faster acquisition, while offering robust ihMT quantification and exquisite image quality of the SC. This opens great possibility for further in vivo characterization of SC physiopathology using MRI. For long δTR (e.g. >800ms here) the whole dataset is considered to form the composite ihMTR images, hence leading to the biased images shown on the bottom-right side of each panel, (d) and (f). As δTR is shortened, selective TR-matched (± δTR) pairs of MT images are automatically selected by the algorithm and the corresponding ihMTR measured in CSF is progressively constrained to 0±ε %. The best TR-matched datasets, automatically chosen by the algorithm (indicated by the blue vertical lines), result from a compromise between the lowest ihMTR in CSF and the highest SNR in the SC, leading to the unbiased ihMTR images shown on the bottom-left side of each panel, (c) and (e). Note that case 1 was acquired with PE set to S-I and a nominal TR of 2500ms, whereas PE was set to A-P and the nominal TR was 3500ms for case 2.

Figure captions:

Figure 1 :

 1 Figure 1: Effect of k-space windowing demonstrated on T2 weighted HASTE SC data. Low spatial resolution of SC images leads to significant Gibb's ringing in the SC area: 2-times zerofilled images(i.e. boxcar window (a)) demonstrate signal oscillations in the SC whereas Hamming windowed data (b) improves MR signal fidelity of SC by removing CSF secondary lobe contamination, though at the expenses of signal blurring. Theoretical expectations of the ideal point spread function (at infinite T2) for Cartesian sampling using both windows are shown in (c).

Figure 2 :

 2 Figure 2: Effect of ECG synchronization on SC data. (a) Reference ECG triggered anatomical scan (T2w FSE) and corresponding MTR (b) and ihMTR (c) map obtained using optimized protocol (ECG synchronization and retrospective data pairing for ihMT). Examples of nontriggered acquisitions acquired on two different volunteers are displayed in (d) -(g). Whereas MTR data do not evidence strong signs of motion artifacts in the cord area (d) and (e), the ihMTR data reconstructed from the same datasets show particularly strong CSF pulsatile artifacts ((f) and (g), green arrows), consistent with increased motion sensitivity of an imagesubtraction based technique. ECG synchronization combined with the retrospective bias correction algorithm recovers exquisite SC image quality (c).

Figure 3 :

 3 Figure 3: Scatter plots of all individual MT weighted signals (single-and dual-frequency saturation) measured during one full ihMT acquisition in CSF, SC and CC. The signals are plotted vs. the NEX index for the non-triggered acquisition in (a), (b) and (c) and vs. the actual TRs for the triggered acquisition in (d), (e) and (f).For each structure, the MR signal was normalized by the average MR signal measured over all instances of the four kinds of MT weighted images. Presented data were acquired with an 8mm slice and phase encoding in the S-I direction. Steep T1 recovery is evidenced for CSF over this TR range in d), and modest recovery occurs in the SC and CC data in (e). Linear regression plots and corresponding slopes are indicated to help appreciate the T1 recovery effects. For CSF, single-and dual-frequency data clouds are overlapping, consistent with no ihMT effect. For SC, non-triggered acquisition (b) yields partially overlapping data clouds, consistent with pulsatile artifact and partial volume from CSF, whereas triggered acquisition shows clear separation of single-and dual-frequency data clouds (e). For CC, good separation of single-and dual-frequency data clouds is observed in both cases, evidencing strong ihMT effect, and illustrating less motion sensitivity. Examples of ROIs used for SCF, SC and CC measurements are shown on the S0 image at the center.

Figure 4 :

 4 Figure4: Principle of the automatic TR-matching algorithm illustrated on two particularly demonstrative datasets. Because of the random nature of ECG rate fluctuations the residual CSF ihMT signal that is observed when no TR-matching filtering is used may be arbitrarily positive, (a) and (d), or negative, (b) and (f). For long δTR (e.g. >800ms here) the whole dataset is considered to form the composite ihMTR images, hence leading to the biased images shown on the bottom-right side of each panel, (d) and (f). As δTR is shortened, selective TR-matched (± δTR) pairs of MT images are automatically selected by the algorithm and the corresponding ihMTR measured in CSF is progressively constrained to 0±ε %. The best TR-matched datasets, automatically chosen by the algorithm (indicated by the blue vertical lines), result from a compromise between the lowest ihMTR in CSF and the highest SNR in the SC, leading to the unbiased ihMTR images shown on the bottom-left side of each panel, (c) and (e). Note that case 1 was acquired with PE set to S-I and a nominal TR of 2500ms, whereas PE was set to A-P and the nominal TR was 3500ms for case 2.

Figure 5 :

 5 Figure5: Demonstration of the efficiency of the TR-matching algorithm for three different nominal TRs: 2500ms, 3500ms and 10000ms (i.e. full-relaxed CSF dataset). CSF signal recovery over the actual TR ranges are strong (a), moderate (b) and very low (c), hence yielding to strong, moderate and negligible CSF signal bias on the corresponding unfiltered ihMTR images (no TR-matching, (d)-(f)), respectively. Conversely, in all cases the best TR-matched ihMTR images ((g)-(i)) show negligible CSF signal (as in the unfiltered full relaxed dataset shown in (f)) hence providing high-fidelity ihMTR maps over the cervical SC. ε = 0.3% was used here. The fully relaxed acquisition yielded ihMTRCSF ≈ 0.45% with no TR-matching in (f).

Figure 6 :

 6 Figure 6:Quantitative evaluation of the TR-matching strategy in terms of ihMTR and SNR efficiency. While the CSF signal bias was particularly strong (a) for short nominal TRs (e.g. 2500ms) and led to instability of the SC ihMTR (relatively large error bars), the TR-matching strategy satisfactorily constrained the CSF to 0±ε % and reduced SC data inter-subject variability for short TRs (b). Note that both SC and CSF ihMTR are particularly stable over the whole TR range in (b). The corresponding ihMT SNR efficiency (c) is a quantitative index of the SNR improvement and consequently shorter acquisition time made possible by the use of short TRs and allows appreciating SNR losses due to data-pair exclusion (see text for more details). All data in (a)-(c) were acquired with S-I PE direction. A target ε of 0.3% was used for all data (see text for more details).

Figure 7 :

 7 Figure 7: Comparison of the ihMTR measured in SC and CC for both S-I and A-P PE directions (a) suggest that S-I direction should be preferred to avoid residual bias leading to strong underestimation of SC metrics. Conversely the CC ihMTR are in very good agreement for both PE directions. Panels (b) and (c) display scatter plots similar to Fig.3(e) and (f) except that the phase encoding was set to A-P with 50% phase oversampling. Linear regression plots and corresponding slopes are indicated to help appreciate the T1 recovery effect and compare with S-I data from Fig.3. A much steeper T1 recovery is observed for the SC ROI in (b) as compared the plot from Fig.3(e) suggesting stronger partial volume effects from CSF. All data were acquired using TR=3500ms, a 5mm slice and a target ε of 0.3%.

Acknowledgments

This study was supported by IRME. V.H.P. received support from A*MIDEX (n°ANR-11-IDEX-0001-02) funded by the French Government "Investissements d'Avenir" program.

The authors thank Elisabeth Soulier and Véronique Gimenez for help with the volunteer management.