Nils Krah
email: nils.krah@creatis.insa-lyon.fr

Catherine Therese

T Quiñones

Jean-Michel M Létang

Simon Rit

C T Quiñones

Scattering proton CT

à la diffusion de documents scientifiques de 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

The use of charged particles in medical imaging has been proposed long ago by [START_REF] Cormack | Representation of a Function by Its Line Integrals, with Some Radiological Applications[END_REF] and [START_REF] Koehler | Proton Radiography[END_REF] and later by [START_REF] Cormack | Quantitative proton tomography: preliminary experiments[END_REF], [START_REF] Hanson | Computed tomography using proton energy loss[END_REF], [START_REF] Hanson | Proton computed tomography of human specimens[END_REF]. While initially discarded in disfavour of x-ray computed tomography (CT), proton CT has seen renewed interest more recently in the context of proton therapy. We refer to [START_REF] Schulte | Proton CT for Improved Stopping Power Determination in Proton Therapy, invited[END_REF], [START_REF] Poludniowski | Proton radiography and tomography with application to proton therapy[END_REF], and [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF] for recent reviews. Proton CT in its common implementation directly provides a relative stopping power (RSP) image of the patient -a quantity needed by treatment planning systems. This is expected to be an advantage over the currently used x-ray planning CT images which need to be translated into relative stopping power by means of an approximate conversion.

In this work, we focus on so-called "single tracking" or "list-mode" proton CT scanners, i.e. systems which record the kinematic properties of individual protons. Typically, these systems measure each proton's position and direction upstream and downstream of the imaged object, as well as some quantity indicative of the energy lost by the proton while traversing the object. Independently of the detector hardware employed to achieve this, the reconstructed quantity is the RSP. We will henceforth refer to this modality as energy-loss proton CT.

The reconstruction problem in proton CT is in principal similar to the one in x-ray CT in that it involves line integrals across the imaged object. A peculiar trait of proton CT is that proton trajectories are stochastic and not straight because protons undergo multiple Coulomb scattering (MCS). This is commonly accounted for by estimating the most likely path (MLP) for each proton and performing line integrals along the resulting curvilinear lines [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF][START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF][START_REF] Erdelyi | A comprehensive study of the most likely path formalism for protoncomputed tomography[END_REF][START_REF] Collins-Fekete | A theoretical framework to predict the most likely ion path in particle imaging[END_REF][START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF].

In addition to energy-loss proton CT, other contrast mechanisms have been proposed, which exploit different types of interaction of protons with a medium. In particular, these are attenuation and scattering proton CT. The former measures the reduction in proton fluence after an object due to inelastic nuclear interactions and reconstructs a map of the nuclear attenuation coefficient [START_REF] Quiñones | Filtered back-projection reconstruction for attenuation proton CT along most likely paths[END_REF]. The latter estimates the angular dispersion of protons due to MCS in the object and reconstructs a parameter describing MCS, e.g. radiation length [START_REF] Bopp | Proton computed tomography from multiple physics processes[END_REF][START_REF] Bopp | Quantitative proton imaging from multiple physics processes: a proof of concept[END_REF][START_REF] Taylor | An experimental demonstration of a new type of proton computed tomography using a novel silicon tracking detector[END_REF]). An appealing practical aspect of scattering proton CT is that the imaging system is simpler than in energy-loss proton CT: it only requires tracking devices (to measure a proton's direction upstream and downstream of the object), but no detector to measure a proton's residual energy. It has also been observed that scattering proton CT might provide higher contrast images than energy-loss proton CT [START_REF] Plautz | 200 MeV Proton Radiography Studies With a Hand Phantom Using a Prototype Proton CT Scanner[END_REF].

At the same time, there is currently very limited literature available about reconstruction algorithms for scattering proton CT which use the MLP. [START_REF] Quiñones | Filtered back-projection reconstruction for attenuation proton CT along most likely paths[END_REF] present a filtered backprojection reconstruction method for attenuation proton CT which incorporates the MLP via distance driven binning, as initially proposed by [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF]. [START_REF] Bopp | Quantitative proton imaging from multiple physics processes: a proof of concept[END_REF] reconstruct scattering proton CT images via iterative reconstruction yet without considering the MLP. [START_REF] Taylor | An experimental demonstration of a new type of proton computed tomography using a novel silicon tracking detector[END_REF] have used a backprojection-then-filtration algorithm [START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojectionthen-filtering approach[END_REF], but provide relatively few details. In this work, we present a distance-driven binning reconstruction for scattering proton CT which takes into account the MLP.

Furthermore, we analyse systematically the potential performance of scattering proton CT and compare it with the more common contrast mechanism , i.e. energy-loss proton CT, as well as attenuation pCT. Such a study is, to our knowledge, also missing in the literature so far. In particular, we investigate artefacts in scattering proton CT and their origin. We analyse the intrinsic lower limits in terms of image noise and finally quantify the typical signal-to-noise ratio, contrast-to-noise ratio, and spatial resolution achievable in scattering proton CT. We perform our study partly based on analytical models and partly based on Monte Carlo simulations.

Materials and Methods

Rationale of scattering proton CT

The stochastic interaction of protons with the Coulomb field of nuclei in a medium results in multiple deviations from the protons' original direction. This phenomenon is referred to as multiple Coulomb scattering (MCS) [START_REF] Gottschalk | On the scattering power of radiotherapy protons[END_REF]. Let the proton beam be aligned with the z-axis, then the quantities of interest are the angular deviations θ x and θ y as well as the lateral deviation t x and t y in the two planes perpendicular to the beam axis, respectively. Henceforth, we will drop the index x or y for simplicity and tacitly assume that angles are projected onto one of the two transversal planes.

We denote with A ≡ θ 2 the variance of the angular distribution in some depth z and define the scattering power as the rate of change of the angular variance with depth,

T (z) = dA(z) dz ⇔ A(z) = z 0 T (z )dz + A 0 . ( 1
)
This definition is in analogy to the (energy) stopping power dE(z)/dz. We will henceforth consider A 0 = 0, without loss of generality. The rational behind scattering proton CT is to measure A or a related quantity and reconstruct the scattering power T of the traversed tissue.

It is important to note, however, that T is not a purely material specific local quantity, but also depends on the proton energy E(z) at depth z which in turn depends on the initial energy E in and the type and amount of tissues upstream of z. In this work, we explore if and how a scattering power map could be reconstructed using analytical reconstruction methods.

Relative scattering power

In this section, we present how analytical methods used in energy-loss pCT could be applied to scattering pCT. To start off, we recall that the proton stopping power dE/dz depends on the proton energy [START_REF] Paganetti | Proton Therapy Physics[END_REF]. In fact, this reciprocal dependency is what leads to the sharp Bragg peak in depth dose profiles and what gives proton therapy its characteristic dosimetric properties. However, the ratio of stopping powers in human tissue relative to that of water is almost independent of energy [START_REF] Arbor | Monte Carlo comparison of x-ray and proton CT for range calculations of proton therapy beams[END_REF]. It can therefore be considered a local tissue property which upon integration along a projection line gives the water equivalent path length making energy-loss pCT a well-defined tomographic problem. We will derive the tomographic problem of scattering pCT based on the concept of relative scattering power.

In this section, we consider as spatial variable the depth z along a certain line of projection under a certain projection angle, starting at z = 0 where the protons enter the object. We will later generalize the variable z to account for the most likely path, but wish to keep notation as simple as possible here.

We use the radiation length X 0 to parametrize the degree of MCS in a medium, so that the scattering power dA/dz = T (X 0 (z), E(z)) depends on the proton energy E(z) at z and on the radiation length X 0 (z) [START_REF] Gottschalk | On the scattering power of radiotherapy protons[END_REF]. The quantity X 0 (z) is purely local, i.e. it only depends on the local tissue properties, while the term E(z) depends on the material traversed by a proton before reaching z and therefore on the projection angle.

We observe that

f : z → A(z) and
(2)

g : z → E(z) (3)
are always monotonous functions because the proton energy always decreases and the angular dispersion always increases. Therefore, the mapping

h ≡ g • f -1 (4)
is always bijective. In more practical terms, E = h(A) is the energy of an initially co-linear proton beam in a medium at a depth where its angular dispersion is A. We can therefore define the scattering power as a function of angular variance, τ (X 0 (z), A(z)) ≡ T (X 0 (z), h(A(z))).

(5)

We note that the function h is not unique because A(z) and E(z) depend on the composition of the medium and the line of projection along which z is taken, unless the medium is homogeneous. With h w the mapping in water, we define the scattering power in water as a function of A,

τ w (A) ≡ T (X 0,w , h w (A(z))) (6)
and the relative scattering power

δ(X 0 (z), A(z)) = τ (X 0 (z), A(z)) τ w (A) . ( 7
)
Combining equations 1 and 7, we obtain the integral expression

G(A out ) = Aout 0 1 τ w (A) dA = zout 0 δ(X 0 (z), A(z))dz. (8)
The left-hand side can be interpreted as the water equivalent scattering path length G, i.e. the thickness of a block of water needed to generate the angular dispersion A out . This is in analogy to the water equivalent path length in energy-loss pCT. Reconstructing δ from G is a typical tomography problem.

Image reconstruction along most likely paths

The true trajectory of a proton through the imaged object is unknown due to the stochastic nature of multiple Coulomb scattering, but the most likely path (MLP) can be estimated using well-established methods [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF][START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF][START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF][START_REF] Collins-Fekete | Pre-treatment patient-specific stopping power by combining list-mode proton radiography and x-ray CT[END_REF][START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF].

In particular, we assume a point source moving on a circular trajectory around the isocenter. For a given source position, i.e. projection angle, we define the orthonormal basis vectors u, v, w, where w defines the source-to-isocenter direction and u and v represent the axes of the detector plane. We denote the MLP of the i-th proton with Γ i and describe it as function of the depth w along w, i.e. Γ i (w) • w = w, as do all MLP models.

In this notation, equation 8 formally becomes

G(A out,i ) = A out,i 0 1 τ w (A) dA = Γ i δ( Γ i (w))dw. (9)
It relates the relative scattering power integrated along the MLP of the i-th proton to the water equivalent scattering thickness of the object along this MLP. It is worth noting a fundamental difference here however: while in energy-loss pCT, an energy is measured for each proton, the angular dispersion A out can only be estimated from a set of proton histories. In practice, the quantity G(A out,i ) is thus not available proton per proton, as the equation suggests, but only for groups of protons. Estimation of A out (and G(A out )) occurs when binning the data, as we will explain in the following. We use the distance-driven binning technique initially proposed for energy-loss pCT [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF] and later adapted for attenuation pCT [START_REF] Quiñones | Filtered back-projection reconstruction for attenuation proton CT along most likely paths[END_REF]. The basic idea of this method is to bin the measured data into a series of virtual projections at different distances from the sources while accounting for the MLP. The subsequent filtered backprojection (FBP) uses for each voxel the projection images which most appropriately matches the voxel's position in the image. In the following, we adapt the method to scattering pCT.

Dependence of relative scattering power on angular variance 2 MATERIALS AND METHODS

First off, let j ∈ J ⊂ Z 2 be a set of spatial indices corresponding to the binning grid of a virtual projection and y ∈ R 2 the spatial location in the projection plane. Then the pixel indicator function h: R 2 → R is defined as

h j (y) = 1 if y ∈ R 2 is in pixel j, 0 otherwise. (10)
The measured quantity of interest is the difference between the entrance and exit direction of each proton, ∆θ i = θ out i -θ in i . We treat MCS in the two directions perpendicular to the beam as independent statistical processes [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] and use the average of the squared difference of both projected angles, [(∆θ u i ) 2 + (∆θ v i ) 2 ]/2, to estimate the variance attributed to a pixel. The distance-driven variance Ãj (w) of the angular distribution of protons crossing pixel j is defined as,

Ãj (w) = i∈I h j (u i (w), v i (w))[(∆θ u i ) 2 + (∆θ v i ) 2 ]/2 i∈I h j (u i (w), v i (w)) , ( 11
)
where we have used the fact that ∆θ = 0 on average. The tilde symbol indicates that à is only an estimate of the variance A.

The distance-driven binned variances Ãj (w) are converted into scattering WEPL G through an empirical lookup table. In particular, we ran a Monte Carlo simulation in GATE/Geant4 at high statistics (10 6 primaries) of a perfectly collimated (θ in = 0) pin-like proton beam impinging onto boxes of water with varying thickness and recorded each proton's exit angle. For each thickness, we estimated the variance A out = θ 2 out and obtained the numerical relationship G LUT (A out ) by interpolation. The distance-driven binned scattering WEPL is then

g j (w) = G LUT ( Ãj (w)). (12)
A distance-driven projection g s j (w) is computed for each source position s and the δ-map is reconstructed using the same algorithm as in [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF], i.e. using the FDK algorithm [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF] adapted with a voxel-specific backprojection selecting the distance w according to the distance between the voxel and the exit detector for each source position.

Dependence of relative scattering power on angular variance

To which extent the reconstructed δ-map reflects the material specific local properties (to be specified further down) is related to how much δ(z, A(z)) depends on A. We recall that the same kind of argument holds for energy-loss where the relative stopping power is approximately independent of E [START_REF] Arbor | Monte Carlo comparison of x-ray and proton CT for range calculations of proton therapy beams[END_REF]. To study this aspect, we analysed the relation δ(A) numerically in different kinds of media. For computational simplicity, we relied on an analytical model rather than recurring to full Monte Carlo simulations.

In particular, we used the analytical model of MCS proposed by [START_REF] Highland | Some practical remarks on multiple scattering[END_REF] and modified by [START_REF] Gottschalk | Multiple Coulomb scattering of 160 MeV protons[END_REF] which accounts for the protons' energy loss while traversing an object. The variance A(z) of the angular distribution in some depth z in a homogeneous medium of radiation length X 0 is given by

A(z) ≡ θ 2 (z) = 1 + 1 9 log 10 z X 0 2 Ω 2 0 z 0 1 p 2 (z )v 2 (z ) 1 X 0 dz , (13)
with Ω 0 = 14.1 MeV an empirical constant and p and v the protons' momentum and velocity, respectively. We note that several other analytical MCS models exist and refer to Gottschalk (2010) for a review.

Evaluating the term

1 p 2 (z)v 2 (z) = (E(z) + E p ) 2 c 2 (E(z) + 2E p ) 2 E 2 (z) , (14)
with E p the proton's rest energy, requires the kinetic energy E(z) as a function of depth. This can be readily obtained numerically from a simple Monte Carlo simulation of protons penetrating the medium of interest or by numerically integrating tabulated stopping powers, e.g. from NIST's PSTAR database [START_REF] Berger | ESTAR, PSTAR, and ASTAR: Computer Programs for Calculating Stopping-Power and Range Tables for Electrons[END_REF], in the continuous slowing down approximation. The scattering power, i.e. the derivative of equation 13, contains two terms,

T (z) = dA dz (z) = 1 + 1 9 log 10 z X 0 2 Ω 2 0 1 p 2 (z)v 2 (z) 1 X 0 (15) + 2 9z 1 + 1 9 log 10 z X 0 Ω 2 0 z 0 1 p 2 (z )v 2 (z ) 1 X 0 dz ≈ 1 + 1 9 log 10 z X 0 2 Ω 2 0 1 p 2 (z)v 2 (z) 1 X 0 . ( 16
)
The second term stems from the chain rule applied to the logarithmic pre-factor and its magnitude is on the order of 1-3% of the first term. We neglect it in the following. The relationship τ (A) for a given medium can be constructed numerically by evaluating equations 13 and 15 and linking T and A by interpolation via z. The relative scattering power of the medium is δ(A) = τ (A)/τ w (A).

At z = 0 (i.e. A = 0) in equation 15, the relative scattering power τ (0)/τ w (0) = T (0)/T w (0) corresponds to the radiation length ratio X 0,w /X 0 . This would also be true, if the energy-loss in any medium were the same as in water and the energy-dependent terms 1/β 2 p 2 cancelled out (disregarding the logarithmic pre-factor). Therefore, we use the ratio X 0,w /X 0 as material specific property to compare relative scattering powers δ with.

We use will the expressions presented in this section to analyse how much δ depends on A in materials and beam energies of interest for scattering proton CT.

Statistical limitations in scattering proton CT

In this section, we derive an expression to quantify the inherent noise of a scattering pCT image and compare it to energy loss pCT. We disregard noise and artefacts created by the discretization of the reconstruction algorithm and concentrate on noise due to uncertainties in the projections only. For simplicity, we consider a homogeneous water cylinder and the noise at its centre. We tacitly assume that all expressions below are evaluated in the central pixel j centre of a projection plane at depth w centre corresponding to the cylinder centre as [START_REF] Schulte | Density resolution of proton computed tomography[END_REF].

As part of the distance driven binning, we estimate the angular variance A from a finite number of proton histories as

à = (1/2N ) i∈I [(∆θ u i ) 2 + (∆θ v i ) 2 ] (see equation 11),
where N is the number of proton histories in a pixel. This estimate will be intrinsically uncertain even if the angles ∆θ are measured with perfect accuracy. The standard error σ 2 à of the variance à estimated from a finite sample of the distribution is

σ 2 Ã = 1 2N µ 4 - 2N -3 2N -1 A 2 ≈ 1 2N 3 - 2N -3 2N -1 A 2 , ( 17
)
where µ 4 is the fourth moment of the distribution of ∆θ and A its variance (Rao, 1973). The normalisation factor is 1/2N instead of 1/N because it accounts for the fact that ∆θ u and ∆θ v are both used and that they are independent measurements. As the true values of A 2 and µ 4 are unknown, we estimate them from a high-statistics Monte Carlo simulation of protons penetrating a water phantom. The simplified expression on the right-hand side is valid under the assumption that the angles ∆θ are normally distributed (µ 4 = 3A 2 ). This is only approximately true because of large angle electromagnetic scattering leading to non-Gaussian tails [START_REF] Gottschalk | On the scattering power of radiotherapy protons[END_REF] and due to nuclear interactions present in measured data if their filtering is imperfect. Following the work of [START_REF] Schulte | Density resolution of proton computed tomography[END_REF] on the intrinsic noise in energy-loss pCT, we use first-order error propagation to calculate the uncertainties of the scattering WEPL,

σ 2 g = σ 2 Ã dG dA A= Ãout 2 = σ 2 Ã τ 2 w ( Ãout ) , ( 18
)
where we have used the fact that G is the primitive function of A by construction (see equation 9). The effect of the ramp filter on the uncertainty σ 2 g in the central pixel j centre in the projections can be described by a factor 1/12a 2 [START_REF] Gore | Statistical limitations in computed tomography[END_REF] and reconstructing the image via FBP involves halving the sum of the P projections distributed over the full 2π circle (in our case) in constant angular increments ∆α = 2π/P , so that

σ 2 δ = 1 2 P p=1 1 12a 2 σ 2 g ∆α 2 = σ 2 g π 2 6a 2 P , ( 19
)
where we have used the fact that σ 2 g is independent of the projection angle in cylindric objects. Combining equations 17, 18, and 19, we obtain

σ 2 δ = 1 2N μ4 - 2N -3 2N -1 Ã2 out 1 τ 2 w ( Ãout ) π 2 6a 2 P . (20)
We wish to express this uncertainty as a function of dose to the centre of the cylinder along the lines of [START_REF] Schulte | Density resolution of proton computed tomography[END_REF]. To this end, it is useful to introduce the proton fluence into equation 20 above. Let Φ in be the initial fluence of a proton beam impinging centrally onto the cylinder, then the downstream fluence behind the cylinder centre and the fluence at the centre, respectively, are

Φ out = Φ in exp (-2κr) (21)
Φ centre = Φ in exp (-κr) , ( 22
)
where r is the cylinder radius and κ the attenuation coefficient due to nuclear interactions.

Here, we approximate proton trajectories as straight lines disregarding MCS. The number N in equation 20 refers to the number of protons crossing an area a × a in the object's centre given by the pixel side length and the slice thickness (we use cubic voxels) and exiting the object downstream, so that

N a 2 = Φ out = Φ centre exp (-κr) . ( 23
)
Inserting equation 23 into equation 20 and letting (2N -3)/(2N -1) ≈ 1 (because N 1), we have

σ 2 δ = μ4 -Ã2 out π 2 τ 2 w ( Ãout )12P a 4 Φ centre exp (-κr) . ( 24
)
The dose D centre accumulated at the cylinder centre during a full pCT scan is linked to the fluence Φ centre by

D centre = P Φ centre S(E centre ) + κγE centre ρ ( 25
)
where S(E centre ) is the proton stopping power for the proton energy at the centre, ρ is the mass density of the cylinder material and γ=0.65 is the fraction of energy transferred to secondary particles during nuclear interactions [START_REF] Schulte | Density resolution of proton computed tomography[END_REF]. Substituting equation 25 into equation 24 yields the final expression:

σ δ = (μ 4 -Ã2 out ) 1/2 τ w ( Ãout ) π 2 (S(E centre ) + κγE centre ) 12a 4 exp (-κr) ρD centre 1/2 . ( 26
)
If the angles ∆θ are normally distributed, μ4 -Ã2 out = 2 Ã2 out , and the equation above simplifies accordingly.

For comparison, we state the analogous expression for the uncertainty of the reconstructed relative electron density (proportional to RSP) in energy-loss pCT as derived in [START_REF] Schulte | Density resolution of proton computed tomography[END_REF] (here with an additional factor 2 in the denominator to account for halving opposite projections along a full 360 degree acquisition trajectory):

σ ηe = σ Eout S(E out ) π 2 (S(E centre ) + κγE centre ) 6a 4 exp (-κr) ρD centre 1/2 (27)
In a water cylinder with 10 cm radius, E centre ≈ 151 MeV for a 200 MeV proton beam and ≈ 263 MeV for a 300 MeV proton beam. We use the same value κ = 0.01 cm -1 as [START_REF] Schulte | Density resolution of proton computed tomography[END_REF] which is inferred from proton on oxygen reactions [START_REF] Malmer | ICRU Report 63. Nuclear Data for Neutron and Proton Radiotherapy and for Radiation Protection[END_REF]. At 151 MeV, the proton stopping power S(E centre ) in water is about 5.4 MeV cm 2 /g [START_REF] Berger | ESTAR, PSTAR, and ASTAR: Computer Programs for Calculating Stopping-Power and Range Tables for Electrons[END_REF] and the scattering power T w is about 6.5 × 10 -5 cm -1 . At 263 MeV, the values are about 3.8 MeV cm 2 /g and 2.3 × 10 -5 cm -1 , respectively.

Monte Carlo simulations

We performed Monte Carlo simulations using Geant4 [START_REF] Agostinelli | GEANT4 -A simulation toolkit[END_REF] version 10.01.p02 via GATE [START_REF] Sarrut | A review of the use and potential of the GATE Monte Carlo simulation code for radiation therapy and dosimetry applications[END_REF] to evaluate the various aspects of scattering proton CT described in the preceding sections. We used the QGSP_BIC physics list, but retained only those protons for the reconstruction which had not undergone any nuclear interaction. Angular deflection in our data is therefore purely caused by MCS. On experimental data, one would impose a threshold-based filtering, commonly referred to as "3-sigma" cuts, to eliminate events with large angular deviation which are likely to be associated with nuclear interactions [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. Furthermore, we set the proton step size limit to 1 mm throughout this study.

Phantoms

We used a homogeneous water cylinder with 20 cm diameter to study the noise properties of scattering proton CT and a homogeneous aluminium cylinder of 5 cm diameter to illustrate how the dependence of δ on A (see section 2.4) results in a cupping artefact. We selected the Gammex 467 phantom (Sun Nuclear, Melbourne, USA) which has a diameter of 33 cm and includes 17 cylindrical tissue equivalent inserts. Finally, we employed a cylindrical water phantom with cylindrical aluminium inserts of 1.5 cm diameter arranged in a spiral around the centre (see figure 8 in the Results section). All phantoms were implemented digitally in GATE by specifying their material composition and dimensions.

For the part of our study based on analytical models, we used a 20 cm thick block of water with a 2 cm thick slab of bone-like material (SB3 Gammex insert material) at variable depths. This phantom was described analytically by its RSP and radiation length X 0 as a function of depth.

Simulation set-up

The simulation set-up consisted of a mono-energetic fan-beam proton source, the phantom at isocentre, and two trackers (one upstream and one downstream of the phantom). The fan-beam source was positioned 5 m upstream from isocentre and its dimensions at isocentre were adapted to the phantom size: 267×5 mm 2 for the water, aluminium, and spiral phantom, and 377×5 mm 2 for the Gammex phantom, where the longer dimension corresponds to the transverse direction and the shorter one to the direction of the rotation axis. The beam energy was either 200 MeV or 300 MeV, as we will indicate in the Results section. The uniform source fluence was set to yield a dose at the phantom centre of 50 mGy for an entire scan, except for the water cylinder phantom where the dose varied as indicated in the Results section. The ideal trackers recorded each proton's position and direction and the down-stream trackers additionally registered each proton's energy for the energy-loss images. A total of 360 projections were acquired in angular increments of 1 degree and binned into 1 × 1 mm 2 pixels for the water, aluminium, and Gammex phantom, and a finer 0.5×0.5 mm 2 pixels for the spiral phantom which was used to assess spatial resolution.

Image quality metrics

Spatial resolution Following [START_REF] Khellaf | 2D directional ramp filter[END_REF], we estimated the spatial resolution using the spiral phantom described in section 2.6.1. In and around each aluminium insert, we collected the pixel values together with their distance to the insert centre, r. This provided us with an isotropically sampled radial edge spread function (ESF), as in [START_REF] Richard | Towards task-based assessment of CT performance: System and object MTF across different reconstruction algorithms[END_REF], to which we fitted the following model based on the error function (erf):

ESF(r) = a 2 1 -erf r -R √ 2σ + 1, (28)
with a, R, and σ fit parameters. The rationale of using the error function is that image blur in proton CT is dominated by MCS which in turn is approximately Gaussian. As measure for the spatial resolution (in line pairs per mm: lp/mm), we used the frequency f 10% at which the modular transfer function of the ESF decayed below 10% of its maximum value, which is analytically related to σ,

f 10% = ln 10 2 1 πσ . ( 29
)
Noise We determined image noise as the variance of the central pixel value over 110 reconstructed images of the cylindrical water phantom, each based on independently simulated Monte Carlo data. We used this value as reference to compare equation 26 with.

Contrast-to-noise and signal-to-noise We calculated the contrast-to-noise ratio (CNR) and signal-to-noise ratio (SNR) of the inserts in the Gammex phantom (see section 2.6.1). Specifically, we defined a circular region-of-interest (ROI) of half the insert diameter and determined the CNR and SNR as:

SNR = µ ROI σ ROI (30) CNR = µ ROI -µ water σ ROI (31)
where µ ROI and σ ROI are the mean value and root-mean-square-error per ROI, and µ water the mean value of the water-like base structure of the phantom.

Results

Dependence of relative scattering power on A

We analysed the dependence of relative scattering power δ on angular variance A as explained in section 2.4. Figure 1 condenses much of the relevant information: the upper four panels (a-d) for a proton energy of 200 MeV and the lower four (e-h) for 300 MeV. The analysis was done in homogeneous material descriptions corresponding to the different insert materials of the Gammex phantom (colour-coded). The dash-dotted lines refer to liquid water.

Panels a and e show the angular variance A(z) of an initially co-linear pin-like proton beam as a function of depth z and panels b and f the scattering power T (z). Notably, A(z) and T (z) become very steep and non-linear towards the protons' range. This obviously corresponds to shallower depths at 200 MeV than at 300 MeV. Panels c and g show the scattering power τ (A) strongly varies as a function of A. Panels d and h show the relative scattering power δ(A) as a function of A. For materials with a density similar or inferior to water, this dependence is relatively weak, while δ decreases considerably with increasing A for bone-like tissues (CB and SB). For example, the relative scattering power of cortical bone (CB2_50) at a variance of 1.5 × 10 -3 rad 2 (representative of 20 cm of water at 200 MeV proton energy), is about 20% lower than its X 0,w /X 0 ratio.

The left panel of figure 2 shows the reconstructed image of the aluminium cylinder and the right panel depicts the radial profile of reconstructed values. There is a clear cupping artefact: relative scattering power decreases towards the centre of the cylinder while it reaches about the characteristic ratio of X 0,w /X 0 = 4.06 for aluminium near the edges. The fact that central pixels are associated with comparatively large angular variances for which δ is more strongly underestimated explains this observation. We remark that the mass density of aluminium (2.7 g/cm 3 ) is about 1.5 times higher than in the most dense bone insert (SB3), making the cupping artefact particularly pronounced.

Figure 3 shows the analytically calculated relative scattering power in a case where protons traverse the water phantom with bone slab inserts (section 2.6.1) for a proton energy of 200 MeV and 300 MeV. To this end, we made the radiation length in equations 13 and 15 a function of depth, X 0 (z). The dashed coloured lines indicate the average value of δ in the slabs which are systematically lower than the ratio X 0,w /X 0 (black dashed line). The presence of the slabs also causes an underestimation of δ in the water subsequently traversed. At 200 MeV, the degree of underestimation varies considerably depending on the depth of the bone slab, while it is almost independent of the insert depth at 300 MeV.

Image quality: noise, SNR, CNR

As a figure of merit for the image noise, we calculated the expected precision of the reconstructed relative scattering power in the central pixel of a water cylinder, as explained in section 2.5 (equation 26). We compare the results with the relative stopping power precision in energy-loss proton CT employing the formalism of [START_REF] Schulte | Density resolution of proton computed tomography[END_REF]

(our equation 27).

The left panel in figure 4 shows the precision as a function of cylinder diameter for a dose of 10 mGy at the cylinder centre. The right panel depicts how relative scattering power precision varies as the inverse square root of the dose. For comparison, we included the intrinsic noise in attenuation proton CT in an otherwise identical setting, as derived in [START_REF] Quiñones | Filtered back-projection reconstruction for attenuation proton CT along most likely paths[END_REF].

Regardless of this, the most striking conclusion from our figure 4 is that the order of magnitude in precision is vastly different in the two modalities (beware the two separate y-axes): relative stopping power precision is about 1.6% for a 20 cm cylinder while it is about 50% in scattering proton CT. This was also visually apparent in the Gammex and spiral phantom images in figures 6 and 8, respectively. The main reason for this is that estimating the variance of a distribution (in scattering proton CT) is less certain than estimating the mean value (the mean residual energy in energy-loss proton CT). Specifically, the factor Ãout /τ w in equation 26 is 20-40 times larger than the factor σ Eout /S in equation 27.

Figure 5 shows the CNR (upper panel) and SNR (centre panel) of the Gammex phantom inserts derived from the Monte Carlo simulated images shown in figure 6. The attenuation proton CT image was reconstructed as described in [START_REF] Quiñones | Filtered back-projection reconstruction for attenuation proton CT along most likely paths[END_REF]. We have arranged the values as a function of mass density to illustrate this dependence. For inserts with a mass density close to that of the water-like support structure of the Gammex phantom (most abundant in the human body), the CNR in scattering proton CT is below 0.5 and an order of magnitude lower than the CNR in energy-loss proton CT. The SNR is 3-4 in scattering proton CT compared to around 200 in energy-loss proton CT. Overall, both CNR and SNR are lower by a factor of 40-50 in scattering proton CT than in energy-loss proton CT. This is mainly due to the higher noise in the former modality and not due to intrinsically different levels of contrast, which can be seen by the following two observations: For a water cylinder of 33 cm diameter imaged with a beam energy of 300 MeV at 50 mGy dose at the centre, the δ-resolution is about 100% while the RSP resolution is about 2.5% (values not shown in figure 4), i.e. a relative factor of 40. On the other hand, the lower panel in figure 5 provides an indication of the intrinsic contrast in each modality showing X 0,w /X 0 and RSP, respectively. They are essentially identical except for bone-like densities where scattering proton CT promises 50% higher contrast.

Accuracy of relative scattering power δ

Figure 7 shows the accuracy of the reconstructed relative scattering power (red) compared to the material specific ratio X 0,w /X 0 (black). The values were derived from the Monte Carlo simulated images shown in figure 6 by averaging over 330 pixels around the centre of each insert. The inserts are sorted according to their density. For bone-like inserts (IB3, B200, CB2_30, CB2_50, SB3), the reconstructed value of δ underestimates the ratio X 0,w /X 0 . Instead, for the inserts AP6, BR_12, and BRN-SR2, the reconstructed δ overestimates X 0,w /X 0 . This is in line with the observations made earlier in section 3.1 regarding figure 1 where δ(A) decreased as a function of A for bone-like materials and increased for some of the others. The upper right panel in figure 8 shows the mean value of δ per insert (aluminium) in the spiral phantom as a function of distance from the phantom centre. The ratio X 0,w /X 0 is increasingly underestimated for inserts closer to the centre. This is expected because on average over all projections, the L N 3 0 0 (1 ) L N 3 0 0 (2 ) L N 4 5 0 A P 6 (1 ) variance A is larger for central inserts than for peripheral ones so that the discrepancy between δ(A) and X 0,w /X 0 is larger.

A P 6 (2 ) B R 1 2 W a te r_ S o lid (1 ) W a te r_ S o lid (2 ) B R N -S R 2 L V 1 (1 ) L V 1 (2 ) IB 3 B 2 0 0 C B 2 _ 3 0 C B 2 _

Spatial resolution

The lower right panel in figure 8 shows the spatial resolution estimated from the spiral phantom images as a function of distance to the phantom centre. The values were obtained by fitting an error function to line profiles positioned radially around each insert as described in section 2.7. An example profile is shown in the centre right panel. Error bars were derived from the fit uncertainty. The plot compares scattering proton CT using distance-driven binning reconstruction, straight line reconstruction, and energy-loss proton CT with distance driven binning reconstruction.

Spatial resolution ranges from 0.7 lp/mm at the centre to 1.4 lp/mm near the edge. Incorporating the MLP improves spatial resolution by almost a factor 3. Spatial resolution in scattering proton CT appears to be slightly inferior compared to energy-loss proton CT near the phantom centre, although this difference could still be explained by the statistical uncertainty of the data.

Discussion

In this work, we investigated scattering proton CT reconstruction. Specifically, we introduced the relative scattering power in analogy to the relative stopping power in energy-loss proton CT and adapted a distance-driven filtered backprojection algorithm to scattering proton CT. It is worth mentioning that there is some resemblance between our algorithm and the noise reconstruction in [START_REF] Rädler | Two-dimensional noise reconstruction in proton computed tomography using distance-driven filtered back-projection of simulated projections[END_REF].

Our results show that the reconstruction algorithm works well and incorporating the MLP in the backprojection improves spatial resolution by almost a factor of three compared to straight line projection in the spiral phantom used in this work (figure 8). Spatial resolution was comparable to that in energy-loss proton CT (figure 8). The rationale of introducing relative scattering power was to obtain a quantity which does not strongly vary as a function of angular variance. Our analysis based on tissue-like substitute materials showed that this is approximately true for some materials similar to or less dense than water. For denser bone-like materials, the relative scattering power decreases with increasing angular variance by up to 20% in a 20 cm thick phantom at 200 MeV and up to about 10% at 300 MeV (figure 1). In objects denser than water, such as an aluminium cylinder, this leads to a cupping artefact similar to beam hardening in conventional x-ray CT (figure 2). When dense objects are inserted into an otherwise water-like phantom, the degree to which the relative scattering power underestimates the material specific ratio X 0,w /X 0 of radiation lengths depends on the depth of the inserts (figure 3).

In Monte Carlo simulated images of the Gammex phantom, the reconstructed relative scattering power in dense inserts was systematically lower than the characteristic radiation length ratio X 0,w /X 0 and systematically higher in adipose, breast, and brain inserts (figure 7), consistent with our theoretical analysis. The image of the spiral phantom showed that the reconstructed relative scattering power underestimates the ratio X 0,w /X 0 (4.06 for aluminium) even more so the deeper the inserts lie in the phantom. At the centre, the reconstructed value was about 20% lower than the theoretical X 0,w /X 0 (see upper right panel in figure 8).

One way to correct for the under estimation in dense material would be to exploit its similarity to beam hardening and adapt strategies typically used in x-ray CT. Another option is to rely on iterative reconstruction and explicitly model the energy-dependence of MCS based on some a priori knowledge, as proposed by [START_REF] Bopp | Quantitative proton imaging from multiple physics processes: a proof of concept[END_REF]. However, the MLP would need to be incorporated into the forward projection matrix to achieve best possible spatial resolution, which is probably more complicated in scattering proton CT than in energy-loss proton CT because the angular variance needs to be estimated from several trajectories at some point during the reconstruction.

We investigated the statistical properties of scattering proton CT and derived an expression for the intrinsic relative scattering power precision (section 2.5). At equal dose, the noise in the reconstructed scattering proton CT images was about 50 times higher than in the energy-loss image, but about 10 times lower than in the nuclear attenuation image. As a consequence, the SNR and CNR was 40-50 times lower in the reconstructed scattering image of the Gammex phantom than in the energy-loss image. The interesting conclusion here is: scattering proton CT promises to yield higher contrast at first sight because the ratio X 0,w /X 0 is larger for some dense materials than their RSP. This has been pointed out e.g. by [START_REF] Plautz | 200 MeV Proton Radiography Studies With a Hand Phantom Using a Prototype Proton CT Scanner[END_REF]. The true limiting factor in scattering proton CT, however, is the higher image noise so that SNR and CNR turn out to be less favourable. The practical advantage of a scattering proton CT scanner is that it does not require any device to measure proton residual energy, such as a calorimeter or range telescope. This might be reason enough in some circumstances to accept the disadvantages compared to energy-loss proton CT.

Helium or even heavier ions have been investigated as alternative particles for ion imaging [START_REF] Parodi | Heavy ion radiography and tomography[END_REF], with the former hypothesised to be optimal in terms of spatial resolution of the images [START_REF] Collins-Fekete | A theoretical framework to predict the most likely ion path in particle imaging[END_REF][START_REF] Piersimoni | Helium CT: Monte Carlo simulation results for an ideal source and detector with comparison to proton CT[END_REF]. Studying these particles in the context of scattering proton CT in detail was beyond the scope of our work, but we share a few considerations here. The reason for potentially better spatial resolution in heavier ion CT is that MLP estimation is less uncertain for the heavier ions than for protons because they suffer less from MCS. If this is true for energy-loss CT, it will also be valid for scattering CT because the distance driven binning method is the same in both cases. Regarding noise properties for helium ions as an example, we refer to equation 26 and note that the prefactor containing the ratio of variance and scattering power is about 30% smaller for helium ions of initially 800 MeV than for protons of initially 200 MeV (equal range) behind a 20 cm water cylinder. At the same time, the stopping power at the cylinder centre (second factor in equation 26) is about a factor 4 larger for helium. The dose contribution due to nuclear interactions would require more detailed considerations, but will not change the order of magnitude of the overall expression. Therefore, we do not expect helium scattering CT to perform significantly better or worse than proton scattering CT in terms of image noise. It is worth mentioning in this context that electrons have also been explored for scattering CT imaging [START_REF] Jansen | Feasibility of track-based multiple scattering tomography[END_REF], but due to their much different mass and energy (order of GeV) we cannot immediately transfer our conclusions to that case.

A limitation of our Monte Carlo simulations is that we applied perfect particle filtering by removing all protons which underwent nuclear interactions in the phantom. In practice, such events can only be filtered out to some extent, typically by applying cuts on the angular and energy distributions [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. As nuclear scattering is associated with large angle scattering, it is likely that this will lead to some sort of systematic misestimation of the angular variance and thus of the reconstructed relative scattering power. Investigating these aspects was beyond the scope of this work.

Furthermore, our analysis of image noise concentrated on the intrinsic statistical limitations of the modality and disregarded experimental uncertainties, in particular the tracking devices. Their impact can be estimated as follows: According to equation 26, image noise is approximately (within the Gaussian approximation of MCS) proportional to the variance of the deviation angles ∆θ. Compared to an ideal simulation, uncertainties in the measured entrance and exit angle widen the distribution of deviation angles ∆θ and increase the effectively measured variance by σ 2 exp ≡ σ 2 entry + σ 2 exit , where σ entry and σ exit are the standard errors which quantify the tracker uncertainty. Image noise will increase proportionally. Measurement uncertainties will be important mainly in thin objects where the intrinsic angular variance due to MCS is small. The expected uncertainty of the angle measurement can be estimated based on the tracker geometry, material and the proton energy (Krah et al., 2018a). Using the Phase-II proton CT scanner prototype reported in [START_REF] Sadrozinski | Development of a Head Scanner for Proton CT[END_REF] as example, we calculated that noise in the scattering image of a water cylinder recorded at 200 MeV would increase by 20% for a cylinder diameter of 2 cm, but by less than 3% for a diameter of 10 cm or more. It is likely that heterogeneous phantom geometry and composition will contribute to image noise, as has been shown for energy-loss proton CT [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF]. The initial beam divergence does not impact image noise because the reconstruction uses the differences between exit and entry angles. Overall, the intrinsic image noise estimated in this work should be reasonably representative of the noise in a real scattering proton CT image. In the energy-loss proton CT image of a 10-20 cm thick phantom, on the other hand, contributions to image noise due to the beam's initial energy spread and the detection process are together of a similar magnitude as the intrinsic noise due to energy straggling considered in this work [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF].

The main reason for the development of energy-loss proton CT is that it can directly provide an RSP map of a patient needed by proton therapy treatment planning systems (TPS). Scattering proton CT as described in this work does not provide this quantity, but rather a map of the material's (relative) radiation length (X 0,w /X 0 ). At the same time, Monte Carlo TPS need more information about the tissue properties as input to the underlying physics models than RSP, e.g. elemental composition in Geant4 [START_REF] Agostinelli | GEANT4 -A simulation toolkit[END_REF] or directly radiation length in Fred [START_REF] Schiavi | Fred: A GPU-accelerated fast-Monte Carlo code for rapid treatment plan recalculation in ion beam therapy[END_REF]. These additional quantities can be obtained through tabulated conversion from another quantity [START_REF] Kanematsu | Relationship between electron density and effective densities of body tissues for stopping, scattering, and nuclear interactions of proton and ion beams[END_REF]. Alternatively, scattering proton CT could be combined with energy-loss proton CT, or even with x-ray based CT [START_REF] Vilches-Freixas | Deriving the mean excitation energy map from dual-energy and proton computed tomography[END_REF], to provide complementary information to the TPS. The accuracy of such a combined direct measurement would have to be at least as good as the one associated with a modelbased conversion. Investigating implementation details and whether the inferior image quality in scattering proton CT precludes its use in a clinical context was not the scope of this work.

We finally comment on the operating mode of the proton CT scanner: For both contrast mechanisms, operating a proton CT scanner in list-mode is demanding because the detection hardware needs to be extremely fast to keep up with the particle fluence of a therapeutic accelerator, even if operated at very low intensity [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF]. As a technically simpler alternative, "integrated-mode" approaches have been studied for energy-loss proton CT and among them those where the imaged object is irradiated with pencil beams [START_REF] Rescigno | A pencil beam approach to proton computed tomography[END_REF][START_REF] Farace | Pencil beam proton radiography using a multilayer ionization chamber[END_REF][START_REF] Bentefour | Concept of proton radiography using energy resolved dose measurement[END_REF]Krah et al., 2018b) and the detector records the mean energy loss or residual range per pencil beam. In this sense, the averaging operation occurs physically as part of the detection process while it is done numerically during distance driven binning in list-mode operation. Indeed, [START_REF] Meyer | Comparative Monte Carlo study on the performance of integrationand list-mode detector configurations for carbon ion computed tomography[END_REF] found comparable RSP accuracy with both operation modes for carbon ion CT, but we would expect a similar assessment for protons. Integrated mode scattering proton CT would require a detection device capable of measuring the angular variance (or another suitable quantity) per pencil beam directly instead of estimating it numerically during distance driven binning in list-mode. Whether new detector technology, e.g. high density scintillating glass [START_REF] Wilkinson | High density scintillating glass proton imaging detector[END_REF], might be suitable for this purpose and what accuracy is to be expected would require a dedicated study. Regardless of the contrast mechanism, integrated-mode scanners provide lower spatial resolution than list-mode scanners because reconstruction can only be performed along the average MLP of all particles in one beam and not each proton's MLP individually (Krah et al., 2018a).

Conclusions

In this work, we investigated a contrast mechanism of proton CT imaging which exploits the effect of multiple Coulomb scattering of protons in the imaged object. We showed that analytical reconstruction is possible by introducing a relative scattering power and developed a distancedriven filtered backprojection reconstruction, based on an algorithm previously proposed for energy-loss proton CT, which makes use of the most likely path concept. We found that the properties of the relative scattering power can lead to an underestimation of relative scattering power in image regions denser than water and to artefacts which are somewhat similar to beam hardening in x-ray CT imaging. We showed that spatial resolution in scattering proton CT was similar to energy-loss proton CT. Intrinsic noise was observed to be substantially higher than in energy-loss proton CT (by a factor of 40-50 at the centre of a 20 cm water cylinder). Our analysis revealed that signal-to-noise and contrast-to-noise is likewise worse than in energy-loss by at least a factor of 40 despite the higher contrast in dense tissues, mainly because of the higher noise. The practical advantage of scattering over energy-loss proton CT is that it only needs tracker devices, but no energy sensitive detectors such as calorimeters. It could therefore be interesting as a complementary modality in certain settings despite its inferior performance.

Figure 1 :

 1 Figure1: Relationships between angular variance, scattering power, and relative scattering power for 200 MeV (upper group of panels) and 300 MeV (lower panels). Coloured lines refer to Gammex insert materials; the black dash-dotted line always refers to water. Colored dashed lines in panels d and h indicate the radiation length ratio X 0,w /X 0 . All plots were generated analytically as explained in section 2.4.

Figure 2 :

 2 Figure 2: Left: Reconstructed relative scattering power image of the aluminium cylinder phantom. Right: Profile of relative scattering power as a function of distance from the cylinder centre. The interesting detail is the cupping artefact, i.e. the decrease of reconstructed values towards the centre.

Figure 3 :

 3 Figure 3: Relative scattering power as a function of depth in a 20 cm thick water phantom into which a bone-like slab of 2 cm is inserted at different depths. Beam energies: 200 MeV (left) and 300 MeV (right).

Figure 4 :

 4 Figure4: Left: Relative scattering/stopping power precision at the center of a water cylinder as a function of its diameter for scattering and energy-loss proton CT, respectively. Beware of separate y-axes. Right: Precision in a cylinder of 20 cm diameter as a function of dose at the cylinder centre for scattering (equation 26), energy-loss (equation 27), and attenuation proton CT (calculated as in[START_REF] Quiñones | Filtered back-projection reconstruction for attenuation proton CT along most likely paths[END_REF]).

Figure 5 :Figure 6 :

 56 Figure 5: Upper and middle panels: CNR and SNR of the Gammex phantom inserts obtained via Monte Carlo simulation. Data are presented as a function of the insert density. Lower panel: Material specific radiation length ratio X 0,w /X 0 and RSP.

Figure 7 :

 7 Figure 7: Mean relative scattering power δ values of the Gammex phantom inserts extracted from the reconstructed image shown in the centre panel in figure 6 (red dots). Red shaded areas indicate the standard error of the mean value. Black lines indicate the material specific radiation length ratio X 0,w /X 0 .

Figure 8 :

 8 Figure 8: Left column: Reconstructed images of the spiral phantom: scattering proton CT in straight line and MLP reconstruction (upper and center panel, respectively) and energy-loss proton CT in MLP reconstruction. Right column: Mean reconstructed value as a function of insert distance from centre (upper panel), example radial profile across an insert edge (centre panel), and spatial resolution as a function of distance from the phantom centre (lower panel).

Acknowledgements

The work of Nils Krah was partially supported by the project DIC20161236452 of the Fondation pour la Recherche Médicale and partially by funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 753370. This work was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). We also gratefully acknowledge the support from the CNRS/IN2P3 Computing Center (Lyon/Villeurbanne -France), for providing a significant amount of the computing resources needed for this work.


Edouard Pineau

Sébastien Razakarivony

Thomas Bonald

† Safran

Unsupervised ageing detection of mechanical systems on a causality graph

Multivariate time series (MTS) have specific features that complicate their analysis: interactions in space and time between the MTS components, variable length, absence of trivial alignment between samples and high dimensionality. Hence, finding a representation of MTS from which we can extract meaningful information is a challenging task. In general, specific assumptions are needed to obtain a valuable representation. In this paper, we assume that a dataset of MTS samples has an underlying causal structure that we can exploit to represent samples. Our contribution is a new representation framework that consists of first finding the overall causality graph G in a studied dataset and then mapping each sample onto G to obtain a causality-based representation. Since causality is an important feature underlying MTS data, we claim and show that representating samples on G is meaningful. We name this method Sequence-to-Graph (Seq2Graph). We apply Seq2Graph on health monitoring tasks, using two MTS datasets coming from ageing mechanical systems.

I. INTRODUCTION

Nowadays, more and more data is packaged as multivariate time series (MTS), for example industrial records, physiological data, vehicles sensors, etc. from which we may want to extract information in order to take decision from data. An important preliminary step for data information mining involves finding consistent and meaningful features that represent the observed samples. For MTS data, numerous features exist, for example moments or patterns. Each gives a certain amount of information about the data. Depending on the objective of the representation, certain features will be more relevant than others. In this paper, we propose to represent MTS with a particular feature: the causality.

The standard definition of causality used for time series data is the Granger causality [START_REF] Granger | Investigating causal relations by econometric models and cross-spectral methods[END_REF]. It consists in evaluating, for each couple of variables (X (i) , X (j) ) of a MTS sample X, if the variable X (i) (the cause) is useful to forecast the variable X (j) (the consequence). Hence, Granger causality is a feature that describes the causal dependencies underlying data. Causality can be represented as a graph G, as in Figure 1.

When observed samples are generated from a unique mechanical system S, it is common to assume a unique causality structure G for the whole dataset X . For example, G can be the skeleton of an articulated body or represent the causal relationships (statistical or physical) between sensors arranged in an engine, from which we observe samples. An illustration Corresponding author: edouard.pineau@telecom-paris.fr Figure 1: Graphical representation of the causality underlying a bivariate time series X. The edges are weighted with a set of edge's weights W X that depends on X. In this paper, we claim that comparing MTS samples by comparing their causality weights is relevant. is given in Figure 2. The causal structure G is then shared between all observed samples. When the assumption of a unique graphical structure for all samples is valid, like in the aforementioned examples, a causality graph G can be extracted directly from the studied dataset X , using an appropriate statistical method. G is then an abstract representation of the system S (e.g. a mechanical system) that generated all the observed samples.

Yet, Granger causality inference only aims at discovering and interpreting causalities between observed variables at the scale of the dataset (i.e. one graph for all samples). In this paper, we propose to relegate the inference of the causal graph G as a preliminary task preceding the representation inference of individual MTS samples. Our assumption is the following: since G is the graphical model of the system S that generated all the observed samples (e.g. a mechanical system), G is also a natural meaningful latent structure on which each data sample can be represented. In particular, each sample X can have its own set of edge's weights W X on G, hence its own causality-based representation. We can then compare samples by comparing their causality graph edge's weights. We claim and show in this paper that such a representation is meaningful for mechanical system's state monitoring.

We divide the problem into two subproblems. First, the extraction of the causal structure G using sparse MTS modeling approach [START_REF] Tank | Neural granger causality for nonlinear time series[END_REF]. Second, the choice and training of a relational neural network F φ parametrized by φ that takes MTS samples as input and maps them onto G [START_REF] Pineau | Seq2var: multivariate time series representation with relational neural networks and linear autoregressive model[END_REF] to obtain a causalitybased representation. Our main contribution is an unsupervised multivariate time series representation framework based on Granger causality, its implementation with neural networks Figure 2: Under the assumption that a unique causal graph G represents the causal structure underlying a mechanical system (here a plane), we claim that, if the causality assumption is relevant, the respective set of edge's weights W X , W Y and W Z associated to the MTS samples X, Y and Z can characterize the state of the system. Here, following the color code, W Z should be closer from W X than from W Y . We note that a maintenance was effectuated between Y and Z to restore the system's state. and its application to health monitoring. We name it Sequenceto-Graph (Seq2Graph). Section II introduces the problem. Section III details the learning procedure. Section IV gives related work. Section V illustrates the interest of Seq2Graph framework with two experiments of unsupervised mechanical systems health monitoring.

II. SEQ2GRAPH

In this section, we give the main steps of Seq2Graph. Training and implementation details are given in the next section.

Let X ⊆ R d×T be a finite set of d-dimensional multivariate time series (MTS), each indexed over the discrete time range t = 1, . . . , T . We want to represent each sample X ∈ X on a causality graph G.

A. Problem formulation in linear setting

We first assume that for each time series X, there exists a tensor W X ∈ R K×d×d such that

X t = K k=1 W X k X t-k (1)
Model ( 1) is a K order linear vector autoregressive model (L-VAR), with sample-wise parameter W X . We consider that W X represents time series X under linear assumption [START_REF] Ancona | Radial basis function approach to nonlinear granger causality of time series[END_REF]. We therefore can compare different samples X by comparing their parameter W X . In particular, the Granger causality is easily defined in VAR model (1): W X .,i,j = 0 means the absence of causality from variable X (i) to variable X (j) [START_REF] Eichler | Causal reasoning in graphical time series models[END_REF]. The non-zeros entries of W X are the weights of the causalities underlying sample X.

Each parameter W X can be estimated by maximum likelihood (MLE), in closed form or with a generic likelihood term (for example using a Kalman filter [START_REF] Casals | The likelihood of models with varying parameters[END_REF]). Yet, when the dataset X is large, when the samples X are high-dimensional or when W X has consistency or sparsity constraints, estimating every individual W X by MLE can be expensive or untractable. As a tractable alternative, we use [START_REF] Pineau | Seq2var: multivariate time series representation with relational neural networks and linear autoregressive model[END_REF]: a relational neural network (RelNN) [START_REF] Kipf | Neural relational inference for interacting systems[END_REF] F φ with parameters φ is trained to transform time series samples X directly into parameter W X , i.e. F φ (X) = W X . Hence, the representation problem becomes an encoder-decoder representation learning scheme, with F φ the encoder and the VAR model (1) the decoder.

B. Representation on a common causality graph G

As mentioned above, in VAR model (1), W X .,i,j = 0 means the absence of causality from variable X (i) to variable X (j) . Yet, we want to insure that all samples X are represented on the same causality graph G, i.e. all W X should share the same zeros.

In consequence, we assume that each tensor W X has three underlying components: a dataset-level component W ∈ R K×d×d , a sample-level component P X and dataset-level binary adjacency A ∈ {0, 1} d×d , such that:

W X = Ā W + P X (2)
where is the Hadamard product and Ā the adjacency A extended to match the dimensionality of W . W and A are shared between all samples X. The entries of the sparse tensor Ā W are the edge features of the causal graph G. The entries of the tensor Ā W X are the adjustment of the graph to match the properties of sample X.

We note that the decomposition of a regression problem into sample-level P X and dataset-level W is known as random coefficient regression (RCR) [START_REF] Bo Muthén | Random coefficient regression[END_REF]. The particularity of our approach is the shared sparsity given by adjacency A.

In practice, we first infer a sparse tensor W from the whole dataset X . Then we define a graph adjacency A from W such that ∀(i, j) ∈ 1, d :

A i,j = 1 { K k=1 | Wk,i,j |>0} . (3)
We then build and train a neural network P φ with parameters φ that directly and efficiently outputs the adjustment P X in the RCR. Hence, the previously introduced neural network F φ is defined as F φ (X) := W + P G φ (X), where P G φ (X) = Ā P φ (X). The details of the learning procedure and the implementation are given in Section III.

The Seq2Graph framework is illustrated in Figure 3.

C. Generalized L-VAR

To obtain a richer representation learning framework, we extend linear approach to a generalized linear [START_REF] Ashworth | Generalized linear models[END_REF] vector autoregressive model (GL-VAR) [START_REF] Tank | Neural granger causality for nonlinear time series[END_REF]. We remind that generalized linear models extend linear regression by relating the linear transform and the response variable with a nonlinear G is built such that it explains the mean dynamical behavior of the dataset. 2 is the representation inference. X is the prediction of the sample X. The adjustments P G φ (X) are built such that they explain the specific causality underlying sample X, along the edge of the mean parameter W . Here, K = 2. link function. Hence, let g = {g θj } d j=1 be a set of shallow neural networks with parameters θ = {θ j } d j=1 , g θj : R l → R, such that for each sample X ∈ X we can find a tensor

W X ∈ R K×d×d×l such that ∀j ∈ 1, d : X (j) t = g θj K k=1 W X k,.,j X t-k (4)
with l ∈ N * . ( 4) is a neural generalized linear version of a vector autoregressive (VAR) model, called generalized-linear VAR (GL-VAR) [START_REF] Tank | Neural granger causality for nonlinear time series[END_REF]. We remark the the nonlinearities are shared between all samples, like the common causal structure G.

D. Why is causality-based representation relevant?

The causality graph G is a main parameter of the model describing the system's dynamics (Eq. ( 1) or ( 4)). Hence, it is meaningful since it describes an important part of the physics underlying the system. Comparing two samples X and Y by comparing their respective weights W X and W Y on the graph G is therefore a substitute to the comparison of their physical hidden state.

We note that we can compare time series X and Y by comparing their respective set of causal weights W X and W Y with any relevant distance, for example the

l 2 distance W X -W Y 2 .

III. SEQ2GRAPH TRAINING

A. Inference of G with group-lasso training

We search a sparse W from RCR model ( 2) and θ = {θ j } d j=1 the set of parameters of the neural link functions (in the nonlinear case ( 4)). We therefore search a solution to the following mean-squared regression with group-lasso regularization:

min W ,θ E X∼X   d j=1 T t=K+1 X (j) t -g θj K k=1 Wk,.,j X t-k 2 2   + λ d i,j=1 W.,i,j F + γ d j=1 θ j (5)
where . F is the Frobenius norm. The group-lasso penalty for W associated with coefficient λ encourages each W.,i,j F to be null, meaning that all causal links from X (i) to X (j) would be cut. λ ∈ R + controls the speed and intensity of the pruning. Regularization θ j compensates the effect of the group-lasso to avoid the (theoretical) situation where W goes to zero while parameters θ tend to infinite sensitivity. Like in the linear case, we estimate causality adjacency A with

A i,j = 1{ K k=1 Wk,i,j > 0}.

B. Causality adjustment learning

In this section, we train a neural network that directly infers the causality adjustment of RCR model for a given sample.

We assume that the parameters {g, W } has been learned on the whole dataset (see Section III-A) and are fixed. We now build and train a relational neural network P φ to infer each sample-wise adjustment. We note P G φ = Ā P φ the adjustment constrained to live only on the edges of the inferred causality graph G. Ā ∈ R K×d×d×l is the adjacency matrix A of G expanded to match the dimensions of tensor P G φ . Then the problem to solve is:

min φ E X∼X d j=1 T t=K+1 X (j) t -g θj K k=1 W + P G φ (X) k,.,j X t-k 2 2 + ηΩ P G φ (X) (6)
Such learned inference network is therefore both meaningful (since built to represent the dynamics underlying the data) and sparse (since constrained on G). η is a parameter controlling the intensity of the penalty function Ω. A standard penalization Ω would be the l 2 norm on the parameters φ, under the assumption that the parameters of P G (X) are normally distributed around 0. A more general penalty (yet equivalent when minimized), that we use for Seq2Graph training is Ω(P G (X), P reg ) = P G (X) -P reg 2 + P reg 2 . P reg is a parameter trained during the optimization [START_REF] Eichler | Causal reasoning in graphical time series models[END_REF]. We have found that such penalty helps to obtain more consistent representation.

C. Implementation details a) Neural network representation inference: For the representation inference function P G φ , we use a relational neural network (RelNN) [START_REF] Santoro | A simple neural network module for relational reasoning[END_REF]. The RelNN embeds pair of variables. Its adaptation for time series data is taken from [START_REF] Kipf | Neural relational inference for interacting systems[END_REF], where the RelNN takes a MTS as input and embeds all pairs of variables into a binary space. In [START_REF] Pineau | Seq2var: multivariate time series representation with relational neural networks and linear autoregressive model[END_REF], an equivalent RelNN embeds pairs of variables as vectors specialized for linear causality, where it proved to be expressive, noise resistant and able to generalize over the notion of linear causality.

b) Sparsity inducing training for the dataset-level causality-graph inference: Problem ( 5) is first optimized with stochastic gradient descent. Then, we chose proximal gradient descent (PGD) [START_REF] Parikh | Proximal algorithms[END_REF] as fine-tuning optimization procedure in order to achieve true zeros in W as in [START_REF] Tank | Neural granger causality for nonlinear time series[END_REF]. If we had a target sparsity, we could stop the PGD when the level is achieved. In our experiment, we do not have (or assume to not have) the true sparsity level. We propose to monitor the impact of the sparsity on prediction and chose the maximal sparsity that does not degrade the prediction capacity of the model. See experiments for illustration. c) Multi-multivariate time series: There are cases where the d components of a MTS are multidimensional. For example, if the MTS has d variables situated in a 2D space (see Experiment V-A), hence each variable is represented by a 2D time series. More formally, the problem extends from X ⊆ R d×T to X ⊆ R d×m×T , i.e. ∀X ∈ X X t ∈ R d×m , with m the dimension of individual time series variables. The approach presented in our paper adapts to this general case by replacing W X ∈ R K×l×d by W X ∈ R K×l×d×1 . The additional dimension in W X enables to consider the m time series of each variable as a whole.

IV. RELATED WORK

This section outlines main related work for the two key concepts of our paper: MTS representation and causality inference.

A. Time series representation

Finding interesting and relevant features from time series data has a long-range history. A particularity of time series is that they have no explicit and general features [START_REF] Xing | A brief survey on sequence classification[END_REF]. Hence, unsupervised MTS representation requires strong assumptions to be relevant.

The most common assumption is that similar samples have similar shapes and patterns up to a warping alignment [START_REF] Esling | Time-series data mining[END_REF] and can be represented closely using bag-of-patterns methods [START_REF] Senin | Sax-vsm: Interpretable time series classification using sax and vector space model[END_REF]. This simple assumption gives good results for many time series tasks as long as the sought representation is contained in the observed patterns [START_REF] Wu | Random warping series: A random features method for timeseries embedding[END_REF]. Yet, for multivariate time series, it is limited since it generally does not take into account the interactions between variables. In [START_REF] Schäfer | Multivariate time series classification with weasel+ muse[END_REF], they propose to had the interplay of features in different dimension of MTS to obtain a more relevant bag-of-patterns representation.

More recently, the usage representation learning methods based on neural networks have emerged, giving new representation learning models for MTS data. In [START_REF] Lei | Similarity preserving representation learning for time series analysis[END_REF], they propose a method such that the distance between the learned representations is the dynamic time warping (DTW) distance. A powerful family of neural representation learning methods is the autoencoders (AE) and related methods [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF]. The principle is to train an encoder and a decoder simultaneously, such that the encoded signal of each sample can be correctly decoded. If so, it means that the encoded signal contains the essential information. The Sequential Autoencoders (SAE) [START_REF] Cho | Learning phrase representations using rnn encoder-decoder for statistical machine translation[END_REF], [START_REF] Malhotra | Timenet: Pre-trained deep recurrent neural network for time series classification[END_REF], [START_REF] Zhang | Sequence-to-sequence acoustic modeling for voice conversion[END_REF] is the most popular adaptation of AE for time series data. A richer SAE proposed in [START_REF] Vincent Fortuin | Som-vae: Interpretable discrete representation learning on time series[END_REF] is based on the joint learning of a discrete variational autoencoder (VAE) [START_REF] Van Den Oord | Neural discrete representation learning[END_REF], a self-organizing map (SOM) latent space [START_REF] Kohonen | Self-organized formation of topologically correct feature maps[END_REF] and a Markov transition model [START_REF] Charles | Practical markov chain monte carlo[END_REF]. Yet, the latter model only treats discrete representation of time series. In [START_REF] Franceschi | Unsupervised scalable representation learning for multivariate time series[END_REF], an unsupervised scalable time series representation (USTR) is proposed using the notion of triplet loss. The assumption is that a MTS sample is closer to one of its subsamples (positive sampling) than to a randomly chosen sample of the dataset (negative sampling). This model achieves very good results in classification downstream task. Moreover, it can handle many types of data since the assumptions underlying the model are weak.

Closer to our work, [START_REF] Kipf | Neural relational inference for interacting systems[END_REF] proposes the neural relational inference (NRI), a VAE that transforms time series into a binary relational graph, trained as a variable selection method for neural time series model. Despite high interest and impressive results, NRI is limited to binary representations. [START_REF] Pineau | Seq2var: multivariate time series representation with relational neural networks and linear autoregressive model[END_REF] extends the binary relationships of NRI to a linear Granger causality latent space, with a model called Seq2VAR. Their model is applied to causality detection in contexts where NRI fails (noisy environment, floating intensity of the causalities). Our approach is built upon Seq2VAR, by adding non-linear causality and dataset-level regularization.

In the current paper, we use SAE, USTR and Seq2VAR as comparative models in the experiment part. The patternbased representation are not used since in both experiments we search latent (hidden) signal that is not present the observed patterns.

B. Granger causality inference from time series data

The global Granger causality problem is a statistical test for all pairs of variables where null hypothesis is "absence of causality". In its most common form, Granger causality (GC) is defined in linear models [START_REF] Ancona | Radial basis function approach to nonlinear granger causality of time series[END_REF], with the non zero entries of the parameter of a sparse linear vector autoregressive (L-VAR) model [START_REF] Richard A Davis | Sparse vector autoregressive modeling[END_REF], [START_REF] Han | Transition matrix estimation in high dimensional time series[END_REF], obtained with a penalization on L-VAR parameters during training.

For the nonlinear case, finding a useful latent representation is far from obvious. Among standard approaches for nonlinear GC detection in time series, we find kernel methods [START_REF] Ancona | Radial basis function approach to nonlinear granger causality of time series[END_REF] or information-theory [START_REF] Solo | On causality and mutual information[END_REF]. These methods find the existence of causality between variables, i.e., the causality graph. Yet, these methods do not give information on the importance of the causality where it exists. In Seq2Graph, as shown in the following, we require continuous (real) features on the edges of the causality graph. In [START_REF] Yang | Detecting nonlinear causality in multivariate time series with sparse additive models[END_REF], the authors extend VAR to nonlinear causality with the linear combination of B-spline basis functions. The linear operator, once fitted with appropriate penalization, contains continuous causal information. Another nonlinear extension of L-VAR causality is presented in [START_REF] Tank | Neural granger causality for nonlinear time series[END_REF]. The authors build the nonlinearity with neural networks. They find causality by applying a penalization on the input convolutional kernel. We use this last work for our causalitybased dataset-level regularization.

V. EXPERIMENTS

We propose two experiments to illustrate the interest of Seq2Graph. Both are based on multivariate time series data generated from mechanical systems. For the first, we use synthetic data with controlled causal structure to show the interest of the regularization with global causality graph G with a sparse random coefficient model [START_REF] Bengio | Representation learning: A review and new perspectives[END_REF]. For the second, we use a NASA dataset to assess Seq2Graph on real representation task. In all experiments, we compare Seq2Graph to three time series representation methods: a sequential autoencoder (SAE) [START_REF] Malhotra | Timenet: Pre-trained deep recurrent neural network for time series classification[END_REF], to the unsupervised scalable time series representation (USTR) [START_REF] Franceschi | Unsupervised scalable representation learning for multivariate time series[END_REF] and to the sequence-to-VAR (Seq2VAR) [START_REF] Pineau | Seq2var: multivariate time series representation with relational neural networks and linear autoregressive model[END_REF].

For all experiments, we use PyTorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF]. All models were trained with Adam optimizer [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] with learning rate 5.10 -4 . For the sample-level representation training, we added a learning rate scheduler with step size 100 and a decay factor = 0.9. All experiments computations on a GPU Nvidia Quadro K4000. The hyperparameters are given in Table I.

λ γ η Experiment V-A 10 -3 - 10 -5 Experiment V-B 5 × 10 -3 5 × 10 -3 5 × 10 -4
Table I: Hyperparameters for our experiments.

Code to reproduce the experiments can be found at https: //github.com/anonym-conf-submission/Seq2Graph.

A. Interacting Newtonian system 1) Dataset: We simulate samples from a 10-ball-springs system, consisting of the simultaneous trajectories of 10 identical balls of unit mass in a 2-dimensional space, each ball being connected to some others by springs (the rate of connection is 56%). This system has a natural bidirectional causal structure: each ball's trajectory acts as a cause for changes in the trajectory of the neighbor balls, and conversely. Using the previously introduced notations, we have d = 10 (10 balls) and m = 2 (in a 2-dimensional space, see implementation details). System dynamics follows Newton's law of motion. We assume that the system is ageing and is regularly restored. All samples share a common graph graphical structure G whose adjacency is the interaction matrix formed by the springs.

We simulate a synthetic dataset of 15000 samples (trajectories), 5000 for train, 5000 for validation and 5000 for test. Each trajectory is 49 time-steps-long (T = 49). For each batch b of 50 samples, a constant ageing factor α b ∼ U([0.9, 1]) is applied to the system: at each sample X whose index is s ∈ 0, 50 (within the batch b of 50 samples), we randomly choose a spring (i, j) and multiply its rigidity by α s b , i.e. an exponential ageing coefficient with respect to sample index. Every 50 samples, we restore the state of the system and another ageing factor is sampled and applied to the next batch of 50 samples. For some trajectories, α b = 1, i.e. there is no ageing: the initial hidden causality graph has binary adjacency and remains the same along the life of the system. When α b < 1, the initial graph is deteriorating during along the life (observed through 50 samples) of the system, until restoration. 2) Model: For this first experiment, we assume that the model is linear, i.e. that functions g θj are identities. The objective is to illustrate the impact of representing data on the same causality graph G. We determine by cross-validation that K = 2. The level of sparsity is determined by the quality of the prediction for different levels of sparsity of the causality graph, as shown in Figure 4. The prediction is almost invariant until a sparsity of about 56%. We note that we find back the true adjacency with [START_REF] Richard A Davis | Sparse vector autoregressive modeling[END_REF].

3) Metrics and results: We assess the quality of the representation inference function P G φ . We test if we can represent the ageing of the system with respect to a reference healthy sample X ref (first sample of a batch) picked in the validation set. We then build the test ageing curve

X -→ K k=1 (P G φ (X ref ) -P G φ (X)) k 2 2 (7)
for all samples X ∈ X test . This ageing curves are presented in Figure 5 and Table II We see that Seq2Graph outperforms both SAE, USTR and Seq2VAR for unsupervised representation learning, when meaningful information is fully contained in the causality. In particular, SAE and USTR completely miss the consistent ageing information, as expected from pattern-based approaches. Although consistent, Seq2VAR seems to suffer when the causalities become lower. In fact, lowering the causality improves the difficulty to capture them, hence prevents to find the trend hidden in causality. Adding a common causal structure W in Seq2Graph naturally helps the identification and the consistency of the sample representations.

B. NASA turbofan degradation simulation dataset 1) Dataset: NASA public Commercial Modular Aero-Propulsion System Simulation dataset (C-MAPSS) is a tool for simulation of realistic large commercial turbofan engine data [START_REF] Saxena | Turbofan engine degradation simulation data set[END_REF]. An engine degradation simulation was carried out using C-MAPSS, under different conditions and different faults. We use the FD001 dataset which contains 100 time series recorded at sea level with one fault mode for each (degradation of the high-pressure compressor, a fundamental turbofan piece). The time series are the output of the turbine-engine system that takes a fuel flow as input and outputs 21 variables, whose 13 are not constant (we only keep these 13 variables). Time series are 206 time-steps long on average. Each time series is the recording of a turbine engine going to failure. The engine is operating normally at the start of each time series and develops a fault of unknown initial magnitude in its first moments. We only know that the impact of this fault on the system grows in magnitude until system failure.

For the results of the paper, we split the dataset: the first 60 time series are train set, the 15 next are validation set and the last 25 are test set. We extract from these time series subtrajectories of length 25, with a rolling window with stride 5 to make our dataset. Hence, as for previous experiment, we have several batches of samples. A batch corresponds to the life of the system from start to failure. At the end of each batch, the engine is restored and another batch of samples is recorded.

2) Model: Using the previously introduced notations, d = 13. All samples share a common (unknown) structure which is the turbine engine mechanics. We assume that this structure can be represented by a sparse causality matrix. For this experiment each link function g θj is a MLP with 2 hidden layers of 4 channels. We determine by cross-validation that K = 2. The maximal sparsity is determined by the quality of the prediction for different levels of sparsity, as shown in Figure 6. The prediction is almost invariant until a sparsity of about 75%. We solve (6) using the previously found causal graph. Contrary to the previous experiment, the system is not isolated since the observed variables are the response to an unobserved command (the fuel flow). Yet, the variables can effectively interact with each other and statistical causality still makes sense as a representation assumption.

3) Metrics and results: We assess the quality of the learned P G φ by extracting ageing information from the representations, like in the previous experiment. Yet, we remind that, this time, we do not know the importance of the initial fault and we do not know the regularity of the ageing nor its importance (the ageing coefficient). What we know is that the 100 engines go to failure and the degradation of the state is monotonic until restoration. We propose a two-step process to validate Seq2Graph for NASA, by showing our ability to infer from data representation the imminence of a failure.

First, we build an ageing indicator assuming that is a relative position compare to a healthy sample. We pick a healthy sample X ref (first sample of a batch) in the validation set and build the ageing curve [START_REF] Esling | Time-series data mining[END_REF] for all X ∈ X valid . We compute a failure threshold τ valid that must indicate when an engine goes to failure. We set τ valid to the maximal threshold that ensures turbine engine failure detection for all validation batches, that is:

τ valid = min X∈X val,def K k=1 (P G φ (X ref ) -P G φ (X)) k 2 2 (8)
where X val,def is the set of validation samples preceding the engine failure. We note that τ valid has no safety margin, i.e. any threshold above τ valid misses at least one engine failure in the validation set. It is possible to add a margin by lowering τ valid . Second, we build the test ageing curve [START_REF] Esling | Time-series data mining[END_REF] with the samples of X test . We apply the detection test using τ valid (represented by the horizontal dotted line in Figure 7. As a first assessment, we see in Figure 7 that the estimated ageing curves built from SAE, Seq2VAR and Seq2Graph are almost monotonic inside each batch (between two vertical orange lines). We recall that monotony is the only ground truth information we have on the ageing of the system. The fact that SAE, Seq2VAR and Seq2Graph unveils monotonic signal means the ageing information is present both in patterns and values (SAE) and in causality (Seq2VAR and Seq2Graph). We do not find consistent representations with USTR. We also observe that the batch's ageing curves do not begin at the same value (dashed horizontal lines in Figure 7), whatever the method. It is partly imputed to the fact that the mechanical faults are located at the beginning of each batch and that they vary in intensity. Hence, the inferred first samples of each batch do not have to be equal.

We now compare the ability of the different MTS representations to relevantly detect failures. In Figure 8, we see the proportion of alarm at different time steps before actual failure happens, built from the estimated ageing curves illustrated in Figure 7. First, we note that all models detect almost 100% of failures before it happens. Second, we want detection of the coming failures to be reasonably early to avoid false alarms. If curves cross threshold too early, the MTS representation is useless. Figure 8 shows that Seq2Graph is the most consistent in early detection with no alarms far from failure, due to the consistency of the extracted monotonic signal. On the contrary, SAE always finds early failures, hence is less confident in failure detection.

We conducted two experiments on datasets samples from ageing mechanical systems. We have shown that using a common causality graph G is relevant to representation data samples. We show that we can find continuous signal hidden in data (ageing) from the causality based representation. These experiments validate our assumptions.

VI. CONCLUSION

In this paper, we have presented a multivariate time series (MTS) representation framework under the assumption that Granger causality contains relevant information about data. We have proposed a two-step approach, based on neural networks. First, the global causality graph G is found using the weights of a (generalized) linear autoregressive model trained with a group-Lasso penalization (to obtain sparsity). Second, a relational neural network is trained to efficiently infer the representation of each sample, constrained by the mean causality graph G. Our approach was motivated by health monitoring problems. We therefore illustrated Seq2Graph on two tasks of mechanical system health monitoring.

Figure 3 :

 3 Figure 3: Seq2Graph. 1 is the dataset-level causal graph G inference. X is the set of predictions for the whole dataset X .G is built such that it explains the mean dynamical behavior of the dataset. 2 is the representation inference. X is the prediction of the sample X. The adjustments P G φ (X) are built such that they explain the specific causality underlying sample X, along the edge of the mean parameter W . Here, K = 2.

Figure 4 :

 4 Figure 4: Top: Prediction MSE. Black line is the true sparsity. Bottom: Causal graphs for different sparsity levels. The third figure is the inferred causal graph, which matches the ground truth.

Figure 5 :

 5 Figure 5: Unsupervised estimation of ageing curves for the 10 first batches of the test set. From top to bottom: USTR [9], SAE [17], Seq2VAR [22], Seq2Graph. Orange curve is the ground truth.

Figure 6 :

 6 Figure 6: Top: MSE of the prediction. Black line is the sparsity of the inferred causal graph. Bottom: Causal graphs for different sparsity levels. Fourth image is the inferred causal graph.

Figure 7 :

 7 Figure 7: Unsupervised estimation of C-MAPSS ageing curve with different models on the 7 first test batches. From top to bottom: USTR [9], SAE [17], Seq2VAR [22], Seq2Graph. Orange picks are engine failures and repair. Long red dotted horizontal line is the threshold τ valid . Black dashed horizontal lines are the estimated initial states of each engine, computed as the mean value of the curve on the 10 first samples of each batch.

Figure 8 :

 8 Figure 8: Early alarm on CMAPSS data using MTS representation models USTR, SAE, Seq2VAR (see related work for details) and Seq2Graph. Means and standard deviations are built using all batch's first samples as X ref and several encoders trained with different seeds.

Table II :

 II . Performance of several models plus ours on ageing ball-springs problem. The Ageing score is the correlation between estimated and real ageing curve. *Higher is better. MSE stands for mean squared error and serves only as a sanity check (for MSE-based methods).

	Models	MSE	Ageing score*
	SAE	2.2 × 10 -5	0.09
	USTR	-	0.05
	Seq2VAR	2.3 × 10 -7	0.62
	Seq2Graph	4.4 × 10 -7	0.97

Acknowledgment This work is supported by the company Safran through the CIFRE convention 2017/1317.