Didier Bardel

Michel Perez

Daniel Nelias
email: daniel.nelias@insa-lyon.fr

Thibaut Chaise

Jer Ôme Garnier

Florent Bourlier

EVOLUTION OF MECHANICAL BEHAVIOR OF 6XXX ALUMINIUM ALLOY DUE TO THE PRECIPITATION STATE DURING A THERMO-MECHANICAL PROCESS

The aim of this research is to link the microstructural state and the mechanical properties of an age hardening alloy during a fast heat treatment such as encountered during welding.

A coupled model between precipitation state and mechanical properties is used to predict the yield strength and hardening behavior that can be observed experimentally. The method permits the identification of the kinematic and isotropic contributions in the hardening model. The methodology is applied to a 6061-T6 aluminium alloy which is used in the Jules Horowitz reactor vessel.

The general idea of this methodology is to couple an efficient microstructural model to a mechanical one based on the dislocation theory and ad'hoc experiments. The theoretical background is based on the work of Kampmann and Wagner, known as the KWN model, to account for nucleation, growth/dissolution and coarsening of precipitates. This analysis requires transient thermo-mechanical experimental data. The efficiency of these models and their coupling are shown for a serie 6XXX aluminium alloy which contains β and β precipitates. Ultimately these

models are coupled to a FEA model and allows to predict the distribution of precipitates within each element of the mesh, and subsequently its mechanical behavior.

NOMENCLATURE

INTRODUCTION

The work presented in this paper is a part of the study of early qualification of pressure vessel for a new French experimental reactor. The approximate size of this piece is five meters high with a diameter about one meter. To obtain this piece made of AA6061-T6, several shells are welded together using a welding process: the electron beam. The ultimate goal is to predict the metallurgical state generated by the welding process and to deduce the mechanical properties and residual stresses in the vessel.

To achieve this goal, this work tries to understand the link between microstructural evolution and mechanical properties. In a first part, the implementation of a microstructural model is presented. This implementation is based on the works of Hillert [START_REF] Hillert | Diffusioncontrolled lengthening of widmansttten plates[END_REF]. This model is validated using measurement of Small Angle Neutron Scattering (SANS) and Transmission Electron Microscopy (TEM). These experimental campaigns provide informations about the distribution of precipitates present in the alloy for several heat treatment representative of welding. For SANS, an ellipsoidal distribution is using to have a good fitting of raw diffusion spectra. The microstructural state of the alloy for rapid heat treatments can then be provided by finite element simulations. In a second part, a modelisation is presented to link the precipitation state and mechanical properties. This modelisation is validated with cyclic tensile tests performed on various precipitation states.

Microstructural modelisation

The KWN model allows to predict the different microstructural phenomena that can occur during various heating kinetics encountered in a welding process: nucleation, growth, and coarsening. This model is generally implemented for spherical precipitates ( [2] [START_REF] Simar | Sequential modeling of local precipitation, strength and strain hardening in friction stir welds of an aluminum alloy 6005a-t6[END_REF] ) but for the aluminium alloys this hypothesis is not acceptable because precipitates have often large aspect ratio. Indeed, for the 6000 series the Heat Affected Zone (HAZ) is composed by needle shaped precipitates (β and β ) [START_REF] Myhr | Modelling of the microstructure and strength evolution in al-mg-si alloys during multistage thermal processing[END_REF]. Here, Classical Nucleation and Growth Theories (CNGTs) have been adapted to model accurately the precipitation of needle shaped particles. Hillert proposed a modification of the Zener growth equation considering a constant tip radius but without coarsening treatment. Here, an adaptation of this work is presented to have a diffusion/coarsening stage. This adapted expression is implemented in a 'Lagrangian-like' model class management software (PreciSo) to obtain the distribution density of precipitates. Thus, the evolution of precipitate length can be written as follows for a constant aspect ratio [START_REF] Bardel | Evolution of precipitation state and strength for welding of a 6061 aluminium alloy[END_REF]:

           dl dt = 3D Mg 4R p X 0 Mg -X i Mg α m/p X p Mg -X i Mg dl dt = 3D Si 4R p X 0 Si -X i Si α m/p X p Si -X i Si X i Mg x X i Si y = K s .exp ξ ξ -2/3 4V p mol γ 3RkT (1)
Here, as in [START_REF] Myhr | Modelling of the age hardening behaviour of al-mg-si alloys[END_REF] a fictitious phase is proposed with β -β nature. This equivalent phase allows to obtain a good representation of the microstructural state for several thermal loads representative of welding. These simulations are compared with both TEM and SANS measurements applied on tensil test samples [START_REF] Bardel | Evolution of precipitation state and strength for welding of a 6061 aluminium alloy[END_REF].

Mechanical tests

Tensile tests have been carried out at room temperature on specimens previously heated up to peak temperatures of 473, 573, 673, 773 and 833 K and with various heating rates (0.5, 5, 15, 50, 200 K/s) in order to measure their mechanical properties. The thermal loadings are representative of the thermal histories encountered in Electron Beam welding [START_REF] Maisonnette | Influences mcaniques et mtallurgiques de procds haute temprature sur un alliage daluminium 6061-t6[END_REF]. The experimental device that is used here has been specifically developed at INSA Lyon to meet these kind of requirements [START_REF] Maisonnette | Effects of heat treatments on the microstructure and mechanical properties of a 6061 aluminium alloy[END_REF].

The precipitation state is then used as an input parameter of a mean field micro-mechanical models that predict the yield strength of the alloy.

σ y = σ 0 + ∑ k i .C i SS 2/3 + M b ∑ N i l i 2 √ 3Γ ∑ i N i F i N T OT 3/2 (2)
These results are successfully compared to tensile tests experiments performed on various heating conditions, representing the heat-affected zone of a weld. In this paper the influence of several assumptions on the yield stress modeling have been reviewed and new equations are used [START_REF] Bardel | Evolution of precipitation state and strength for welding of a 6061 aluminium alloy[END_REF]. The use of a full distribution approach allows to improve the physical meaning of the fitted parameters. Here, only one fitted parameter is used and the identified Taylor's factor is perfectly compatible with recent works on aluminium alloy [START_REF] Starink | A model for the yield strenght of overaged al-zn-mg-cu alloys[END_REF].

Then, for hardening curves, classical elastoplastics laws can be used such as Amstrong-Frederic models [START_REF] Jiang | Characteristics of the armstrong-frederick type plasticity models[END_REF]. In these laws, the yield stress is predicted by the previous predicted approch. Other coefficients values are dependant on the simulated precipitates mean radius and volume fraction.

CONCLUSION AND PERSPECTIVES

In this paper, a coupled physical precipitation/yield stress modelisation is presented. In the first part, an adjustment strategy is proposed for a single fictitious phase to fit microstructural evolution obtained by SANS. Then, yield stress prediction are validated with tensile tests. This modelisation provides a good fitting in a wide range of thermals histories and can be used in finite element software with phenomenological plasticty models such as Amstrong Frederic's.

The prospect of this work is the identification of viscous part of mechanical behaviour. Then, all the tools will be available for prediction of residual stresses for age hardening aluminium alloys.

ACKNOWLEDGMENT

The authors gratefully acknowledge Areva and the French Atomic Energy Commission (CEA) for industrial and financial support. The staff of the D22 beamline of Laue Langevin Institut is thanked for technical help and Alexis Deschamps for his SANS advices.

FIGURE 1 :

 1 FIGURE 1: SIMULATED VS EXPERIMENTAL MEAN RA-DIUS OF NEEDLES.

FIGURE 2 :

 2 FIGURE 2: SIMULATED VS EXPERIMENTAL VOLUME FRACTION OF NEEDLES.

FIGURE 3 :

 3 FIGURE 3: INSTRUMENTATION FOR TENSILE TESTS.

FIGURE 4 :

 4 FIGURE 4: YIELD STRESS PREDICTION VS EXPERI-MENTS.

Proceedings of the ASME 2013 Pressure Vessels and Piping Conference PVP2013 July 14-18, 2013, Paris, France PVP2013-97426

  Ratio of atomics volumes for precipitate and matrix.

	V p mol	Molecular volume of precipitate.
	α m/p σ y	Yield stress.
	σ 0	Yield stress of pure aluminium.
	k i	Element strengthening constant of element 'i'.
	C i
		l	Length of needle.
		D Mg Diffusion coefficient of magnesium.
		D Si	Diffusion coefficient of silicon.
		R p	Tip radius of needle.
		X 0 Mg	Atomic fraction of magnesium in matrix.
		X i Mg	Atomic fraction of magnesium at matrix/precipitate in-
			terface.
		X p Mg	Atomic fraction of magnesium in precipitate.
		X 0 Si X i Si X p Si	Atomic fraction of silicon in matrix. Atomic fraction of silicon at matrix/precipitate interface. Atomic fraction of silicon in precipitate.
		γ	Matrix/precipitate interface energy.
		K s	Solubility product of precipitate.
		ξ	Needle aspect ratio.
		k	Boltzmann constant.


M Bogdan

E Roubin

J.-B Colliat

F Benboudjema

L Stefan

MORPHOLOGICAL MODELING OF CEMENT BASED MATERIALS AND HYDRATION PROCESS

In this work we present a new morphological model for heterogeneous media that is well suited for cement based materials. The main application addressed here is a model for the evolution cement paste with respect to the hydration degree, containing only anhydrous cement, water, and hydrated products. Secondly, those morphologies are projected on FE trusses, which allow us to predict the evolution of any chosen material characteristic (elastic properties and diffusivity) with the evolution of the hydration degree.

1.

INTRODUCTION

Cement based materials modeling has always been a great challenge. For nuclear or construction industry, a reliable predictive approach is still an undergoing challenge. A new approach based on random field's thresholding is presented here, in order to generate microstructures of heterogeneous materials. This mutli-scale morphological modeling framework is particularly adapted to cement based materials since it is capable to represent easily a material with multiple concentrical phases.

Following the work of R.J. Adler [START_REF] Adler | Random Fields and Geometry[END_REF], the idea is to yield phases from random fields' excursion sets. The advantages of such a method are various. First, contrary to most of the numerical hydration models, we will not be dealing with spheres but random shapes of inclusions. Secondly, the creation of such morphologies is much lighter than average, since we will need only three parameters (two for the random field, and one for the threshold) for an entire morphology (instead of four parameters for each inclusion).

Eventually, one can generate various random fields, with chosen characteristics. By thresholding them, various morphology can be created, controlling parameters such as volume fraction or specific surface area, and thus create a heterogeneous model for any material.

The focus is made here on anhydrous cement within a "matrix" of water. By changing the thresholds, a n-concentrical-phases material can be created. A simple analytical hydration model [START_REF] Powers | Studies of the physical properties of hardened Portland cement paste[END_REF] has been chosen, considering only three phases: anhydrous cement, hydration products (including gel water) and free water. The link that has been made between these volume fractions and the thresholds allows us to make evolve the initial morphology, and thus model the hydration process for a cement paste.

Once the morphologies have been created, at any hydration degree, they are projected on a FE mesh with embedded discontinuities [START_REF] Benkemoun | Failure of heterogeneous materials: 3D meso-scale FR models with embended discontinuities[END_REF]. It gives us the possibility to predict properties such as the Young modulus, the Poisson's ratio, the effective diffusivity, or any transport properties, using a sequenced multi-scale approach.

MORPHOLOGY CREATION

In this part we explain shortly how the morphologies are created. With first the generation of random fields with chosen characteristics, then the thresholding, and finally the link between the random fields variables and the geometric characteristics of the excursion sets.

2.1

Random fields The first step in the framework we present here is the generation of random field realisations. We give here only the main ideas and present the main tools used for the generation process. Since the whole process relies on correlated random fields, the main issue was the numerical implementation of their generation.

First, we used the orthogonal decomposition [START_REF] Karhunen | Über lineare Methoden in der Wahrscheinlichkeitsrechnung[END_REF], which allows, for Gaussian correlated fields, a separation between the space and random variables. In a second time, the Karhunen-Loève decomposition [START_REF] Loève | Probability theory[END_REF] is used to determine the space variables by solving an eigenvalues problem. Demonstration can be found in [START_REF] Adler | Random Fields and Geometry[END_REF]. The numerical implementation relies here on a Finite Element Method. Finally, the last step is to generate a field with such space functions. For the random variables, a simple random number generator can be used. However, such multi-dimensional, large size, generation processes are very demanding in memory storage. In this regard, we used here the Turning Band Method [START_REF] Matheron | The intrinsic random functions and their application[END_REF], which allows creating mutlidimensional fields from several one-dimensional ones.

Excursion sets

The second step is the thresholding of those realisations. Hereafter, we explain the principle of excursion sets. We shall consider random fields as random functions: γ(x,ω) : M R N →R k , where M is a N dimensional Euclidian space. In the framework we present here, we are dealing with Gaussian (or Gaussian related) correlated random functions γ, which are defined in a cube of size T, with value in R. An excursion A u of such a function is defined by the points of M where the values of γ are above a threshold u (equation ( 1)):

A u ≡ A u (γ,M) = { x XM : γ(x) ≥ u } ( 1
)
An example is shown on Figure 1, for a one dimensional realisation γ, defined over M. A second example (Figure 2) shows the variety of excursion sets yielded for different values of the threshold u, from sponge (left side) to meatball like (right side), with all the intermediaries.

Characterization

The whole point of using this framework lies in the work of R.J. Adler [START_REF] Adler | Some new random field tools for spatial analysis[END_REF]. He established a probabilistic link between the random fields' quantities (L c and σ) and the threshold on one side, and the geometrical and topological quantities of the excursions sets on the other side.

To be able to describe the morphologies, one needs a measure. It has been chosen here to use the Lipschitz-Killing curvatures (denoted LKC's). In an N dimensional space, there are N+1 such curvatures. So here, for a three dimensional geometry, we have four LKC's. Three of them are geometrical, and the last is topological. They are defined by:

-

L 3 (A u ) is the volume of (A u ) -L 2 (A u ) is half the surface area of (A u ) -L 1 (A u ) is twice the calliper diameter of (A u ) -L 0 (A u )
is the Euler characteristic of (A u ), which is the only topological measure, that can be calculated by:

L 0 (A u ) = #{connected components in A u } -#{"handles" in A u } + #{"holes" in A u } (2)
Moreover, a positive Euler characteristic reflects a disconnected morphology, and a negative one reflects a connection between the components. The link previously mentioned, is so-called probabilistic since its only gives us expectations of the LKC's. However, since we have only three parameters as inputs (correlation length, variance and threshold), we choose to keep only three parameters as outputs: the volume fraction V v , the volumic surface area S and the number of aggregates N.

So in the end, we will have a link that can be expressed as following:

E{L 3 }(u,σ) = V v .T 3 E{L 2 }(u,σ,L c ) = ½.S. T 3 (3) E{L 0 }(u,σ,L c ) = N
where T is the size of the excursion. For further details about those equations, please refer to [START_REF] Adler | Random Fields and Geometry[END_REF] for the explicit formulae, and [START_REF] Adler | Some new random field tools for spatial analysis[END_REF] for the demonstrations.

The only issue here is that those equations are highly non-linear. So most of the time, the previous system of equation does not have solutions. Typically, so far we cannot yield excursions with more than 15% of volume fraction using Gaussian fields, and no more than 30% using chi-square distributions. It still isn't very realistic regarding concrete mixes used nowadays (it would represent cement pastes with w/c ratios higher than 0.7), but it is the best that we can do so far.

A last point worth mentioning, is that excursion sets are not necessarily defined by "all the points of M where values of γ are above a given threshold u". We can also define them by "all the points of M where values of γ are between two thresholds u and v", which will allow us to create several "concentrical" phases.

HYDRATION PROCESS

In this part we present an application to the cement paste hydration process, using a simple analytical model to create initial morphologies, and make them evolve.

Analytical model

The Powers and Brownyard model [START_REF] Powers | Studies of the physical properties of hardened Portland cement paste[END_REF] has been implemented here, but in a simple version. Usually, this model considers five volumic fractions: anhydrous cement, hydration products, gel water, free water and chemical shrinkage. We decided to use only three for this first application, taking into account only anhydrous cement, hydration products (including gel water) and free water . We also decided not to take into account the chemical shrinkage (its volume is taken into account within the volume of water). Considering these assumptions, we also decided that the hydration degree should be able to reach one (due to the high w/c ratio, this hypothesis is realistic) at the end of the process in order not to underestimate the volume of hydration products (20% of the initial volume of anhydrous cement is included in the pores at the final stage). The simplified equations can be resumed as follows:

V anh = (1-p)(1-α) V h = 2.12(1-p)α V w = 1 -V anh -V h
where p is the initial porosity, α the hydration degree, and V anh , V h , V w respectively the volume fractions of anhydrous cement, hydration products and water.

Initial morphology

The geometrical and physical description of the anhydrous cement is simple for a classical Portland cement, without any additions. If we consider a CEM-I cement, the usual average density is ρ c ≈ 3.1 g.cm -3 , with an average mass surface S mass =4000 cm².g -1 . This gives us a specific surface of S sp =12000 cm².cm -3 . Moreover, if we assume a disconnected morphology, we will set a positive Euler characteristic, and finally, we will set the size of the cube, T = 100µm (it is the usual size of a representative element of volume (REV) in such modelisations, as it can be found in [4] or [START_REF] Šmilauer | Microstructured-based micromechanical prediction of elastic properties in hydrating cement paste[END_REF]).

Setting one by one the previous parameters, we can deduce the ratio u/σ² of the random field needed and the correlation length L c . In an ideal way, we should also set the volume fraction, which would lead to setting the variance σ but the strong non-linearity of equation ( 3) does not allow us to do so. So far we are limited to a maximum volume fraction around 30%, which for a cement paste represent a w/c of 0.7. We have therefore to assume arbitrary a variance σ = 1.

With all these parameters, we are now able to generate a realisation of a random field (Figure 3, left), and by thresholding it, we will create an excursion set (Figure 3, right) that will represent anhydrous cement. The cube's complementary will be set as water.

3.3

Hydration process Once an initial morphology is created, we use equation ( 3) to make evolve the threshold according to equation ( 4). So for a given hydration degree, we will know the exact volume fraction of each phase, and thus the thresholds we have to apply to yield the different phases. As we can see on Figure 4, we can follow the evolution of the morphology according to the hydration degree.

APPLCATIONS

We present here some possible applications for this morphological framework, with first the transition from morphology to FE truss, then an application aimed to estimate the Young modulus evolution over the hydration degree increase, and finally the estimation of an effective diffusion coefficient.

Embedded discontinuities

Once the morphologies are created, we need to "project" them on a FE truss, in order to make any kind of calculation. It has been chosen to use not adapted meshes, not to constrain nodes spatial position by the evolving morphology. It means that both mechanical and geometrical properties have to be handled inside some interface elements. These elements that are between two phases are split in two parts, each having different material properties by enhancing them with strain (weak) discontinuities [START_REF] Ortiz | A Finite Element method or localized failure analysis[END_REF]. We choose to use the Embedded Finite Element Method (E-FEM), to achieve this enhancement.

For the initial morphologies, we have two different materials, which mean we have three types of elements: the ones only in anhydrous cement, the ones only in water, and the ones between the two (cut elements). Following, for a given hydration degree, we will have three different materials, and three different cut elements.

To set properly the materials properties, a projection of the different excursion sets are made onto the trusses. That allow us to have a multiphase material, and each phase has its own material properties.

4.2

Young modulus evolution One possible application of this framework is presented here, within a sequenced multiscale approach. For a given cement paste, we can follow a property over hydration. For example, we estimate here the Young modulus evolution. We enter properties at a microscopic scale (each material is then different), and we are computing for a macroscopic property. With the material properties set in Table 1, and by doing first a tri-tension test, and then a pure shear test, for each step of hydration, we can estimate the Young modulus and Poisson's ratio. Timoshenko's bi-phasic beam elements are used, in order to take into account tension, compression and shear.

As we can see on Figure 4, the macroscopic Young modulus is growing continuously over hydration degree. The overall shape of the curve is quite representative for such w/c ratio [121]. A slight raise of the slope can be noticed after α=0.4. The final value is not very representative. However, one have to remember that such a cement will never go further than α=0.8, and that the volume fraction of hydration products is overestimated by all the assumptions made earlier. Besides, the percolation threshold seems too high compared to available experimental results [START_REF] Stefan | Prediction of elastic properties of cement pastes at early ages[END_REF]. The Young modulus values for anhydrous cement and hydration products are weighted mean between all the chemical components (Portlandite, Afm, Aft and C-S-H for a CEM-I cement) obtained by nano-indentation.

Effective diffusion coefficient

Another type of application is to focus on hydrated cement paste, and estimate transport properties. An application is presented here for the diffusion coefficient.

Same way as previously, we can set material parameters for each phases, here effective diffusivity (D e in Table 1). This time we use bi-phasic bars with one degree of freedom. Then by applying unitary gradients in each direction of the REV, we can deduce the diffusivity matrix, which should be symmetric, and isotropic.

Using the Fick's law:

Q = D . C (5)

Figure 1 :Figure 2 :

 12 Figure 1 : Example of a one dimensional excursion set (A u )

Figure 3 :

 3 Figure 3 : Random field (left) and Excursion set (right, representing anhydrous cement)

Figure 4 :

 4 Figure 4: Morphologies for different hydration degrees (... water, ... hydration products, ... anhydrous cement)

Figure 4 :

 4 Figure 4 : Young modulus evolution for a cement paste (w/c = 0.7)

C 1 = 11 1

 111 [1 0 0] ; C 2 = [0 1 0] ; C 3 = [0 0 1] we can compute the resulting flux Q, for each C i , and thus reconstruct matrix D e : D e = o p.10 -There is a numerical limitation for ν (0.1<ν<0.275), due to the beam elements.

Table 1 :

 1 Material properties[START_REF] Stefan | Prediction of elastic properties of cement pastes at early ages[END_REF]

	Property	E (MPa) ν 1	D e ( m².s -1 )
	Anhydrous cement	135000 0.25	1,00.10 -14
	Hydration products	25000 0.25	1,00.10 -11
	water	1	0.27	1,00.10 -07

We have, according to theory, a matrix (and thus a REV) more or less isotropic, with an order of magnitude between diagonal and extra-diagonal terms.

CONCLUSIONS

This work presents a morphological framework for modeling cement paste like materials within a sequenced multi-scale approach. Morphologies are yielded by thresholding realization of random fields. Then, by their projections, we present two different applications relying on non adapted meshes. As we pointed it out, efforts have still to be made on order to have more realistic morphologies (volume fractions), and the analytical model for the hydration process needs to be more accurate, taking into account more phases.


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

Vincenzo Patera

Simon Rit

Angelo Schiavi

Ilaria Rinaldi

Regularised patient-specific

INTRODUCTION

Proton beam therapy offers high selectivity in dose deposition due to the protons' distinct depth dose profile, commonly referred to as Bragg curve. This feature potentially allows for high dose to the tumour while sparing healthy surrounding tissue [START_REF] Durante | Nuclear physics in particle therapy: a review[END_REF]. To fully exploit this advantage, precise knowledge of the relative stopping power (RSP) of the patient tissue is needed to correctly predict the proton range in the treatment planning system (TPS), in the best case on the day of treatment.

MATERIALS AND METHODS

Currently, proton treatment planning is performed based on X-ray computed tomography (CT) data typically acquired some days before the first delivered fraction. The X-ray CT is expressed in Hounsfield Unit (HU) and approximately converted to RSP [START_REF] Schneider | Patient specific optimization of the relation between CT-Hounsfield units and proton stopping power with proton radiography[END_REF]. The uncertainties associated with the conversion process require additional safety margins around the contoured treatment volume [START_REF] Jäkel | Relation between carbon ion ranges and x-ray CT numbers[END_REF][START_REF] Paganetti | Range uncertainties in proton therapy and the role of Monte Carlo simulations[END_REF]. Proton radiography and proton CT have already been investigated a long time ago [START_REF] West | Radiography with 160 MeV protons[END_REF][START_REF] Koehler | Proton Radiography[END_REF][START_REF] Steward | Proton Beam Radiography in Tumor Detection[END_REF][START_REF] Takada | Proton computed tomography with a 250 MeV pulsed beam[END_REF]) and more recently been brought up in the context of proton beam therapy [START_REF] Schneider | Proton radiography as a tool for quality control in proton therapy[END_REF][START_REF] Schulte | Conceptual design of a proton computed tomography system for applications in proton radiation therapy[END_REF]. In analogy to X-ray based transmission imaging, protons are shot through the patient and an image is formed based on the residual energy or residual range measured by a suitable detector positioned behind the patient [START_REF] Parodi | Heavy ion radiography and tomography[END_REF][START_REF] Poludniowski | Proton radiography and tomography with application to proton therapy[END_REF][START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF]. A full tomographic scan directly provides a volumetric RSP map of the patient and therefore bypasses any conversion.

It could thereby allow for reduced safety margins [START_REF] Arbor | Monte Carlo comparison of x-ray and proton CT for range calculations of proton therapy beams[END_REF]. [START_REF] Schneider | Patient specific optimization of the relation between CT-Hounsfield units and proton stopping power with proton radiography[END_REF] proposed an approach based on a single proton radiography expressing the water equivalent thickness (WET) of the patient which is arguably easier to achieve in practice and exposes the patient to a lower dose than a full proton CT scan. The idea is not to circumvent the conversion, but to perform a patient-specific calibration of the CT-RSP curve by comparing a proton radiographic image with a proton digitally reconstructed radiography (pDRR). This optimised curve would then be used for the treatment planning instead of the clinical one.

More recently, the method has been studied in combination with two types of proton imaging set-ups. [START_REF] Doolan | Patientspecific stopping power calibration for proton therapy planning based on single-detector proton radiography[END_REF] used a single plane detector in double scattering beam delivery.

In this case, the beam is gradually pulled back by a range modulator device while the pixel detector behind the patient records the signal created by many protons at once. On the other hand, in single proton tracking, the imaging set-up records information proton by proton and an estimate of individual proton trajectories through the patient can be derived, the so-called most likely path [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] Collins-Fekete | Pretreatment patient-specific stopping power by combining list-mode proton radiography and x-ray CT[END_REF] have investigated the patient-specific calibration in combination with this kind of system.

In this work, we consider a proton imaging set-up which combines pencil beam scanning with a range telescope [START_REF] Rinaldi | Experimental characterization of a prototype detector system for carbon ion radiography and tomography[END_REF][START_REF] Parodi | Heavy ion radiography and tomography[END_REF][START_REF] Farace | Pencil beam proton radiography using a multilayer ionization chamber[END_REF]. It can be employed in clinical practice potentially using available quality assurance (QA) equipment and without further hardware modifications in the treatment room [START_REF] Krah | A comprehensive theoretical comparison of proton imaging set-ups in terms of spatial resolution[END_REF]. It is therefore of interest to develop the mathematical tools necessary to perform a patient-specific CT-RSP calibration in combination with such a set-up. In particular, we develop a regularised optimisation procedure which, as it turns out, is necessary to obtain plausible calibration curves. We perform a Monte Carlo study to assess the accuracy and robustness of the proposed methods. Although we focused on a specific kind of imaging system, most of our methods can be adapted and applied to other proton imaging set-ups as well.

MATERIALS AND METHODS

Patient specific calibration -general principle

The CT-RSP conversion curve is typically modelled as a series of linear segments parametrised by a set of parameter tuples (H α , R α ), as illustrated in the centre of figure 1 [START_REF] Schneider | The calibration of CT Hounsfield units for radiotherapy treatment planning[END_REF][START_REF] Taasti | Inter-centre variability of CT-based stopping-power prediction in particle therapy: Survey-based evaluation[END_REF]. The patient-specific calibration optimises the values of R α by minimising the difference between a pDRR and a measured proton radiography. A brute force implementation would start from an initial guess for the CT-RSP curve, convert a planning CT (expressed in Hounsfield Units [HU]) to a three-dimensional (3D) RSP map, project it to produce the pDRR, compare with the proton radiography (e.g., by calculating the mean squared difference), update the set of R α values, and iterate again. We will present a method to separate the projection operation from the actual optimisation step.

HU RSP

Proton radiography (pR)

Optimize R 1 , R 2 , R 3 ...

Framework of this study

Naturally, an important question to be asked would be: "How accurate does the patient-specific calibration method reproduce the HU to RSP relation of the patient's tissue?" Answering this experimentally on the basis of phantom measurements would require precise knowledge of the phantom's spatial material distribution. For realistically complex anthropomorphic or animal tissue phantoms, this information is difficult to have, although proton radiographic measurements have been performed in phantoms of known material composition in some studies [START_REF] Farace | Pencil beam proton radiography using a multilayer ionization chamber[END_REF], Wohlfahrt et al. 2018). Phantoms which are geometrically less complex than a patient, on the other hand, might lead to unrealistic conclusions. For this reason, we remained entirely within a Monte Carlo framework for our study where we could use a digital patient model as ground truth. It was derived by converting a patient CT into RSP before importing it into Monte Carlo (see section 2.4). The CT-RSP curve used for this conversion was therefore the ground truth reference curve and we evaluated how well it was recovered by the optimisation procedure. At no point do we attempt to directly compare the calibration outcome with some clinical curve. This would require quantifying or modelling the typical discrepancies between the true tissue properties and those predicted by a clinical curve, and it was not the scope of this work.

Proton radiographies

We used Fred, a GPU-accelerated Monte Carlo code developed at University Rome La Sapienza, to simulate proton radiographies. The code has been validated against other MC toolkits, Geant4 (Agostinelli et al. 2003, Allison et al. 2016) and FLUKA [START_REF] Ferrari | FLUKA: A multi-particle transport code (program version[END_REF][START_REF] Böhlen | The FLUKA Code: Developments and Challenges for High Energy and Medical Applications[END_REF])), as well as experimental data from the CNAO proton therapy facility in Italy [START_REF] Schiavi | Fred: A GPU-accelerated fast-Monte Carlo code for rapid treatment plan recalculation in ion beam therapy[END_REF] and, more recently, against measurements from the Krakow proton therapy facility [START_REF] Garbacz | Proton Therapy Treatment Plan Verification in CCB Krakow Using Fred Monte Carlo TPS Tool[END_REF]. Its advantage for our study was the high computation speed.

A patient CT image was imported into Fred and to each voxel, an elemental composition, density, RSP, and radiation length was assigned as explained in section 2.4. No segmentation into tissue compartments was performed. We modelled the range telescope as a homogeneous box of water in which the integrated depth dose profiles were scored in depth increments of 2.3 mm. The geometry used in Fred is shown in figure 2 as well as a proton radiography of the patient used in this study.

The irradiation plan contained a two dimensional array of pencil beams with 10 3 protons per spot and a beam energy of 200 MeV, i.e., high enough so that the protons fully traversed the patient and partly penetrated the range telescope. For simplicity, we used a parallel beam geometry and kept the beam axis aligned with the voxel grid. Furthermore, we chose the beam spot positions to be centred on the voxel surfaces. The initial beam profile was assumed to be a symmetric Gaussian distribution characterised by the full width at half maximum (FWHM) parameter.

For each scanned beam spot, the range telescope recorded a (discretised) Bragg curve from which the water equivalent thickness (WET) was determined using a data processing procedure detailed in [START_REF] Krah | A comprehensive theoretical comparison of proton imaging set-ups in terms of spatial resolution[END_REF]. From the list of beam spot coordinates and WET values, the proton radiography was constructed as 2D image so that each spot corresponded to one pixel in the proton radiography.

Conversion of CT to RSP, density, and elemental composition

Fred interpolates all parameters (density, RSP, radiation length, elemental composition) linearly between the materials provided in a lookup We built an additional reference lookup table (table 2) using the material compositions and HU intervals reported in table 6 in [START_REF] Schneider | Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions[END_REF]. The density was determined using the fit shown in their figure 9. The stopping power and radiation length values were calculated 155 based on the elemental composition exploiting Geant4's functionality (Agostinelli et al. 2003).

156

We used this table for some of the presented results to determine the regulariser weights (see 157 section 2.10.1).

(R α , R α+1 ) if h ∈ (H α , H α+1 ) and between (R α-1 , R α ), if h ∈ (H α-1 , H α ).
175

This motivates the definition of a vector w whose components w α are calculated in the 176 following way:

w α (h) ≡ h -H α-1 H α -H α-1 Ω α (h) + H α+1 -h H α+1 -H α Ω α+1 (h) w 0 (h) = H 1 -h H 1 -H 0 Ω α=1 (h) w N (h) = h -H N -1 H N -H N -1 Ω α=N (h) (1)
with

Ω α (h) = 1 for H α-1 ≤ h < H α 0 otherwise. ( 2
)
The value of w α thus quantifies how much r depends on the corresponding parameter R α .

Note that most components of the vector w are zero except for two (or one if h lies exactly on an interval boundary) and that the sum over all w α is always one. The RSP value r is obtained by contracting w with the vector R containing the parameters R α of the conversion curve:

r = N α=0 w α R α = w • R. (3)
When applying the piecewise barycentric interpolation to all voxels in the CT image, the three dimensional HU matrix, h ijk , is mapped onto a four dimensional coefficient matrix w ijk,α .

In accordance with equation 3, the 3D RSP map is:

r ijk = N α=0 w ijk,α R α = w ijk • R (4)

Proton digitally reconstructed radiography

The pDRR is obtained by projecting the RSP map along lines corresponding to the pencil beams used to acquire the proton radiography. We do not consider projection directions oblique to the CT grid in this work (see section 2.3). Denoting with proj k the projection operation along axis k, the pDRR is calculated as

D ipjp = proj k,ipjp (r ijk ) = proj k,ipjp N α=0 w ijk,α R α = N α=0 proj k,ipjp (w ijk,α ) R α = N α=0 A ipjp,α R α (5) with A ipjp,α = proj k,ipjp (w ijk,α ). ( 6
)
The projection operation is wrapped into the matrix A because it only acts on the coefficients w ijk , but not on the parameters R α , and may therefore be interchanged with the contraction over the R vector. For an ideal pencil beam with FWHM=0 and neglecting scattering, the projection is simply the sum along k multiplied by the voxel size s, i.e., A ip,jp,α = s k w i=ipj=jpk,α . A more realistic projection operation must at least take into account the beam size and possibly the effect of scattering. In this work, we investigated the following three projection models:

(1) Gaussian beam without scattering, (2) Gaussian beam with material and depth dependent scattering, and (3) Monte Carlo based projection.

Gaussian beam without scattering

In this model, the beam is described by a Gaussian transverse profile characterised by the FWHM, which is 5-10 mm for clinical proton beams at the energies used for proton imaging (≈ 200 MeV). The effect of scattering within the patient is neglected and the beam profile is considered to be constant in depth and independent of the pencil beam position. Mathematically, this corresponds to a convolution with a shift invariant Gaussian kernel G ip-i,jp-j in the plane orthogonal to the proton beam. The matrix A ip,jp,α is then calculated as

A ip,jp,α = s k i,j G ip-i,jp-j w ijk,α = s i,j G ip-i,jp-j k w ijk,α . (7)
The convolution needs to be performed only once after summing along k, making the implementation computationally inexpensive.

Gaussian beam model with scattering

This more realistic projection model takes into account the effect of scattering which leads to a gradual widening of the proton beam with depth. The most important contribution is due to MCS, i.e., a succession of many small angle deflections as a proton propagates through the medium. We take as depth dependent beam size a combination of the initial beam size and the contribution due to MCS in the medium,

σ beam (z k ) = σ 2 init + σ 2 mcs (z k ), (8)
because the Gaussian pencil beam can approximately be thought of as a weighted superposition of many pin-like pencil beams, each of which undergoes MCS. We denote with z k the depth of the slices in the CT voxel grid along the beam axis.

We used a discretised version of the analytical formulae from [START_REF] Lynch | Approximations to multiple Coulomb scattering[END_REF] and [START_REF] Gottschalk | Multiple Coulomb scattering of 160 MeV protons[END_REF] to calculate σ 2 mcs as a function z k :

σ 2 mcs (z k ) = Ω 2 0 1 + 0.038 ln z k X 0 s k k =0 z 2 k β 2 (z k )p 2 (z k ) 1 X 0 (z k ) , (9)
with Ω 0 = 13.6 MeV and s the voxel size. The radiation length X 0 depends on the type of material and parametrises the amount of MCS the protons suffer while traversing it. In a heterogeneous medium like a patient, it depends on the three dimensional position within

MATERIALS AND METHODS

the patient volume. The terms β and p are the proton's velocity relative to the speed of light, β = v/c, and its momentum, respectively. They relate to the proton energy as βp = (E 2 +2E E p )/(E+E p )c, where E p is the proton's rest mass, and decrease as a proton continuously slows down while traversing the patient.

Evaluating equation 9 requires some prior knowledge about the material composition within the patient which we obtained by converting the patient CT to an RSP and X 0 map through the default tables implemented in Fred. Alternatively, one could perform a first optimisation of the CT-RSP curve using the Gaussian model without scattering, use the result to evaluate equation 9, and optimise the CT-RSP curve again, this time using the Gaussian beam model with scattering.

We make the following approximations to be able to sum along the beam axis z only in equation 9. For a given pencil beam, we determine the mean radiation length the protons encounter in each slice of the CT volume by calculating the weighted average value of X 0 around the beam center (i p , j p ) in the (i, j) plane at depth z k . We use a two dimensional Gaussian weighting kernel with a standard deviation reflecting the beam spot size. This yields the depth dependent term X 0 (z k ) in equation 9. Similarly, we calculate the water equivalent thickness z WET that the protons belonging to a given pencil beam have traversed on average when penetrating the patient up to some geometric depth z k . We then use the expression [START_REF] Bortfeld | An analytical approximation of depth -dose distributions for therapeutic proton beams[END_REF] to calculate the average proton energy in depth z k . In this formula, R refers to the proton range in water and χ and ξ are fit constants.

E(z WET ) = ((R -z WET )/χ) 1/ξ from
We determine the approximate beam size for each pencil beam (i p , j p ) and in all CT slices k as explained above to build a space variant Gaussian convolution kernel G ipjp;ij (z k ) of width σ beam,ipjp (z k ). The matrix A ipjp,α is calculated as

A ipjp,α = s k i,j G ipjp;ij (z k )w ijk,α . (10)
Fast Fourier transform based convolution algorithms cannot be used here because the kernel is shift variant. Instead, the multiplication and summation operations in equation 10 must be executed explicitly for all HU intervals α. We therefore implemented equation 10 on an Nvidia Titan X GPU card using the OpenCL library. This helped lowering the computation time from several hours per projection down to about 15 sec.

Monte Carlo based projection

We used Monte Carlo to perform the projection in equation 5 to avoid the approximations inherent in the analytical beam model and to include the effect of nuclear interactions on the proton trajectories. To this end, we implemented a dedicated plugin to Fred which iterates through the following scheme: Initially, the patient CT is loaded and the components of the matrix ω ijk,α are calculated. The pencil beams are defined and arranged in the same way as in the proton radiography simulation (see Sec. 2.3). The protons are then transported across the patient CT volume using the physics models of the MC code. After each transport step (before traversing into a new voxel the latest), the coefficient vector ω is extracted from the current voxel and multiplied by the step length projected onto the beam axis z. This information is cumulated along the entire trajectory and for all protons within one pencil beam. The result is the integrated coefficient vector ω associated with one pencil beam and repeating the procedure over the entire array of pencil beams yields the matrix A ipjp,α .

If the number of protons per pencil beam remained constant across the patient volume, the integrated coefficient vector would simply be normalised by this number. In reality, the number of protons fluctuates as a function of depth because (1) secondary protons are created through inelastic scattering or knock-on elastic collisions, (2) protons stop inside the patient because they lose a large fraction of their energy in a nuclear interaction, and (3) some protons are scattered at large Rutherford-like angles and exit the patient laterally without fully traversing it. Our plugin accounts for these effects by applying a depth dependent normalisation and only includes protons which would have reached the detector. In the current implementation of Fred, the plugin can only be executed in a single thread on a CPU which makes the computation time much longer (order of hours) than for the analytical projection model.

Linear minimisation problem

Henceforth, we use a single index l for the pixels in the proton radiography (see section 2.5).

We define the data attachment term of the cost function as

E data ( R) = l P l - N α=0 A lα R α 2 = P -A R 2 with A lα = proj k (w lk,α ), (11)
where 2 refers to the norm. We have introduced compact matrix/vector notation on the right hand side.

The optimal parameters of the CT-RSP conversion (denoted with the ˆsymbol) curve are the

Parameter correlation coefficients

The minimisation problem in equation 12 is mathematically well posed if matrix A has full rank, or, in practical terms, if there are at least as many pixels in the proton radiography as there are R α parameters to be optimised. Even if this criterion is met, the optimised CT-RSP curve would not be plausible from a physics point of view in many cases, as we will show in section 3.1.

The reason is the following: The WET value encoded in a single pixel of the pDRR results from a summation of RSP values in many voxels along and around the proton beam (see section 2.7).

Increasing the value of a parameter R α increases the RSP in some of the contributing voxels and consequently the average WET. At the same time, lowering the value of another parameter R α tends to lower the WET value. Therefore, the R parameters are expected to be correlated. We also note that neighboring R α parameters are in any case anti-correlated because the CT-RSP curve interpolates linearly between them and is continuous by construction. The correlation coefficient of the parameters R α can determined from the covariance matrix of the cost function as follows:

ρ αα = Cov αα √ Cov αα Cov α α with Cov αα = (A T A) -1 αα . (13)
The coefficient values always lie within the [-1, 1] interval and ρ αα = 1 by construction.

MATERIALS AND METHODS

Regularised optimisation

As we will show in the results, correlation among the parameters R α would lead to unrealistic optimised CT-RSP curves. Therefore, we introduce a regularisation term in the cost function to control the variation of the CT-RSP curve. It penalises the intrinsic variation of the patientspecific CT-RSP curve. We define

E reg = α M α=α 0 (d α+1 -d α ) 2 , ( 14
)
with the discrete derivative d α of the curve calculated as

d α = R α -R α-1 H α -H α-1 . ( 15
)
The regularisation may be applied to a subset of HU intervals only (see section 2.10.1) and α 0 and α M refer to the first and last of them, respectively.

Equation 14 can be written as matrix equation,

E reg = α M α=α 0 A reg α α R α 2 = A reg R 2 , ( 16
)
where the matrix components are given as

A reg α α = 1 H α+1 -H α-1 δ α ,α-1 -2δ α ,α + δ α ,α+1 , with δ α ,α = 1 for α = α 0 otherwise. (17)

Regularised cost function

The conversion curve is expected to be rather smooth and monotonic in the regions HU< -100

and HU> 100, but should be allowed to fluctuate slightly in the region of soft tissue, i.e., -100 < HU < 100 [START_REF] Schneider | Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions[END_REF]. For this reason, we applied the regularisation separately to the three intervals.

The regularised cost function is the sum of the data term (equation 11) and three regulariser terms (equation 16). It takes the form

E( R) = E data + γ air E reg,air + γ soft E reg,soft + γ bone E reg,bone . (18)
We have used the extra subscripts "air", "soft", and "bone" as short hand notation to refer to the three distinct HU regions of the conversion curve subject to the regularisation, with the indices α 0 and α M in equation 16 chosen appropriately to confine the selected HU intervals.

We used the python package scipy.linalg.lstsq to minimise the cost function, which is based on singular value decomposition.

The γ factors set the relative weight of the regularisers compared to the data term and their value should be chosen such that the optimised curve contains a plausible level of smoothness.

We adopted the following method to automatically select the weights γ: Starting with an initial guess, we ran the optimisation and then calculated the magnitude of each regulariser term once for the optimised CT-RSP curve and once for a reference curve. We increased the respective weight γ, if the regulariser term of the optimised curve was greater than that of the reference curve, and lowered it in the opposite case. This procedure was repeated iteratively until all regulariser terms of the optimised and reference curve corresponded to within a few percent. In this way, the optimised CT-RSP curve had a comparable amount of variation as the reference.

The reference curve should be chosen (or constructed) such that it includes a level of RSP variation which is thought to be necessary to reflect the tissue variability within one patient.

The relative weights need to be set for each optimisation because the data term to which they refer is linked to a specific proton radiography (and the choice of a region-of-interest -see section 2.12). Note also that the regulariser does not force the absolute RSP values of the optimised and reference curves to resemble each other, but only their degree of smoothness.

Therefore, optimised RSP may still vary strongly among patients even in the case of a very smooth reference curve which suppresses RSP variation within one patient. In this work, we used the ground truth CT-RSP curve as reference to fix the regulariser weights, except for figure 4.

Assessment of range accuracy

The patient-specific CT-RSP calibration optimises the parameters R α by minimising the difference between predicted and measured average residual range for a field of pencil beams which fully penetrate the patient. From the treatment point of view, it is of interest to estimate how accurately the optimised curve predicts the range of lower energy protons stopping inside the patient.

To quantify the range accuracy which would be obtained with an optimised CT-RSP curve, we performed the following kind of Monte Carlo experiment. We defined a volume inside the patient's head and built an irradiation plan containing 10,000 pencil beams with randomised positions and beam energies. The positions were constrained by the lateral extension of the volume and the energy was chosen so that the pencil beams would penetrate the patient at least 20% of its total water equivalent thickness at that point and at most 80%. We scored the integrated depth dose profiles individually for all pencil beams. The same plan was simulated twice: once using the ground truth CT-RSP curve (see sections 2.2 and 2.3) to convert the patient CT to RSP and once using the optimised one. In this way, we obtained 10,000 pairs of integrated depth dose profiles. From each of them, we estimated the proton range by determining the water equivalent depth in which the distal edge decayed below 80 % of the peak value (Paganetti 2012a).

The profiles were linearly interpolated to reach subvoxel resolution. We built histograms out of the 10,000 range differences and determined the mean value as estimate for the range accuracy.

We note that for certain pencil beams directed towards complex structures (e.g. nose and air cavities), the 80% distal fall-off technique may suffer from range mixing and therefore lead to additional variation of the range differences. Manual inspection of a random subsample of profiles led us to conclude that such cases are rare enough not to compromise the results in terms of average range accuracy in the Monte Carlo experiment.

Region of interest in the proton radiography

To perform the CT-RSP curve optimisation, it suffices to acquire a proton radiography of a region of interest of the patient instead of, e.g., the entire head to perform the patient-specific calibration as long as the image contains a sufficient amount of pixels. The region could for example be chosen such to reduce the dose given to healthy tissue. To investigate how much the optimised CT-RSP curve depends on such a choice, we simulated nine proton radiographies capturing differently positioned (but equally sized) rectangular areas (shown in figure 9) and calculated an optimised CT-RSP curve for each of them.

In practice, the region imaged by the proton radiography could be selected in accordance with the treatment volume or, instead, be placed so that the imaged volume does not correspond to the treated one, e.g. to avoid giving dose to an organ at risk distal to the treatment volume during image acquisition. We performed the Monte Carlo experiment (see section 2.11) to investigate whether this would have an impact on the range accuracy. Specifically, we selected a region of interest, ROI exp , and defined the randomised irradiation plans so that the lateral coordinates of all pencil beams lay within it. We used the CT-RSP curve previously optimised based on another (or the same) region of interest proton radiography, ROI opti , as look up table

RESULTS

Figure 3: Left: Optimised CT-RSP curves obtained without (blue) and with (red) regulariser applied. The black curve is the ground truth reference used to simulate the proton radiography.

The two dashed vertical lines indicate the three separately regularised HU regions "air", "soft tissue", and "bone" (see section 2.10.1). Right: Relative difference between the DRRs calculated using the regularised and non-regularised CT-RSP curves, respectively.

in the Monte Carlo. We then determined the range accuracy as explained in section 2.11 for each pair (ROI exp , ROI opti ).

RESULTS

Regularisation

In figure 3 left, we show a patient-specific CT-RSP with and without regularisation used in the optimisation process. The proton radiography was simulated with an 8 mm FWHM beam and we used a rectangular region of interest indicated as black box in the right panel containing 5500 pixels. The projection was performed using the Gaussian beam model with scattering (see section 2.7.2). The patient-specific calibration curve contained 33 (H α , R α ) pairs to be optimised.

It is evident that the variation in the non-regularised curve is unacceptable, while it is kept at a plausible level when employing regularisation. The right panel shows the relative difference between two pDRR images: one calculated using the regularised CT-RSP curve and the other one using the non-regularised version. The two DRRs agree to better than 3% although the two conversion curves themselves differ by much more. This underlines that the non-regularised solution, despite its unrealistic variation, does indeed optimise the residual range similarly well as the regularised one.

The left panel of figure 4 shows optimised calibration curves for which the regulariser weights γ in equation 18 were determined based on two different reference CT-RSP relationships (see section 2.4). The regulariser terms have consequently enforced a degree of smoothness in the optimised curves similar to each reference curve. At the same time, both optimised curves approximately recovered the RSP values of the ground truth curve. This demonstrates that the regulariser does not affect absolute RSP, but only the RSP variability within the curve (see section 2.10.1).

The right panel shows optimised CT-RSP curves for which different HU intervals were used.

The five HU thresholds in the yellow curve were [-1000, -100, 0, 100, 1376]. The γ weights were determined based on the dashed reference curve and consequently the optimised curves are smoother than the ground truth curve. It is interesting to note that the yellow curve does not exhibit any unrealistic variation despite the fact that no regularisation was applied.

The patient-specific calibration curves shown in figure 5 are based on proton radiographies simulated with different beam sizes, where "pin like" means FWHM=0 mm. We used the Gaussian beam model with scattering to perform the projection of the pDRR. The optimisation was regularised as explained in section 2.10. No regularisation was used for the green curve. The lower panel shows the relative difference in percent between the ground truth CT-RSP and the optimised ones. We have omitted the data point for HU = -1000 because dividing by a very small ground truth RSP value (0.012) makes the relative difference large and meaningless.

Figure 6 illustrates the parameter correlation matrix ρ αα (equation 13) for three different beam sizes. We have used the Gaussian beam model with scattering (see section 2.7.2) to calculate the pDRR. The central diagonal contains values of 1 by construction, while the two diagonals off-set to the left and right are close to -1 because R α parameters of neighbouring HU intervals of the conversions curve tend to be anti-correlated. Non-zero correlation coefficients further away from the central diagonal occur because RSP information from many voxels is mixed during the projection operation (see section 2.9). Indeed, they are much more prominent in the left panel corresponding to an 8 mm FWHM beam than in the right panel, which refers to a pin-like beam.

In figure 7, we present optimised calibration curves obtained from proton radiographies simulated with different scattering processes activated. The lower panel shows the relative difference in percent between optimised and ground truth CT-RSP curves. "Full physics" includes all processes currently available in Fred, i.e., MCS and nuclear elastic and inelastic scattering; for "No nuclear", only MCS was included; and "no nuclear, no MCS" corresponds to a pure continuous slowing down simulation in which protons travel on straight lines through the voxel geometry. Energy loss and straggling was left at default settings in all cases. The Gaussian beam model with scattering was used to calculate the pDRR except for "no nuclear, no MCS"

where the Gaussian beam model without scattering was employed. The FWHM=8 mm curves were regularised (see section 2.10), while the pin-like curves were not.

Figure 8 left shows optimised CT-RSP curves obtained with the three projection models (see section 2.7). In all cases, the proton radiography was simulated with 8 mm FWHM pencil beams and the optimisation was regularised. The Monte Carlo projection (see section 2.7.3) was performed with 1000 protons per pencil beam. The right panel shows the relative difference between the pDRR obtained with the Gaussian projection model with scattering and the Monte Carlo projection.

RESULTS

RESULTS

Figure 9: Left: regions of interest in the proton radiography used for the optimisation. Right: Optimised CT-RSP obtained using these regions of interest.

Range accuracy and region of interest

In the figure 9, we show optimised CT-RSP curves (right panel) using the different regions of interest in the proton radiography illustrated in the left panel. Colours in all panels match.

The optimisation was performed using the Gaussian projection model with scatter (see section 2.7.2) and employing the regularisation separately in three HU regions (see section 2.10.1).

The curves vary by a few percent among each other, but deviate no more than 4% from the ground truth curve. There appears to be some systematic under estimation in the range above 800 HU associated with bony tissue.

Figure 10 presents the results of the Monte Carlo range experiment explained in section 2.11.

The upper left panel shows an example distribution of range differences obtained from 10, 000 pencil beams. The blue distribution refers to the difference between predicted and true range of the pencil beams when using an optimised CT-RSP. The mean relative range difference quantifies the overall average range accuracy while the RMSE is indicative of the range precision. The yellow distribution was generated by repeating the Monte Carlo experiment twice, both times with the same ground truth curve, and it reflects the intrinsic range variations due to statistics in the depth dose profiles.

The upper right panel of figure 10 shows the mean relative range difference obtained when performing the Monte Carlo experiment in a region of interest, ROI exp , while optimising the CT-RSP curve based on another, ROI opti (see section 2.12). The diagonal elements in the graphic correspond to the case when the imaged area corresponds to the treated volume. On the other hand, for combinations such as ROI exp = 1 and ROI opti = 4 (see figure 9, left), there is no overlap between the imaged and the treated volume. The lower panels show the RMSE for the control experiment (left) and the optimised vs. ground truth runs (right). The fact that the RMSE tends to be slightly higher when the Monte Carlo experiment was performed in a relatively heterogeneous region (e.g. ROI 9) is in line with the remark in section 2.11 concerning the robustness of the 80% fall-off test. Overall, the RMSE associated with the optimised CT-RSP curves is less than 0.3% higher than the RMSE in the control run.

Figure 11 shows the mean relative range accuracy and the RMSE (as error bars) resulting from the Monte Carlo experiment using the optimised CT-RSP curves based on DRRs calculated with the three beam models (see section 2.7). We performed both, the optimisation and the Monte Carlo experiment, in the same nine regions of interest shown in the left panel of figure 9.

The RMSE is only slightly larger than the RMSE of the control run (ground truth vs. ground truth; not shown), i.e. it is mainly due to Monte Carlo statistics and variation of the range retrieved by means of the 80% fall-off criterion (see section 2.11). Overall, the Gaussian beam

DISCUSSION

model with fixed beam size seems to introduce a slight overestimation of the proton range. The other two beam models, on average, reproduce the ground truth range faithfully to within 0.2-0.5%. Whether or not the minimal positive bias observed with the variable FHWM model is statistically significant appears difficult to establish based on the nine data points.

DISCUSSION

In the presented work, we introduced methods for the patient-specific calibration of the CT-RSP conversion curve based on proton radiography, with a special regard to proton imaging set-ups which combine pencil beam scanning with a position-insensitive range telescope. Our results demonstrate that the optimisation of a piecewise linear conversion curve needs to be regularised because it otherwise leads to unrealistic variations in the optimised curve. The advantage of a quadratic norm is that the regulariser terms can be easily combined with the data attachment term keeping the optimisation computationally very feasible. With the proposed procedure to select the relative weights of the regulariser terms, the variation in the optimised curve corresponds to that observed in the stoichiometric curve used as reference.

Our results show that the main reason for the need to regularise the cost function is linked to the correlation among the parameters R α of the CT-RSP conversion curve, which in turn is strongly influenced by the beam size (figure 6). In a wider phantom than the one considered in this work, scattering becomes more important leading to increased correlation and potentially poorer optimisation results unless more strongly regularised. In case of a pin-like pencil beam, for example, the non-regularised optimisation yielded almost the same result as the regularised one (figure 5). This explains why Collins-Fekete et al. ( 2017), who investigated the patient-specific CT-RSP calibration based on proton radiographies acquired with single tracking set-ups, did not encounter similar problems: the integration in their case is performed along an estimate of individual proton trajectories with very thin uncertainty envelopes around them and the degree of parameter correlation is expected to be even smaller than with a pin-like pencil beam in our case. On the other hand, the set-up used by [START_REF] Doolan | Patientspecific stopping power calibration for proton therapy planning based on single-detector proton radiography[END_REF] detects many protons at the same time, as in our case, yet delivered in a passive field rather than by pencil beam scanning.

Parameter correlation is therefore expected to be even stronger because the effective beam size is larger. The authors did not report any variation in the optimised CT-RSP curves, such as seen in figure 3, although they did not regularise the optimisation. The reason is probably that they overrode HU values to a single value in each organ making the phantom geometry much more homogeneous. They furthermore excluded heterogeneous regions prone to range mixing.

We remark that such a region-wise homogenous RSP map might potentially over-simplify the patient geometry.

We note that the way of parametrising the CT-RSP curve has an impact on the correlation between parameters and thus on the need to regularise the optimisation. Figure 4 shows that no regularisation might be required when using relatively few HU intervals. On the other hand, the CT-RSP curve has then less granularity to resolve intra-patient tissue variability. In practice, a reasonable option would be to choose the HU intervals in accordance with the parametrisation of the clinical lookup table, possibly with some additional intervals for more granularity. The way of interpolating between HU intervals is also expected to affect the optimisation outcome. In our case, neighbouring R α values are intrinsically anti-correlated because the linear segments are assumed to form a continuous conversion curve. This is seen from the negative entries below and above the diagonal of the correlation matrix in figure 6. Alternative parametrisations include a piecewise linear interpolation without the constraint of continuity at the intervals boundaries or a piecewise constant conversion curve which assigns the same RSP value to all voxels whose HU values fall in the same interval. In the latter case, the R α parameters would be correlated only through the projection operation and consequently less variation would be expected in the optimised conversion curve even without regularisation. We remark that such a piecewise 5. CONCLUSION constant conversion curve assigns the same RSP value to all voxels whose HU values fall in the same interval and is therefore similar to segmenting the CT image into homogeneous material regions as in [START_REF] Doolan | Patientspecific stopping power calibration for proton therapy planning based on single-detector proton radiography[END_REF]. Investigating the impact of different kinds of interpolation was beyond the scope of this work.

Among the three projection models used in our work, the Monte Carlo based one (section 2.7.3) takes the scattering processes into account most accurately through the physics models implemented in the code. In line with this, differences between the pDRR obtained with the Gaussian projection model with scattering and the Monte Carlo projection are mainly visible in heterogeneous regions and along the outer perimeter of the patient's head (figure 8, right) where proton scattering has the greatest impact on the projected image. In the left panel of figure 8, the resulting optimised CT-RSP curve differs from the one obtained with the analytical beam model in its shape, but both lie within 3-5% of the ground truth curve. In terms of range accuracy (figure 11), our results suggest that the simple Gaussian beam model with a fixed beam size might lead to an overestimation of up to about 1%. The other two models reproduce the ground truth range expectation to 0.2-0.5% accuracy on average. Among the two, the Monte Carlo beam model seems to be slightly preferable over the variable FWHM Gaussian model, although we cannot deduce whether this observation is statistically significant. At the moment, the latter is much faster (< 1 min) to calculate thanks to our GPU implementation compared to the former (order of hours) which currently runs on single CPU only, although these technical issues could be improved in the future.

Remapping the HU values h in the CT voxels to coefficient vectors w (section 2.6) allowed us to separate the projection operation from the actual optimisation. This is particularly useful when using a cost function that requires iterative minimisation. We emphasise that the remapping can be used with proton radiographies acquired by any kind of set-up. For example, when using a single tracking proton scanner [START_REF] Collins-Fekete | Pretreatment patient-specific stopping power by combining list-mode proton radiography and x-ray CT[END_REF], the projection operation to construct the pDRR (section 2.7) would be performed by integrating the w ijk,α matrices along the estimated most likely path of every detected proton.

In our Monte Carlo experiment, the range accuracy achievable with patient-specific conversion curves optimised using the methods presented here was well below 1%. Whether or not the imaged region-of-interest is chosen in accordance with the treated volume did not systematically impact the range accuracy. This statement only refers to the optimisation process itself and does not take into account possible local variability in the tissue composition.

We chose to perform our study entirely in Monte Carlo because it allowed us to define a ground truth CT-RSP conversion curve. This approach is idealised in one regard: In reality, no one-to-one correspondence exists between the HU and RSP of the patient tissue. Our ground truth scenario is therefore an approximation because we generate the RSP map through a one-toone conversion process. Studying the impact of the tissue composition on range accuracy could be achieved by modelling the expected tissue variation [START_REF] Yang | Theoretical variance analysis of single-and dual-energy computed tomography methods for calculating proton stopping power ratios of biological tissues[END_REF][START_REF] Möhler | Range prediction for tissue mixtures based on dual-energy CT[END_REF][START_REF] Wohlfahrt | Dual-energy CT based proton range prediction in head and pelvic tumor patients[END_REF]) and perturbing the converted RSP value in each voxel. We emphasise that this was not the scope of this work. It should also be noted that the patient-specific CT-RSP calibration based on proton radiography cannot solve the issue of tissue variability within one patient because it intrinsically relies on a conversion curve. Only a full volumetric proton CT image would provide a voxel-wise RSP value and thereby circumvent the conversion process.

CONCLUSION

We developed a series of methods to perform a patient-specific calibration of the HU to RSP conversion curve based on proton radiography. We introduced a regulariser term in the cost function and showed that this is necessary to avoid unrealistic variation in the optimised curve.

The main reason is the correlation between the parameters of the piecewise linearly interpolated CT-RSP curve and the effect is more pronounced with larger beam size. We performed Monte

Carlo experiments and obtained a range accuracy of better than 0.5%. Our results further indicate that from the point of view of the optimisation procedure, no systematic loss of range accuracy is probably to be expected if the imaged region of the patient does not correspond to the treatment volume. This work underlines that it is possible to perform a patient-specific calibration of the CT-RSP conversion curve on the basis of a single proton radiographic projection acquired with a detector typically available in a treatment facility.

Figure 1 :

 1 Figure 1: Optimisation scheme for the patient-specific CT-RSP calibration. The plot in the center is a schematic illustration of the CT-RSP curve defined as set of linear interpolations.

Figure 2 :

 2 Figure 2: Left: Geometry used in Monte Carlo to simulate proton radiographies. The detector is idealised as solid block of water and the dose scored in 2.3 mm depth intervals. Right: Example of a simulated proton radiography.

  figure1). We use the symbol to indicate vector notation, e.g. R = (R 0 , R 1 . . . , R N ).

  piecewise linear character of the CT-RSP conversion curve to remap the CT 169 numbers [in HU] in the CT image onto a set of coefficients. This allows us to state the optimi-170 sation as a linear minimisation problem. Let h be the HU value in one voxel of the patient CT image. Its conversion to r depends on a parameter tuple (H α , R α ) if and only if h lies in the HU intervals left or right of H α , i.e., h ∈ (H α-1 , H α+1 ). Specifically, r is the linear barycentric interpolation between the parameter pair

  components of the vector ˆ R which minimises the objective function, i.e., ˆ R = argmin E( R) . (12) Note that the HU intervals H α are implicitly encoded in the matrix A and not subject to optimisation, i.e. the parameters ˆ R are optimal for a given discretisation of the CT-RSP curve. Optimising H α and R α alongside each other would mean to optimise both A and R in equation 11 -a much harder problem to solve, if solvable at all. In practice, the intervals H α can be chosen to yield a desired degree of granularity (e.g. smaller intervals in soft tissue region, fewer in bone region), or in accordance with a clinical calibration curve already in use in the facility. In this work, we used the same intervals as in the ground truth curve except for one dedicated test where different intervals where used in the optimization than in the proton radiography simulation (figure 4 in section 3.1).

Figure 4 :

 4 Figure4: Left: Optimised CT-RSP curves using two different reference curves to determine the regulariser weights γ. Right: Optimised CT-RSP curves using different HU-intervals for the parameterisation. The regulariser weights γ were determined in accordance with the smoothness of the reference curve based on[START_REF] Schneider | Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions[END_REF] (see section 2.4). For the curve with five HU thresholds only, these were placed at HU=-1000, -100, 0, 100, 1376. No regularisation was applied to this curve.

Figure 5 :

 5 Figure 5: CT-RSP curves obtained from proton radiographies simulated with different beam sizes. The lower panel shows the relative difference in percent between the optimised and the ground truth CT-RSP curve.

Figure 6 :

 6 Figure 6: Visualisation of the parameter correlation matrix ρ αα (equation 13) for three different beam sizes.

Figure 7 :

 7 Figure 7: CT-RSP curves obtained from proton radiographies simulated with different physics processes activated. The lower panel shows the relative difference in percent between the optimised and the ground truth CT-RSP curve.

Figure 8 :

 8 Figure 8: Left: CT-RSP curves obtained using three different projection models (see section 2.7). The lower panel shows the relative difference in percent between the optimised and the ground truth CT-RSP curve. Right: Relative difference between pDRR images obtained with the Gaussian projection model with scattering and the Monte Carlo based projection.

Figure 10 :

 10 Figure 10: Upper left: example distribution of relative range differences obtained from the Monte Carlo experiment (see section 2.11). Upper right: Mean relative range difference between true and predicted range using different regions of interest for the optimisation and the Monte Carlo experiment (see section 2.11). Lower panels: RMSE in control experiment (left) and between true and predicted range (right).

Figure 11 :

 11 Figure 11: Mean relative range accuracy obtained in the Monte Carlo experiment (see section 2.11). The CT-RSP curves were optimised for the nine regions of interest shown in figure 9 left and using DRRs calculated with three different beam models (see section 2.7).

  table according to the CT number in each voxel, i.e., it treats each tissue as a composition of two out of the given set of reference materials.

	Fred uses
	RSP, density, and radiation length to calculate electromagnetic energy loss and multiple Coulomb
	scattering (MCS), and the elemental composition to obtain cross sections of nuclear interactions
	(Schiavi et al. 2017). When simulating the proton radiographies, we used Fred's default set
	of materials reported in table 1 which is essentially table 4 in Schneider et al. (1996) with an
	additional column containing the calculated radiation length X 0 . Following the arguments in
	Sec. 2.2, we defined the relationship between CT numbers and RSP in this table as ground
	truth conversion curve. Some materials in the table refer to tissues which do not appear in a
	patient head and which would be disregarded when using the table in a TPS. For the purpose of
	providing a ground truth information, we kept them to have a relatively granular lookup table
	in the simulation.

Table 1 :

 1 Material table from[START_REF] Schneider | The calibration of CT Hounsfield units for radiotherapy treatment planning[END_REF] used to simulate the proton radiographies. It contains the CT-RSP relationship which we used as ground truth knowledge.

		O	C	H	N	Elemental fraction in % Ca P Na Mg S	Cl	K	Fe	I	ρ [g/cm 3 ] [HU] CT No. RSP	X 0 [cm]
	Air	23.5	0	0	76.5 0	0	0	0	0	0	0	0	0	0.0012	-1000	0.001062 30862.6
	Lung infl.	74.9	10.5 10.3	3.1	0	0.2	0.2 0	0.3 0.3 0.2 0	0	0.26	-741	0.257998 140.478
	Adip. tissue	27.8	59.8 11.4	0.7	0	0	0.1 0	0.1 0.1 0	0	0	0.95	-70	0.978975 43.3744
	Yellow marrow	23.1	64.4 11.5	0.7	0	0	0.1 0	0.1 0.1 0	0	0	0.98	-42	1.01303	42.5351
	Water	88.81 0	11.19 0		0	0	0	0	0	0	0	0	0	1	0	1	36.0823
	Breast	52.7	33.2 10.6	3		0	0.1	0.1 0	0.2 0.1 0	0	0	1.02	3	1.02898	37.8775
	Red marrow	43.9	41.4 10.5	3.4	0	0.1	0	0	0.2 0.2 0.2 0.1 0	1.03	14	1.04102	38.0317
	GI tract	75.1	11.5 10.6	2.2	0	0.1	0.1 0	0.1 0.2 0.1 0	0	1.03	23	1.02804	35.7108
	Lymph	83.2	4.1	10.8	1.1	0	0	0.3 0	0.1 0.4 0	0	0	1.03	28	1.02701	35.1365
	Testis	69.4	16.9 10.6	2.2	0	0.2	0.2 0	0.1 0.2 0.2 0	0	1.04	32	1.04104	35.7064
	Brain	71.2	14.5 10.7	2.2	0	0.4	0.2 0	0.2 0.3 0.3 0	0	1.04	37	1.04	35.3985
	Thyroid	74.5	11.9 10.4	2.4	0	0.1	0.2 0	0.1 0.2 0.1 0	0.1 1.05	40	1.04496	34.9107
	Muscle	71	14.3 10.2	3.4	0	0.2	0.1 0	0.3 0.1 0.4 0	0	1.05	42	1.04401	35.0533
	Kidney	72.4	13.2 10.3	3		0.1	0.2	0.2 0	0.2 0.2 0.2 0	0	1.05	43	1.04496	34.9957
	Lung defl.	74.9	10.5 10.3	3.1	0	0.2	0.2 0	0.3 0.3 0.2 0	0	1.05	44	1.04401	34.785
	Ovary	76.8	9.3	10.5	2.4	0	0.2	0.2 0	0.2 0.2 0.2 0	0	1.05	45	1.04601	34.7622
	Eye lens	64.6	19.5 9.6	5.7	0	0.1	0.1 0	0.3 0.1 0	0	0	1.07	50	1.06005	34.9361
	Liver	71.6	13.9 10.2	3		0	0.3	0.2 0	0.3 0.2 0.3 0	0	1.06	53	1.05396	34.6543
	Spleen	74.1	11.3 10.3	3.2	0	0.3	0.1 0	0.2 0.2 0.3 0	0	1.06	54	1.05396	34.5251
	Heart	74.5	11	10.2	3.3	0	0.1	0.1 0	0.2 0.3 0.2 0.1 0	1.06	55	1.053	34.5819
	Skin	64.5	20.4 10	4.2	0	0.1	0.2 0	0.2 0.3 0.1 0	0	1.09	75	1.08401	34.2994
	Cartilage	74.4	9.9	9.6	2.2	0	2.2	0.5 0	0.9 0.3 0	0	0	1.1	98	1.08097	32.4815
	Skel. spong.	36.7	40.4 8.5	2.8	7.4	3.4	0.1 0.1 0.2 0.2 0.1 0.1 0	1.18	260	1.15605	29.3503
	Skel. sacrum	43.8	30.2 7.4	3.7	9.8	4.5	0	0.1 0.2 0.1 0.1 0.1 0	1.29	413	1.23801	25.3857
	Vert. col. (D6,L3) 43.7	28.7 7	3.8	11.1 5.1	0	0.1 0.2 0.1 0.1 0.1 0	1.33	477	1.26696	24.0933
	Femur	36.8	34.5 7	2.8	12.9 5.5	0.1 0.1 0.2 0.1 0	0	0	1.33	499	1.26895	23.8711
	Ribs (2nd, 6th)	43.6	26.3 6.4	3.9	13.1 6	0.1 0.1 0.3 0.1 0.1 0	0	1.41	595	1.32907	22.0085
	Vert. col. (C4)	43.6	26.1 6.3	3.9	13.3 6.1	0.1 0.1 0.3 0.1 0.1 0.1 0	1.42	609	1.33693	21.7318
	Humerus	36.9	31.4 6	3.1	15.2 7	0.1 0.1 0.2 0	0	0	0	1.46	683	1.37006	20.8816
	Ribs (10th)	43.4	23.5 5.6	4		15.6 7.2	0.1 0.1 0.3 0.1 0.1 0	0	1.52	763	1.41299	19.6246
	Cranium	43.5	21.2 5	4		17.6 8.1	0.1 0.2 0.3 0	0	0	0	1.61	903	1.48007	17.9969
	Mandible	43.5	19.9 4.6	4.1	18.7 8.6	0.1 0.2 0.3 0	0	0	0	1.68	1006	1.53401	16.9642
	Cortical bone	43.5	15.5 3.4	4.2	22.5 10.3 0.1 0.2 0.3 0	0	0	0	1.92	1376	1.71398	14.0583

Table 2 :

 2 Lookup table according to[START_REF] Schneider | Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions[END_REF] with additional columns containing RSP (at 200 MeV) and radiation length X 0 .

	O	C	H	N	Elemental fraction in % Ca P Na Mg	S	Cl	K	Ar	ρ [g/cm 3 ]	CT No. [HU]	RSP	X 0 [cm]
	23.2	0	0	75.5	0	0	0	0	0	0	0	1.3	0.000605	-1050	0.000537909	60521.5
	74.9	10.5	10.3	3.1	0	0.2	0.2	0	0.3	0.3	0.2	0	0.052695	-950	0.0529346	693.293
	19.8	68.1	11.6	0.2	0	0	0.1	0	0.1	0.1	0	0	0.907347	-120	0.934753	46.3824
	30.8	56.7	11.3	0.9	0	0	0.1	0	0.1	0.1	0	0	0.943934	-82	0.967335	43.3455
	41.1	45.8	11	1.5	0	0.1	0.1	0	0.2	0.2	0	0	0.97083	-52	0.989457	40.9452
	50.9	35.6	10.8	2.2	0	0.1	0	0	0.2	0.2	0	0	0.997725	-22	1.01295	38.9291
	57.8	28.4	10.6	2.6	0	0.1	0	0	0.2	0.2	0.1	0	1.02462	8	1.03652	37.2415
	72.3	13.4	10.3	3	0	0.2	0.2	0	0.2	0.2	0.2	0	1.03056	19	1.0359	35.7234
	62.2	20.7	9.4	6.2	0	0	0.6	0	0.6	0.3	0	0	1.09681	80	1.09214	34.0022
	35.5	45.5	9.5	2.5	4.5	2.1	0.1	0	0.1	0.1	0.1	0	1.08746	120	1.0832	33.8555
	36.3	42.3	8.9	2.7	6.4	3	0.1	0	0.1	0.1	0.1	0	1.13485	200	1.11908	31.2418
	37.2	39.1	8.2	2.9	8.3	3.9	0.1	0	0.1	0.1	0.1	0	1.19409	300	1.16414	28.6198
	38	36.1	7.6	3	10.1	4.7	0.1	0.1	0.2	0.1	0	0	1.25333	400	1.20957	26.3764
	38.7	33.5	7.1	3.2	11.7	5.4	0.1	0.1	0.2	0	0	0	1.31257	500	1.25601	24.5068
	39.4	31	6.6	3.3	13.2	6.1	0.1	0.1	0.2	0	0	0	1.37181	600	1.30145	22.8421
	40	28.7	6.1	3.5	14.6	6.7	0.1	0.1	0.2	0	0	0	1.43104	700	1.34623	21.3871
	40.5	26.5	5.6	3.6	15.9	7.3	0.1	0.2	0.3	0	0	0	1.49028	800	1.38998	20.0717
	41.1	24.6	5.2	3.7	17	7.8	0.1	0.2	0.3	0	0	0	1.54952	900	1.43531	18.9582
	41.6	22.7	4.9	3.8	18.1	8.3	0.1	0.2	0.3	0	0	0	1.60876	1000	1.4818	17.9454
	42	21	4.5	3.9	19.2	8.8	0.1	0.2	0.3	0	0	0	1.668	1100	1.52574	17.0139
	42.5	19.4	4.2	4	20.1	9.2	0.1	0.2	0.3	0	0	0	1.72724	1200	1.57143	16.2022
	42.9	17.9	3.9	4.1	21	9.6	0.1	0.2	0.3	0	0	0	1.78648	1300	1.61658	15.4529
	43.2	16.5	3.6	4.2	21.9	10	0.1	0.2	0.3	0	0	0	1.84572	1400	1.66118	14.7594
	43.5	15.5	3.4	4.2	22.5	10.3	0.1	0.2	0.3	0	0	0	1.90496	1500	1.7082	14.1711

Acknowledgements

The work of Nils Krah was supported by funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 753370. The work of Ilaria Rinaldi was supported partially by a research fellowship "Visiting Professor Sapienza" from the University of Rome. This work was performed within the framework the project DIC20161236452 of the "Fondation pour la Recherche Médicale" (FRM), of the SIRIC LYriCAN Grant INCa INSERM DGOS 12563 and 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 ANR. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan X Pascal GPU used for this research.